deseq2_fun.R 7.75 KB
Newer Older
1
#' DESEQ2 Differential Analysis function.
Etienne Rifa's avatar
Etienne Rifa committed
2
3
#'
#'
Rifa Etienne's avatar
Rifa Etienne committed
4
#' @param data a phyloseq object (output from decontam or generate_phyloseq)
5
#' @param output Output directory
Rifa Etienne's avatar
Rifa Etienne committed
6
7
#' @param rank Taxonomy rank to merge features that have same taxonomy at a certain taxonomic rank (among rank_names(data), or 'ASV' for no glom)
#' @param column1 Column name of factor to test (among sample_variables(data))
8
9
10
11
#' @param verbose Verbose level. (1: quiet, 3: verbal)
#' @param comp Comma separated list of comparison to test. Comparisons are informed with a tilde (A~C,A~B,B~C). If empty, test all combination
#'
#'
Rifa Etienne's avatar
Rifa Etienne committed
12
#' @return Returns list with table of features and plots for each comparison. Exports CSV files with significant differentialy abondant ASV.
Etienne Rifa's avatar
Etienne Rifa committed
13
#'
Etienne Rifa's avatar
Etienne Rifa committed
14
#' @import futile.logger
Etienne Rifa's avatar
Etienne Rifa committed
15
#' @import phyloseq
Etienne Rifa's avatar
Etienne Rifa committed
16
17
18
19
#' @importFrom BiocGenerics estimateSizeFactors
#' @importFrom DESeq2 DESeq
#' @importFrom DESeq2 results
#' @importFrom DESeq2 resultsNames
Etienne Rifa's avatar
Etienne Rifa committed
20
#' @importFrom gridExtra grid.arrange
Etienne Rifa's avatar
Etienne Rifa committed
21
22
23
#' @importFrom ggpubr ggtexttable
#' @importFrom ggpubr ttheme
#' @importFrom ggpubr text_grob
Rifa Etienne's avatar
Rifa Etienne committed
24
#' @importFrom plotly ggplotly
Etienne Rifa's avatar
Etienne Rifa committed
25
26
#' @import ggplot2

Etienne Rifa's avatar
Etienne Rifa committed
27
28
29
30
31
32
#'
#' @export


# Decontam Function

33
34
deseq2_fun <- function(data = data, output = "./deseq/", column1 = "", verbose = 1,
                       rank = "Species", comp = ""){
Etienne Rifa's avatar
Etienne Rifa committed
35
36
37
38
39
40
41
42
43
44
45
46


  if(verbose == 3){
    invisible(flog.threshold(DEBUG))
  } else {
    invisible(flog.threshold(INFO))
  }


  if(!dir.exists(output)){
    dir.create(output, recursive = TRUE)
  }
Etienne Rifa's avatar
Etienne Rifa committed
47

Etienne Rifa's avatar
Etienne Rifa committed
48
  if(rank != 'ASV'){
Etienne Rifa's avatar
Etienne Rifa committed
49
    flog.info('Glom rank...')
Etienne Rifa's avatar
Etienne Rifa committed
50
51
52
53
54
55
56
    data.glom <- tax_glom(data, taxrank=rank)
  } else {
    data.glom <- data
  }

  # save.image("debug.rdata")
  # quit()
Etienne Rifa's avatar
Etienne Rifa committed
57
58
59
  flog.info('Defining comparison...')
  if(comp==""){
    fun <- paste('combinaisons <- utils::combn(na.omit(unique(sample_data(data.glom)$',column1,')),2) ',sep='')
Etienne Rifa's avatar
Etienne Rifa committed
60
    eval(parse(text=fun))
Etienne Rifa's avatar
Etienne Rifa committed
61
    print(combinaisons)
Etienne Rifa's avatar
Etienne Rifa committed
62
63
64
65
66
67
68
69
70
71
72
73
  }else{
    comp_list <- unlist(strsplit(comp,","))
    # combinaisons <- combn(comp_list,2)
    combinaisons <- matrix(, nrow = 2, ncol = length(comp_list))
    for (i in 1:length(comp_list)){
      tmp <- unlist(strsplit(comp_list[i],"\\~"))
      # cbind(combinaisons,tmp)
      combinaisons[1,i] <- tmp[1]
      combinaisons[2,i] <- tmp[2]
    }
  }

Etienne Rifa's avatar
Etienne Rifa committed
74
  flog.info('Done...')
Etienne Rifa's avatar
Etienne Rifa committed
75

Rifa Etienne's avatar
Rifa Etienne committed
76
  outF = list()
Etienne Rifa's avatar
Etienne Rifa committed
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
  '%!in%' <- function(x,y)!('%in%'(x,y))
  for (col in (1:ncol(combinaisons))){
    fun <- paste('tmp <- sample_data(data)$',column1,sep='')
    eval(parse(text=fun))
    if((combinaisons[1,col] %!in% tmp) || combinaisons[2,col] %!in% tmp){
      flog.warn(paste(combinaisons[1,col],' not in sample_data. Next;'),sep='')
      next
    }
    flog.info(paste('Combinaison ',combinaisons[1,col], ' ' , combinaisons[2,col],sep=''))
    fun <- paste('nb_cond1 <- nsamples(subset_samples(data.glom, ',column1, ' == "',combinaisons[1,col],'"))',sep='')
    eval(parse(text=fun))
    print(nb_cond1)
    fun <- paste('nb_cond2 <- nsamples(subset_samples(data.glom, ',column1, ' == "',combinaisons[2,col],'"))',sep='')
    eval(parse(text=fun))
    print(nb_cond2)

    if(nb_cond1 < 2){
      flog.warn(paste('Condition: ',combinaisons[1,col],' only have ',nb_cond1,' sample(s). Skip...',sep=''))
      tab0 = data.frame(matrix(ncol = 14, nrow = 0))
      write.table(tab0, file = paste(output,'/signtab_',column1,'_',paste(combinaisons[,col],collapse="_vs_"),'.csv',sep=''),quote=FALSE,sep="\t", row.names=FALSE)
      next
    } else if (nb_cond2 < 2){
      flog.warn(paste('Condition: ',combinaisons[2,col],' only have ',nb_cond2,' sample(s). Skip...',sep=''))
      tab0 = data.frame(matrix(ncol = 14, nrow = 0))
      write.table(tab0, file = paste(output,'/signtab_',column1,'_',paste(combinaisons[,col],collapse="_vs_"),'.csv',sep=''),quote=FALSE,sep="\t", row.names=FALSE)
      next
    } else if (nb_cond1 < 2 & nb_cond2 < 2){
      tab0 = data.frame(matrix(ncol = 14, nrow = 0))
      write.table(tab0, file = paste(output,'/signtab_',column1,'_',paste(combinaisons[,col],collapse="_vs_"),'.csv',sep=''),quote=FALSE,sep="\t", row.names=FALSE)
      flog.warn(paste('Condition: ',combinaisons[1,col],' only have ',nb_cond1,' sample(s). Skip...',sep=''))
      flog.warn(paste('Condition: ',combinaisons[2,col],' only have ',nb_cond2,' sample(s). Skip...',sep=''))
      next
    } else {
      flog.warn(paste('Condition: ',combinaisons[1,col],' contains ',nb_cond1,' samples.',sep=''))
      flog.warn(paste('Condition: ',combinaisons[2,col],' contains ',nb_cond2,' samples.',sep=''))
    }



    fun <- paste('tmp <- subset_samples(data.glom, ',column1, ' %in% c("',combinaisons[1,col],'","',combinaisons[2,col],'"))',sep='')
    eval(parse(text=fun))

    # if(nsamples(tmp)<4){
    # 	flog.warn('Less than 4 samples')
    # 	next;
    # }

    tmp <- prune_taxa(taxa_sums(tmp) >= 1, tmp)
Etienne Rifa's avatar
Etienne Rifa committed
125
    # return(tmp)
Etienne Rifa's avatar
Etienne Rifa committed
126
127

    flog.info('DESeq2...')
Etienne Rifa's avatar
Etienne Rifa committed
128
    fun <- paste('dseq <- phyloseq_to_deseq2(tmp, ~ ',column1,')',sep='')
Etienne Rifa's avatar
Etienne Rifa committed
129
130
131
132
133
    eval(parse(text = fun))

    gm_mean = function(x, na.rm=TRUE){
      exp(sum(log(x[x > 0]), na.rm=na.rm) / length(x))
    }
Etienne Rifa's avatar
Etienne Rifa committed
134
135
    geoMeans = apply(DESeq2::counts(dseq), 1, gm_mean)
    dseq2 = estimateSizeFactors(dseq, geoMeans = geoMeans)
Etienne Rifa's avatar
Etienne Rifa committed
136

Etienne Rifa's avatar
Etienne Rifa committed
137
138
    # flog.info('DESeq2...')
    # print(dseq2)
Etienne Rifa's avatar
Etienne Rifa committed
139

Etienne Rifa's avatar
Etienne Rifa committed
140
141
    dseq3 = DESeq2::DESeq(dseq2, test="Wald", fitType="parametric")
    flog.debug(show(dseq3))
Etienne Rifa's avatar
Etienne Rifa committed
142

Etienne Rifa's avatar
Etienne Rifa committed
143
    res = results(dseq3, cooksCutoff = FALSE, contrast = c(column1,combinaisons[1,col] , combinaisons[2,col]))
Etienne Rifa's avatar
Etienne Rifa committed
144
145
    flog.debug(show(res))

Rifa Etienne's avatar
Rifa Etienne committed
146
147
    # outlist1 <- list()

Etienne Rifa's avatar
Etienne Rifa committed
148
149
    flog.info('Plot...')
    alpha = 0.05
Etienne Rifa's avatar
Etienne Rifa committed
150
    if(length(which(res$padj < alpha)) > 0){
Etienne Rifa's avatar
Etienne Rifa committed
151
152
153
154
      sigtab = res[which(res$padj < alpha), ]
      # sigtab = res
      tabOUT = cbind(taxon_id = row.names(res),as(res, "data.frame"), as(tax_table(data)[rownames(res), ], "matrix"))
      colnames(res)[1]=resultsNames(dseq3)[2]
Rifa Etienne's avatar
Rifa Etienne committed
155
      # save.image("debug.rdata")
Etienne Rifa's avatar
Etienne Rifa committed
156
      write.table(tabOUT, file = paste(output,'/signtab_',column1,'_',paste(combinaisons[,col],collapse="_vs_"),'.csv',sep=''),quote=FALSE,sep="\t", row.names=FALSE)
Etienne Rifa's avatar
Etienne Rifa committed
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171

      if(rank != 'ASV'){
        fun <- paste('x = tapply(sigtab$log2FoldChange, sigtab$',rank,', function(x) max(x))',sep='')
        fun2 <- paste('sigtab$',rank,' = factor(as.character(sigtab$',rank,'), levels=names(x))')
      } else {
        fun <- paste('x = tapply(sigtab$log2FoldChange, row.names(sigtab), function(x) max(x))',sep='')
        fun2 <- paste('sigtab$',rank,' = factor(as.character(row.names(sigtab)), levels=names(x))')
      }
      eval(parse(text=fun))
      x = sort(x, TRUE)
      eval(parse(text=fun2))

      fun <- paste('p <- ggplot(sigtab,
			aes(x=',rank,',
			y=log2FoldChange,
Sebastien Theil's avatar
Sebastien Theil committed
172
			color=',rank,')) + geom_point(size=6) + theme(axis.text.x = element_text(angle = -90, hjust = 0, vjust=0.5), legend.position = "none")',sep='')
Etienne Rifa's avatar
Etienne Rifa committed
173
174
      eval(parse(text=fun))

Rifa Etienne's avatar
Rifa Etienne committed
175
176
177
      Ftable = sigtab[,c("baseMean","log2FoldChange","stat","pvalue","padj",rank)]
      ggtable <- ggtexttable(Ftable,theme = ttheme("mOrange"),rows=NULL)

Etienne Rifa's avatar
Etienne Rifa committed
178
      pdf(file=paste(output,'/deseq2_',column1,"_", paste(combinaisons[,col],collapse="_vs_"), '.pdf',sep=''),width=15,height=16)
Etienne Rifa's avatar
Etienne Rifa committed
179
      grid.arrange(p,ggtable,top=text_grob(paste('Combination ',combinaisons[1,col], ' VS ' , combinaisons[2,col],sep=''), size=20))
Rifa Etienne's avatar
Rifa Etienne committed
180
181
182
183
184
185
      invisible(dev.off())


      outF[[paste(combinaisons[,col],collapse="_vs_")]] = list(plot = p, table = Ftable)


Etienne Rifa's avatar
Etienne Rifa committed
186
187
188
189
190
    } else{
      flog.info(paste('No significant results for comparison ',combinaisons[1,col], ' ' , combinaisons[2,col], sep=''))
      tab0 = data.frame(matrix(ncol = 14, nrow = 0))
      # print(length(c(resultsNames(deseq)[2], colnames(res), colnames(tax_table(data)))))
      # print(c(resultsNames(deseq)[2], colnames(res), colnames(tax_table(data))))
Etienne Rifa's avatar
Etienne Rifa committed
191
      colnames(tab0) <- c(resultsNames(dseq3)[2], colnames(res), colnames(tax_table(data)))
Etienne Rifa's avatar
Etienne Rifa committed
192
193
194
195
      write.table(tab0, file = paste(output,'/signtab_',column1,'_',paste(combinaisons[,col],collapse="_vs_"),'.csv',sep=''),quote=FALSE,sep="\t", row.names=FALSE)

    }
  }
Sebastien Theil's avatar
Sebastien Theil committed
196
  return(outF)
Etienne Rifa's avatar
Etienne Rifa committed
197
  flog.info('Done...')
Etienne Rifa's avatar
Etienne Rifa committed
198
}