Commit a9fcc852 authored by Etienne Rifa's avatar Etienne Rifa
Browse files

deseq2 fix

beta_light allows to plot other ordination axis
close #13
parent fbb2a3cc
Pipeline #40813 passed with stage
in 10 seconds
...@@ -134,25 +134,26 @@ deseq2_fun <- function(data = data, output = "./deseq/", column1 = "", verbose = ...@@ -134,25 +134,26 @@ deseq2_fun <- function(data = data, output = "./deseq/", column1 = "", verbose =
geoMeans = apply(DESeq2::counts(dseq), 1, gm_mean) geoMeans = apply(DESeq2::counts(dseq), 1, gm_mean)
dseq2 = estimateSizeFactors(dseq, geoMeans = geoMeans) dseq2 = estimateSizeFactors(dseq, geoMeans = geoMeans)
flog.info('DESeq2...') # flog.info('DESeq2...')
print(dseq2) # print(dseq2)
dseq3 = DESeq2::DESeq(dseq2, test="Wald", fitType="parametric") dseq3 = DESeq2::DESeq(dseq2, test="Wald", fitType="parametric")
flog.debug(show(dseq3)) flog.debug(show(dseq3))
res = results(dseq3, cooksCutoff = FALSE) res = results(dseq3, cooksCutoff = FALSE, contrast = c(column1,combinaisons[1,col] , combinaisons[2,col]))
flog.debug(show(res)) flog.debug(show(res))
alpha = 0.05
# outlist1 <- list() # outlist1 <- list()
flog.info('Plot...')
alpha = 0.05
if(length(which(res$padj < alpha)) > 0){ if(length(which(res$padj < alpha)) > 0){
# sigtab = res[which(res$padj < alpha), ] sigtab = res[which(res$padj < alpha), ]
sigtab = res # sigtab = res
sigtab = cbind(row.names(sigtab),as(sigtab, "data.frame"), as(tax_table(data)[rownames(sigtab), ], "matrix")) tabOUT = cbind(taxon_id = row.names(res),as(res, "data.frame"), as(tax_table(data)[rownames(res), ], "matrix"))
colnames(sigtab)[1]=resultsNames(dseq3)[2] colnames(res)[1]=resultsNames(dseq3)[2]
# save.image("debug.rdata") # save.image("debug.rdata")
write.table(sigtab, file = paste(output,'/signtab_',column1,'_',paste(combinaisons[,col],collapse="_vs_"),'.csv',sep=''),quote=FALSE,sep="\t", row.names=FALSE) write.table(tabOUT, file = paste(output,'/signtab_',column1,'_',paste(combinaisons[,col],collapse="_vs_"),'.csv',sep=''),quote=FALSE,sep="\t", row.names=FALSE)
if(rank != 'ASV'){ if(rank != 'ASV'){
fun <- paste('x = tapply(sigtab$log2FoldChange, sigtab$',rank,', function(x) max(x))',sep='') fun <- paste('x = tapply(sigtab$log2FoldChange, sigtab$',rank,', function(x) max(x))',sep='')
...@@ -174,7 +175,7 @@ deseq2_fun <- function(data = data, output = "./deseq/", column1 = "", verbose = ...@@ -174,7 +175,7 @@ deseq2_fun <- function(data = data, output = "./deseq/", column1 = "", verbose =
Ftable = sigtab[,c("baseMean","log2FoldChange","stat","pvalue","padj",rank)] Ftable = sigtab[,c("baseMean","log2FoldChange","stat","pvalue","padj",rank)]
ggtable <- ggtexttable(Ftable,theme = ttheme("mOrange"),rows=NULL) ggtable <- ggtexttable(Ftable,theme = ttheme("mOrange"),rows=NULL)
pdf(file=paste(output,'/deseq2_',column1,'.pdf',sep=''),width=15,height=16) pdf(file=paste(output,'/deseq2_',column1,"_", paste(combinaisons[,col],collapse="_vs_"), '.pdf',sep=''),width=15,height=16)
grid.arrange(p,ggtable,top=text_grob(paste('Combination ',combinaisons[1,col], ' VS ' , combinaisons[2,col],sep=''), size=20)) grid.arrange(p,ggtable,top=text_grob(paste('Combination ',combinaisons[1,col], ' VS ' , combinaisons[2,col],sep=''), size=20))
invisible(dev.off()) invisible(dev.off())
......
...@@ -10,6 +10,7 @@ ...@@ -10,6 +10,7 @@
#' @param ord0 Currently supported method options are: c("DCA", "CCA", "RDA", "CAP", "DPCoA", "NMDS", "MDS", "PCoA") #' @param ord0 Currently supported method options are: c("DCA", "CCA", "RDA", "CAP", "DPCoA", "NMDS", "MDS", "PCoA")
#' @param output The output file directory. #' @param output The output file directory.
#' @param tests Whether to compute tests or not (TRUE/FALSE) #' @param tests Whether to compute tests or not (TRUE/FALSE)
#' @param axes Axes to plot (c(1,2))
#' @param verbose Verbose level. (1: quiet, 2: print infos, 3: print infos + debug) #' @param verbose Verbose level. (1: quiet, 2: print infos, 3: print infos + debug)
#' #'
#' @return Return specific plots and tests in list and output them in the output directory. #' @return Return specific plots and tests in list and output them in the output directory.
...@@ -25,7 +26,7 @@ ...@@ -25,7 +26,7 @@
# Decontam Function # Decontam Function
diversity_beta_light <- function(psobj, rank = "ASV", col = NULL, cov = NULL, dist0 = "bray", ord0 = "MDS", output="./plot_div_beta/", tests = TRUE, verbose = 2) { diversity_beta_light <- function(psobj, rank = "ASV", col = NULL, cov = NULL, dist0 = "bray", ord0 = "MDS", output="./plot_div_beta/", axes = c(1,2), tests = TRUE, verbose = 2) {
if(verbose == 3){ if(verbose == 3){
invisible(flog.threshold(DEBUG)) invisible(flog.threshold(DEBUG))
...@@ -69,10 +70,10 @@ diversity_beta_light <- function(psobj, rank = "ASV", col = NULL, cov = NULL, di ...@@ -69,10 +70,10 @@ diversity_beta_light <- function(psobj, rank = "ASV", col = NULL, cov = NULL, di
flog.info('Plot ...') flog.info('Plot ...')
resBeta <- list() resBeta <- list()
if(!is.null(cov)){ if(!is.null(cov)){
p1 <- plot_samples(data_rank, ordinate(data_rank, ord0, dist0), color = col, shape = cov1[length(cov1)] ) + p1 <- plot_samples(data_rank, ordinate(data_rank, ord0, dist0), color = col, shape = cov1[length(cov1)], axes = axes ) +
theme_bw() + ggtitle(glue::glue("{ord0} + {dist0}")) + stat_ellipse() + scale_shape_manual(values = 0:10) theme_bw() + ggtitle(glue::glue("{ord0} + {dist0}")) + stat_ellipse() + scale_shape_manual(values = 0:10)
}else{ }else{
p1 <- plot_samples(data_rank, ordinate(data_rank, ord0, dist0), color = col) + p1 <- plot_samples(data_rank, ordinate(data_rank, ord0, dist0), color = col, axes = axes ) +
theme_bw() + ggtitle(glue::glue("{ord0} + {dist0}")) + stat_ellipse() theme_bw() + ggtitle(glue::glue("{ord0} + {dist0}")) + stat_ellipse()
} }
# plot(p1) # plot(p1)
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment