internal.R 51.3 KB
Newer Older
stevenn's avatar
stevenn committed
1

2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

## Modified version of expand.grid
expand.grid2.list <- function(listInput)
{
  n = length(listInput)
  if(is.list(listInput) && n>1)
  {
    l1 = listInput[[1]]
    l2 = listInput[[2]]
    res = c()
    
    for(i in l1){
      for(j in l2){ 
        res = rbind(res,paste(i,j,sep = "-"))
      }
    }
    listInput[[1]] = res
    listInput = listInput[-2]
    if(length(listInput)>1 && is.list(listInput)) res = expand.grid2.list(listInput)
  }
  else res = listInput
  return(res)
}


Stevenn Volant's avatar
Stevenn Volant committed
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
## Function for the rdp format
getval <- function(annotation_rdp, interest, threshold_annot){
  annotation_rdp = unlist(strsplit(annotation_rdp,"\t"))
  annotation = c(annotation_rdp[1])
  for(level in interest){
    idlevel=which(annotation_rdp == level)
    if(length(idlevel)>0){
      if(as.numeric(annotation_rdp[idlevel+1]) >= threshold_annot){
        annotation = c(annotation, gsub("\"", "", annotation_rdp[idlevel-1]))
      }
      else annotation = c(annotation, "NA")
    }
    else annotation = c(annotation, "NA")  
  }
  return(annotation)
}

## Read rdp file
read_rdp <- function(filename, threshold_annot)
{
  interest=c("phylum", "class", "order", "family", "genus")
  conn <- file(filename,open="r")
  linn <-readLines(conn)
  tab=t(sapply(1:length(linn), function(i) getval(linn[i], interest, threshold_annot)))
  close(conn)
  
  if(!TRUE%in%duplicated(tab[,1])) rownames(tab)=tab[,1];tab=tab[,-1]
  colnames(tab) = c("Phylum","Class","Order","Family","Genus")
  
  return(tab)
}






CheckCountsTable <- function(counts)
  {
    Error = NULL
    Warning = NULL
    numTest = FALSE%in%sapply(counts,is.numeric)
    if(ncol(counts)<=1){Error = "The number of columns of the counts table must be at least 2" }
    if(nrow(counts)<=1){Error = "The number of rows of the counts table must be at least 2" }
    if(numTest){Error = "The counts table must contain only numeric values" }
    if(!numTest)
    {
      if(0%in%colSums(counts)){Error = "At least one of the column of the counts table is 0" }
      if(min(counts)<0){Error = "The counts table must contain only positive values" }
    }
    if(TRUE%in%sapply(counts,is.na)){Warning = "NA values are considered as 0 is the counts table"; counts[sapply(counts,is.na)]=0}
    
  
    return(list(Error=Error,Warning=Warning,counts=counts))
  }
  
  CheckTaxoTable <- function(taxo)
stevenn's avatar
stevenn committed
84
  {
Stevenn Volant's avatar
Stevenn Volant committed
85
86
    Error = NULL
    Warning = NULL
svolant's avatar
svolant committed
87
    if(ncol(taxo)<1){Error = "The number of columns of the taxonomy table must be at least 1" }
Stevenn Volant's avatar
Stevenn Volant committed
88
89
90
91
92
93
94
95
96
    if(nrow(taxo)<=1){Error = "The number of rows if the taxonomy table must be at least 2" }
    if(TRUE%in%is.numeric(taxo)){Error = "The taxonomy table must contain only character" }

    for(i in 1:ncol(taxo))
    {
      level = levels(taxo[,i])
      nb = length(level)
      if(nb==1 && level=="NA"){ Error = "At least one column contains only NA"}
    }
stevenn's avatar
stevenn committed
97
    
Stevenn Volant's avatar
Stevenn Volant committed
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
    return(list(Error=Error,Warning=Warning))
  }
  
  PercentAnnot <- function(counts,taxo)
  {
    Error = NULL  
    tmp = table(rownames(counts)%in%rownames(taxo))
    Percent = tmp["TRUE"]/sum(tmp)
    if(is.na(Percent)) Percent = 0
    if(Percent==0){Error = "Counts table and annotation do not matched" }
       
    return(list(Error=Error,Percent=Percent))
  }
  
  
  GetDataFromBIOM <-function(dataBIOM)
  {
    ## Counts table
stevenn's avatar
stevenn committed
116
    counts = biom_data(dataBIOM)
stevenn's avatar
stevenn committed
117
118
    counts = as.matrix(counts)
    counts = as.data.frame(counts)
Stevenn Volant's avatar
Stevenn Volant committed
119
120
121
122
    CheckCounts = CheckCountsTable(counts)
    counts = CheckCounts$counts
    
    ## Taxonomy table
stevenn's avatar
stevenn committed
123
    taxo = as.data.frame(observation_metadata(dataBIOM))
Stevenn Volant's avatar
Stevenn Volant committed
124
125
126
127
128
129
    CheckTaxo = CheckTaxoTable(taxo)
    
    ## Pourcentage of annotation
    tmp = PercentAnnot(counts,taxo)
    
    return(list(counts=counts,taxo=taxo,CheckCounts=CheckCounts,CheckTaxo=CheckTaxo,Percent=tmp$Percent,CheckPercent=tmp$Error))
stevenn's avatar
stevenn committed
130
  }
stevenn's avatar
stevenn committed
131
132
  
  
stevenn's avatar
stevenn committed
133
134
135
  GetDataFromCT <-function(dataC,dataT)
  {
    
Stevenn Volant's avatar
Stevenn Volant committed
136
    ## Counts table
stevenn's avatar
stevenn committed
137
    counts = dataC
Stevenn Volant's avatar
Stevenn Volant committed
138
139
140
141
142
143
144
145
146
147
148
    CheckCounts = CheckCountsTable(counts)
    counts = CheckCounts$counts
    
    ## Taxonomy table
    taxo = as.data.frame(dataT)
    CheckTaxo = CheckTaxoTable(taxo)
    
    ## Pourcentage of annotation
    tmp = PercentAnnot(counts,taxo)
    
    return(list(counts=counts,taxo=taxo,CheckCounts=CheckCounts,CheckTaxo=CheckTaxo,Percent=tmp$Percent,CheckPercent=tmp$Error))
stevenn's avatar
stevenn committed
149
150
151
  }
  
  GetInteraction2 <- function(target)
stevenn's avatar
stevenn committed
152
  { 
stevenn's avatar
stevenn committed
153
154
155
156
    res=c()
    namesTarget = colnames(target)[2:ncol(target)]
    k=1
    for(i in 1:(length(namesTarget)-1))
stevenn's avatar
stevenn committed
157
    { 
stevenn's avatar
stevenn committed
158
159
160
161
162
      for(j in (i+1):length(namesTarget))
      { 
        res[k] = paste(namesTarget[i],":",namesTarget[j],sep="")
        k = k+1
      }
stevenn's avatar
stevenn committed
163
    }
stevenn's avatar
stevenn committed
164
165
    
    return(res)
stevenn's avatar
stevenn committed
166
167
168
169
  }
  


Amine  GHOZLANE's avatar
Amine GHOZLANE committed
170
  ## Print the contrasts
stevenn's avatar
stevenn committed
171
172
173
174
175
176
177
178
179
180
  PrintContrasts <- function (coefs, contrasts,contnames) 
  {
    contrasts = as.matrix(contrasts)
    out <-""
    
    for (i in 1:ncol(contrasts)) 
    {
      contrast <- contrasts[,i]
      contrast <- paste(ifelse(contrast > 0, "+ ", ""), contrast, sep = "")
      contrast <- gsub("( 1)|(1)", "", contrast)
181
      out = paste(out,paste("<b>",contnames[i], "</b> <br/>", paste(contrast[contrast != 0], coefs[contrast != 0], collapse = " ", sep = " ")),"<br/>")
stevenn's avatar
stevenn committed
182
183
184
185
    }
    return(out)
    
  }
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
186
187

  
stevenn's avatar
stevenn committed
188
189
  
  ## Get the counts for the selected taxonomy
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
190
  GetCountsMerge <- function(input,dataInput,taxoSelect,target,design)
stevenn's avatar
stevenn committed
191
  {
192
    ## Init
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
193
194
    counts= NULL
    CheckTarget = FALSE
Stevenn Volant's avatar
Stevenn Volant committed
195
196
    CT_noNorm = NULL
    normFactors = NULL
197
198
    FeatureSize = NULL

Amine  GHOZLANE's avatar
Amine GHOZLANE committed
199
    ## Counts and taxo tables
stevenn's avatar
stevenn committed
200
201
    CT = dataInput$counts
    taxo = dataInput$taxo
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
202
203
204
205
        
    ## Select cols in the target
    labels = target[,1]
    ind = which(colnames(CT)%in%labels)
stevenn's avatar
stevenn committed
206
    
207
208
209
    ## Get the feature size for the normalisation
    Size_indcol = which(toupper(colnames(CT))%in%"SIZE")
    if(length(Size_indcol)==1) FeatureSize = CT[,Size_indcol]
stevenn's avatar
stevenn committed
210
    
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
211
212
213
214
215
    if(length(ind)==length(labels))
    { 
      CT = CT[,ind]
      
      ## Order CT according to the target
216
      CT = OrderCounts(counts=CT,labels=labels)$CountsOrder
Stevenn Volant's avatar
Stevenn Volant committed
217
      CT_noNorm = CT
218
      RowProd = sum(apply(CT_noNorm,1,prod))
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
219
220
221
      
      ## Counts normalisation
      dds <- DESeqDataSetFromMatrix(countData=CT, colData=target, design=design)
222
223
224
225
226
      
      ## Normalisation with or without 0
      if(input$AccountForNA || RowProd==0) dds = estimateSizeFactors(dds,locfunc=eval(as.name(input$locfunc)),geoMeans=GeoMeansCT(CT))
      if(!input$AccountForNA && RowProd!=0) dds = estimateSizeFactors(dds,locfunc=eval(as.name(input$locfunc)))
       
Stevenn Volant's avatar
Stevenn Volant committed
227
      normFactors = sizeFactors(dds)
228
        
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
229
230
231
      ordOTU = order(rownames(taxo))
      indOTU_annot = which(rownames(CT)%in%rownames(taxo))
      counts_annot = CT[indOTU_annot[ordOTU],]
232

Amine  GHOZLANE's avatar
Amine GHOZLANE committed
233
234
      if(taxoSelect=="OTU") counts = counts_annot
      else{
Stevenn Volant's avatar
Stevenn Volant committed
235
236
237
        taxoS = taxo[ordOTU,taxoSelect]
        counts = aggregate(counts_annot,by=list(Taxonomy = taxoS),sum)
        rownames(counts)=counts[,1];counts=counts[,-1]
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
238
239
240
      }
      
      ## Ordering the counts table according to the target labels 
241
242
243
      tmpOrder = OrderCounts(counts,normFactors,labels)
      counts = tmpOrder$CountsOrder
      normFactors = tmpOrder$normFactorsOrder
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
244
245
      CheckTarget = TRUE
    }
Stevenn Volant's avatar
Stevenn Volant committed
246
    return(list(counts=counts,CheckTarget=CheckTarget,normFactors=normFactors, CT_noNorm=CT_noNorm))
stevenn's avatar
stevenn committed
247
  }
248
249
250
251
252
253
254
255
256
257
258
  
  ## Get the geometric mean of the counts (0 are replaced by NA values)
  GeoMeansCT <- function(CT)
  {
    CT=as.matrix(CT)
    CT[which(CT<1)]=NA
    gm = apply(CT,1,geometric.mean,na.rm=TRUE)
    return(gm)
  }
  
  
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
259
  ## Order the counts 
260
  OrderCounts <- function(counts,normFactors=NULL,labels)
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
261
262
263
  {
    n = length(labels)
    CountsOrder = counts
264
    normFactorsOrder = normFactors
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
265
266
267
268
    for(i in 1:n)
    {
      ind = which(labels[i]==colnames(counts))
      CountsOrder[,i] = counts[,ind]
269
      if(!is.null(normFactors)) normFactorsOrder[i] = normFactors[ind]
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
270
271
    }
    colnames(CountsOrder) = labels
272
    return(list(CountsOrder=CountsOrder,normFactorsOrder = normFactorsOrder))
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
273
274
275
  }
  
  
stevenn's avatar
stevenn committed
276
  ## Get the dds object of DESeq2
Stevenn Volant's avatar
Stevenn Volant committed
277
  Get_dds_object <- function(input,counts,target,design,normFactorsOTU,CT_noNorm)
stevenn's avatar
stevenn committed
278
279
  {
    dds <- DESeqDataSetFromMatrix(countData=counts, colData=target, design=design)
280
    sizeFactors(dds) = normFactorsOTU
stevenn's avatar
stevenn committed
281
    dds <- estimateDispersions(dds, fitType=input$fitType)
Stevenn Volant's avatar
Stevenn Volant committed
282
    dds <- nbinomWaldTest(dds)
283
284
    countsNorm = counts(dds, normalized = TRUE)
    return(list(dds = dds,raw_counts=counts,countsNorm=countsNorm,target=target,design=design,normFactors = normFactorsOTU,CT_noNorm=CT_noNorm))
stevenn's avatar
stevenn committed
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
  }

  ## Get the design according to the input
  GetDesign <- function(input)
  {
    InterVar = input$InterestVar
    Interaction = input$Interaction2
    alltmp = c(InterVar,Interaction)
    design = as.formula(paste("~", paste0(alltmp, collapse= "+")))
    return(design)
  }
  


  ## Diagnostic Plots
  Plot_diag <- function(input,resDiff)
  {
Stevenn Volant's avatar
Stevenn Volant committed
302
    
stevenn's avatar
stevenn committed
303
304
    VarInt = input$VarInt
    dds = resDiff$dds
305
306
    counts = resDiff$raw_counts
    if(input$CountsType=="norm") counts = resDiff$countsNorm
stevenn's avatar
stevenn committed
307
    target = resDiff$target
308
    normFactors = resDiff$normFactors
Stevenn Volant's avatar
Stevenn Volant committed
309
    CT_noNorm = resDiff$CT_noNorm
stevenn's avatar
stevenn committed
310
    group = as.data.frame(target[,VarInt])
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
311
    rownames(group) = rownames(target)
Stevenn Volant's avatar
Stevenn Volant committed
312
    res = NULL
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
313
    
Stevenn Volant's avatar
Stevenn Volant committed
314
315
    if(ncol(group)>0 && nrow(counts)>0)
    { 
Stevenn Volant's avatar
Stevenn Volant committed
316
317
      colors = rep(c("#1f77b4","#aec7e8","#ff7f0e","#ffbb78", "#2ca02c","#98df8a","#d62728","#ff9896","#9467bd","#c5b0d5","#8c564b",
                     "#c49c94","#e377c2","#f7b6d2","#7f7f7f", "#c7c7c7","#bcbd22","#dbdb8d","#17becf","#9edae5"),ceiling(nrow(target)/20))
Stevenn Volant's avatar
Stevenn Volant committed
318
319
320
321
      
      if(input$DiagPlot=="barplotTot") res = barplotTot(input,counts,group = group, col=colors)
      if(input$DiagPlot=="barplotNul") res = barPlotNul(input,counts, group = group, col=colors)
      if(input$DiagPlot=="densityPlot") res = densityPlotTot(input,counts, group = group, col=colors)
322
      if(input$DiagPlot=="DispPlot") res = plotDispEsts(dds)
Stevenn Volant's avatar
Stevenn Volant committed
323
      if(input$DiagPlot=="MajTax") res = majTaxPlot(input,counts, group = group, col=colors)
324
      if(input$DiagPlot=="SfactorsVStot") res = diagSFactors(input,normFactors,resDiff$raw_counts) 
Stevenn Volant's avatar
Stevenn Volant committed
325
      if(input$DiagPlot=="pcaPlot") res = PCAPlot_meta(input,dds, group,  type.trans = input$TransType, col = colors)
Stevenn Volant's avatar
Stevenn Volant committed
326
      if(input$DiagPlot=="pcoaPlot") res = PCoAPlot_meta(input,dds, group, col = colors) 
327
      if(input$DiagPlot=="clustPlot") res = HCPlot(input,dds,group,type.trans=input$TransType,counts,col=colors)
Stevenn Volant's avatar
Stevenn Volant committed
328
    }
329
    
Stevenn Volant's avatar
Stevenn Volant committed
330
    return(res)
stevenn's avatar
stevenn committed
331
332
  }

Amine  GHOZLANE's avatar
Amine GHOZLANE committed
333
  
334
  HCPlot <- function (input,dds,group,type.trans,counts,col = c("lightblue", "orange", "MediumVioletRed", "SpringGreen")) 
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
335
336
  {
    
Stevenn Volant's avatar
Stevenn Volant committed
337
338
    res = NULL
    
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
339
    ## Get the counts
340
341
    if (input$DistClust == "euclidean" && type.trans == "VST") counts <- assay(varianceStabilizingTransformation(dds))
    if (input$DistClust == "euclidean" && type.trans == "rlog") counts <- assay(rlogTransformation(dds))
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
342
343
344
345
346
347
    
    ## Get the group of leaf
    group = apply(group,1,paste, collapse = "-")    
    nb = length(unique((group)))
    
    ## Get the dendrogram
Stevenn Volant's avatar
Stevenn Volant committed
348
349
350
351
    if(input$DistClust!="sere") dist = vegdist(t(counts), method = input$DistClust)
    if(input$DistClust=="sere") dist = as.dist(SEREcoef(counts))
    hc <- hclust(dist, method = "ward.D")
    
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
352
353
354
    dend = as.dendrogram(hc)
    
    ## Get the type of dendrogram
Stevenn Volant's avatar
Stevenn Volant committed
355
    type <- input$typeHculst
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
356
357
    
    dend <- set(dend, "labels_cex", input$cexLabelDiag)
358
    if(input$colorHC) labels_colors(dend)<-col[as.integer(as.factor(group))][order.dendrogram(dend)]
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
359
360
    if(type=="hori") 
    { 
361
      par(cex=input$cexTitleDiag,mar=c(6,6,4,5))
Stevenn Volant's avatar
Stevenn Volant committed
362
      res = plot(dend, main = "Cluster dendrogram",xlab = paste(input$DistClust,"distance, Ward criterion",sep=" "),cex=input$cexLabelDiag)
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
363
364
365
    }  
    if(type!="hori")
    {
366
      par(cex=input$cexTitleDiag,mar=c(6,6,4,5))
Stevenn Volant's avatar
Stevenn Volant committed
367
      res = circlize_dendrogram(dend, labels_track_height = 0.2, dend_track_height = .3, main = "Cluster dendrogram",xlab = paste(input$DistClust,"distance, Ward criterion",sep=" "))
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
368
    }
Stevenn Volant's avatar
Stevenn Volant committed
369
    return(res)
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
370
371
372
373
374
375
376
377
378
379
  }
  
  
  ## Color for the horizontal dendro
  colLabdendo <- function(n,group) {
    
    group = apply(group,1,paste, collapse = "-")
    
    nb = length(unique((group)))
    namesGrp = names(group)
stevenn's avatar
stevenn committed
380

Amine  GHOZLANE's avatar
Amine GHOZLANE committed
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
    if (is.leaf(n)) {
      a <- attributes(n)
      labCol <- rainbow(nb)[as.integer(as.factor(group))[which(namesGrp == a$label)]]
      attr(n, "nodePar") <- c(a$nodePar, lab.col = labCol)
    }
    return(n)
  }
  
  ## Diagnostic Plots Eigen value
  Plot_diag_Eigen <- function(input,resDiff)
  {
    colors = c("dodgerblue","firebrick1","MediumVioletRed","SpringGreen")
    VarInt = input$VarInt
    dds = resDiff$dds
    counts = resDiff$counts
    target = resDiff$target
    group = as.data.frame(target[,VarInt])
    
    ## If more than 4 levels for one factor
    maxFact =max(sapply(group,FUN=function(x) length(levels(x))))
    if(maxFact>=4) colors = rainbow(maxFact) 
    
    PCAPlot_meta(input,dds, group,  type.trans = input$TransType, col = colors, plot = "eigen") 
  }
  
  Plot_diag_pcoaEigen = function(input,resDiff)
  {
    colors = c("SpringGreen","dodgerblue","black","firebrick1")
    VarInt = input$VarInt
    dds = resDiff$dds
    target = resDiff$target
    group = as.data.frame(target[,VarInt])
    rownames(group) = rownames(target)
    PCoAPlot_meta(input,dds, group, col = colors, plot = "eigen") 
  }
  
  
  
stevenn's avatar
stevenn committed
419
420

  ## barplot total
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
421
  barplotTot <- function(input,counts, group, cex.names = 1, col = c("lightblue","orange", "MediumVioletRed", "SpringGreen")) 
stevenn's avatar
stevenn committed
422
  {
Stevenn Volant's avatar
Stevenn Volant committed
423

stevenn's avatar
stevenn committed
424
    ncol1 <- ncol(group) == 1
Stevenn Volant's avatar
Stevenn Volant committed
425
    par(cex=input$cexTitleDiag,mar=c(6,6,4,5))
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
426
    barplot(colSums(counts), cex.names = cex.names, main = "Total mapped read count per sample", ylab = "Total mapped read count", 
stevenn's avatar
stevenn committed
427
428
429
            ylim = c(0, max(colSums(counts)) * 1.2), density = if (ncol1) {NULL}
            else {15}, 
            angle = if (ncol1) {NULL}
Stevenn Volant's avatar
Stevenn Volant committed
430
            else {seq(0,160,length.out =nlevels(group[, 2]))[as.integer(group[, 2])]}, col = col[as.integer(group[, 1])], las = 2)
stevenn's avatar
stevenn committed
431
    legend("topright", levels(group[, 1]), fill = col[1:nlevels(group[,1])], bty = "n")
Stevenn Volant's avatar
Stevenn Volant committed
432
    if (!ncol1)  legend("topleft", levels(group[, 2]), density = 15,col = 1, angle = seq(0,160,length.out =nlevels(group[, 2]))[1:nlevels(group[, 2])], bty = "n")
stevenn's avatar
stevenn committed
433
434
435
436
437
  
  }


  ## barplot Nul 
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
438
  barPlotNul <-function (input,counts, group, cex.names = 1, col = c("lightblue","orange", "MediumVioletRed", "SpringGreen")) 
stevenn's avatar
stevenn committed
439
440
441
  {
    
    percentage <- apply(counts, 2, function(x) {sum(x == 0)}) * 100/nrow(counts)
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
442
    percentage.allNull <- (nrow(counts) - nrow(removeNulCounts(counts))) * 100/nrow(counts)
stevenn's avatar
stevenn committed
443
444
    ncol1 <- ncol(group) == 1
    
Stevenn Volant's avatar
Stevenn Volant committed
445
    par(cex=input$cexTitleDiag,mar=c(6,6,4,5))
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
446

stevenn's avatar
stevenn committed
447
448
449
450
    barplot(percentage, las = 2, col = col[as.integer(group[,1])], 
            density = if (ncol1) {NULL}
            else {15}, 
            angle = if (ncol1) {NULL}
Stevenn Volant's avatar
Stevenn Volant committed
451
            else {seq(0,160,length.out =nlevels(group[, 2]))[as.integer(group[, 2])]},
stevenn's avatar
stevenn committed
452
453
454
455
456
457
            cex.names = cex.names, ylab = "Proportion of null counts", 
            main = "Proportion of null counts per sample", 
            ylim = c(0, 1.2 * ifelse(max(percentage) == 0, 1, max(percentage))))
    
    abline(h = percentage.allNull, lty = 2, lwd = 2)
    legend("topright", levels(group[, 1]), fill = col[1:nlevels(group[,1])], bty = "n")
Stevenn Volant's avatar
Stevenn Volant committed
458
    if (!ncol1) legend("topleft", levels(group[, 2]), density = 15, col = 1, angle = seq(0,160,length.out =nlevels(group[, 2]))[1:nlevels(group[, 2])], bty = "n")
stevenn's avatar
stevenn committed
459
460
461
462
  }


  ## Plot density
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
463
  densityPlotTot <-function (input,counts, group, col = c("lightblue", "orange", "MediumVioletRed", "SpringGreen")) 
stevenn's avatar
stevenn committed
464
465
466
467
  {
    
    counts <- removeNulCounts(counts)
    ncol1 <- ncol(group) == 1
Stevenn Volant's avatar
Stevenn Volant committed
468
    par(cex=input$cexTitleDiag,mar=c(8,5,4,5))
stevenn's avatar
stevenn committed
469
470
    plot(density(log2(counts[, 1] + 1)), las = 1, lwd = 2, main = "Density of counts distribution", 
         xlab = expression(log[2] ~ (raw ~ count + 1)), 
471
         ylim = c(0, max(apply(counts, 2, function(x) {max(density(log2(x + 1))$y)})) * 1.05), 
stevenn's avatar
stevenn committed
472
         lty = if (ncol1) {1}
473
         else{rep(seq(1:6),ceiling(nlevels(group[, 2])/6))[as.integer(group[, 2])[1]]}, 
stevenn's avatar
stevenn committed
474
475
476
477
478
479
         col = col[as.integer(group[, 1])[1]])
    
    for (i in 2:ncol(counts)) 
    {
      lines(density(log2(counts[, i] + 1)), col = col[as.integer(group[,1])[i]], lwd = 2, 
            lty = if (ncol1) {1}
480
            else{rep(seq(1:6),ceiling(nlevels(group[, 2])/6))[as.integer(group[, 2])[i]]})
stevenn's avatar
stevenn committed
481
482
    }
    legend("topright", levels(group[, 1]), lty = 1, col = col[1:nlevels(group[,1])], lwd = 2, bty = "n")
483
    if (!ncol1) legend("topleft", levels(group[, 2]), lty = rep(seq(1:6),ceiling(nlevels(group[, 2])/6))[1:nlevels(group[, 2])], col = 1, lwd = 2, bty = "n")
stevenn's avatar
stevenn committed
484
485
486
487
488
    
  }


  ## Table of maj. taxo
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
489
  majTab <- function(input,counts,n)
stevenn's avatar
stevenn committed
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
  {
    seqnames <- apply(counts, 2, function(x) {
      x <- sort(x, decreasing = TRUE)
      names(x)[1:n]
    })
    seqnames <- unique(unlist(as.character(seqnames)))
    sum <- apply(counts, 2, sum)
    counts <- counts[seqnames, ]
    sum <- matrix(sum, nrow(counts), ncol(counts), byrow = TRUE)
    p <- round(100 * counts/sum, digits = 3)
    return(p)
  }


  ## Plot maj. taxo
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
505
  majTaxPlot <-function (input,counts, n = 3, group, cex.names = 1, col = c("lightblue",  "orange", "MediumVioletRed", "SpringGreen")) 
stevenn's avatar
stevenn committed
506
  {
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
507
    p = majTab(input,counts,n)
stevenn's avatar
stevenn committed
508
509
510
    maj <- apply(p, 2, max)
    seqname <- rownames(p)[apply(p, 2, which.max)]
    ncol1 <- ncol(group) == 1
Stevenn Volant's avatar
Stevenn Volant committed
511
    par(cex=input$cexTitleDiag,mar=c(6,6,4,5))
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
512
    x <- barplot(maj, col = col[as.integer(group[, 1])], main = "Proportion of mapped reads from\nmost expressed sequence",
stevenn's avatar
stevenn committed
513
                 ylim = c(0, max(maj) * 1.2), cex.main = 1, 
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
514
                 cex.names = cex.names, las = 2, ylab = "Proportion of mapped reads", 
stevenn's avatar
stevenn committed
515
516
517
                 density = if (ncol1) {NULL}
                 else {15}, 
                 angle = if (ncol1) {NULL}
Stevenn Volant's avatar
Stevenn Volant committed
518
                 else {seq(0,160,length.out =nlevels(group[, 2]))[as.integer(group[, 2])]})
stevenn's avatar
stevenn committed
519
520
521
    
    legend("topright", levels(group[, 1]), fill = col[1:nlevels(group[,1])], bty = "n")
    if (!ncol1) legend("topleft", levels(group[, 2]), density = 15, col = 1, 
Stevenn Volant's avatar
Stevenn Volant committed
522
                       angle = seq(0,160,length.out =nlevels(group[, 2]))[1:nlevels(group[, 2])], bty = "n")
stevenn's avatar
stevenn committed
523
    
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
524
    for (i in 1:length(seqname)) text(x[i], maj[i]/2, seqname[i], cex=input$cexLabelDiag, srt = 90, adj = 0)
stevenn's avatar
stevenn committed
525
526
527
528
  }
  

  ## plot SERE Coefs
Stevenn Volant's avatar
Stevenn Volant committed
529
530
531
532
533
534
535
#   SEREplot<-function(input,counts) 
#   {
#     sere = SEREcoef(counts)
#     hc <- hclust(as.dist(sere), method = "ward.D")
#     plot(hc, las = 2, hang = -1, xlab = "SERE distance, Ward criterion",main = "Cluster dendrogram\non SERE values")
#     
#   }
stevenn's avatar
stevenn committed
536
537
538
539
540
  
  
  ## Get the SERE COEF
  SEREcoef<-function(counts)
  {
Stevenn Volant's avatar
Stevenn Volant committed
541
542
543
544
    counts = as.matrix(counts)
    sere <- matrix(0, ncol = ncol(counts), nrow = ncol(counts))
    for (i in 1:(ncol(counts)-1)) {
      for (j in (i+1):ncol(counts)) {
stevenn's avatar
stevenn committed
545
546
547
        sere[i, j] <- sigfun_Pearson_meta(counts[, c(i, j)])
      }
    }
Stevenn Volant's avatar
Stevenn Volant committed
548
    sere=sere+t(sere)
stevenn's avatar
stevenn committed
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
    colnames(sere) <- rownames(sere) <- colnames(counts)
    sere <- round(sere, digits = 3)
    
    return(sere) 
  }
  

  ## function for the SERE coef
  sigfun_Pearson_meta <- function(observed) {
    laneTotals <- colSums(observed)
    total <- sum(laneTotals)
    fullObserved <- observed[rowSums(observed) > 0, ]
    fullLambda <- rowSums(fullObserved)/total
    fullLhat <- fullLambda > 0
    fullExpected <- outer(fullLambda, laneTotals)
    fullKeep <- which(fullExpected > 0)
    oeFull <- (fullObserved[fullKeep] - fullExpected[fullKeep])^2/fullExpected[fullKeep]
    dfFull <- length(fullKeep) - sum(fullLhat != 0)
    sqrt(sum(oeFull)/dfFull)
  }
stevenn's avatar
stevenn committed
569
570


stevenn's avatar
stevenn committed
571
  ## Plots of size factors
572
  diagSFactors<-function (input,normFactors,counts) 
stevenn's avatar
stevenn committed
573
  {
Stevenn Volant's avatar
Stevenn Volant committed
574
575
    geomeans <- exp(rowMeans(log(counts)))
    samples <- colnames(counts)
Stevenn Volant's avatar
Stevenn Volant committed
576
577
      par(cex=input$cexTitleDiag,mar=c(6,6,4,5))
      plot(normFactors, colSums(counts), pch = 19, las = 1,cex = ifelse(input$addLabelSFact,0,input$cexLabelDiag),
stevenn's avatar
stevenn committed
578
579
           ylab = "Total number of reads", xlab = "Size factors", 
           main = "Diagnostic: size factors vs total number of reads")
Stevenn Volant's avatar
Stevenn Volant committed
580
      if(input$addLabelSFact) text(normFactors,colSums(counts),labels = samples,cex=input$cexLabelDiag)
Stevenn Volant's avatar
Stevenn Volant committed
581
      abline(lm(colSums(counts) ~ normFactors + 0), lty = 2, col = "grey")
stevenn's avatar
stevenn committed
582
  }
stevenn's avatar
stevenn committed
583

Amine  GHOZLANE's avatar
Amine GHOZLANE committed
584
585
586
587
588
  
  ### PCoA
  PCoAPlot_meta <-function (input,dds, group_init,col = c("SpringGreen","dodgerblue","black","firebrick1"), plot = "pcoa") 
  {
    cval=c()
Stevenn Volant's avatar
Stevenn Volant committed
589
    time_set = 0
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
590
591
    # Set of shape
    shape=c(19,17,15,18)
Stevenn Volant's avatar
Stevenn Volant committed
592
    
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
593
594
    ## Var of interest
    VarInt  = input$VarInt
Stevenn Volant's avatar
Stevenn Volant committed
595
    
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
    ## Group
    group = as.character(apply(group_init,1,paste, collapse = "-"))
    
    ## Keep only some sample 
    val = c()
    for(i in 1:length(VarInt))
    { 
      Tinput = paste("input$","Mod",VarInt[i],sep="")
      expr=parse(text=Tinput)
      ## All the modalities for all the var of interest
      val = c(val,eval(expr))
    }
    if(length(VarInt)>1) Kval = apply(expand.grid(val,val),1,paste, collapse = "-")
    else Kval = val
    ind_kept = which(as.character(group)%in%Kval)
Stevenn Volant's avatar
Stevenn Volant committed
611
    
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
612
613
614
615
616
    ## Get the group corresponding to the modalities
    group = group[ind_kept]
    nb = length(unique((group)))
    group = as.factor(group)
    
Stevenn Volant's avatar
Stevenn Volant committed
617
618
619
    if(nlevels(group)!=0)
    { 
      ## Get the norm data
620
621
      counts.norm = as.data.frame(round(counts(dds)))
      if(input$CountsType=="norm") counts.norm = as.data.frame(round(counts(dds, normalized = TRUE)))
Stevenn Volant's avatar
Stevenn Volant committed
622
623
624
625
      # was removed
      counts.norm = counts.norm[,ind_kept]
  
      ## Get the distance
Stevenn Volant's avatar
Stevenn Volant committed
626
627
628
      if(input$DistClust!="sere") dist.counts.norm = vegdist(t(counts.norm), method = input$DistClust)
      if(input$DistClust=="sere") dist.counts.norm = as.dist(SEREcoef(counts.norm))
      
Stevenn Volant's avatar
Stevenn Volant committed
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
      ## Do PCoA
      pco.counts.norm = dudi.pco(d = dist.counts.norm, scannf = FALSE,nf=3)
      
      ## Get eigen values
      eigen=(pco.counts.norm$eig/sum(pco.counts.norm$eig))*100
      
      ## xlim and ylim of the plot
      min = min(pco.counts.norm$li)
      max = max(pco.counts.norm$li)
      
      ## get condition set
      condition_set=val[which(val %in% unique(group_init$condition))]
      time_set=val[which(val %in% unique(group_init$time))]
      
      ## Colors
      if(length(col)<length(condition_set) * length(time_set))# && !input$colorgroup)
      {
        col = rainbow(length(condition_set) * length(time_set))
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
647
      }
Stevenn Volant's avatar
Stevenn Volant committed
648
649
650
651
652
653
654
655
656
657
658
      #else if(length(col)<length(condition_set) * length(time_set) && input$colorgroup){
      #  col = rep(col[1:length(condition_set)], length(time_set))
      #}
      if (length(time_set) == 1 && length(condition_set) <= 4){
        cval = apply(expand.grid(condition_set,time_set),1,paste, collapse = "-")
        cval = sort(cval)
      }
      
      # to reactivate
      #pco.counts.norm$li = pco.counts.norm$li[ind_kept,]
      if (plot == "pcoa"){
Stevenn Volant's avatar
Stevenn Volant committed
659
        par(cex=input$cexTitleDiag,mar=c(6,6,4,5))
Stevenn Volant's avatar
Stevenn Volant committed
660
661
        ## Plot axis, label and circles
        plot(pco.counts.norm$li[1:2], xlab=paste("PC1 : ",round(eigen[1],1),"%") , ylab=paste("PC2 : ",round(eigen[2],1),"%"),
Stevenn Volant's avatar
Stevenn Volant committed
662
             xlim=c(min+0.25*min,max+0.25*max), ylim=c(min-0.1,max+0.1), cex.axis=1, cex.lab=1,lwd=2, type="n",main='Principal Coordinates Analysis ')
Stevenn Volant's avatar
Stevenn Volant committed
663
664
665
666
667
668
669
670
671
        # Set different shapes
        if(input$labelPCOA == "Group"){
          if(!is.null(cval)){
            for (i in 1:length(cval)){
              points(pco.counts.norm$li[which(group==cval[i]),1:2],pch=shape[i],col=col[i], cex=input$cexpoint)
            }
            s.class(dfxy = pco.counts.norm$li, fac = group, col = col, label = levels(group),
                    add.plot = TRUE, cpoint = 0, cell=input$cexcircle, clabel=input$cexLabelDiag,  cstar = input$cexstar)
          }else s.class(dfxy = pco.counts.norm$li, fac = group, col = col, label = levels(group),
672
                        add.plot = TRUE, cpoint = input$cexpoint, cell=input$cexcircle, clabel=input$cexLabelDiag,  cstar = input$cexstar)
Stevenn Volant's avatar
Stevenn Volant committed
673
674
675
676
677
678
679
680
681
682
        }  
        else{
          s.label(pco.counts.norm$li, clabel = input$cexLabelDiag,boxes=FALSE, add.plot = TRUE)
          s.class(dfxy = pco.counts.norm$li, fac = group, col = col, label = levels(group), add.plot = TRUE, cpoint = 0, clabel = 0, cstar = input$cexstar, cell=input$cexcircle)
        }
      }else{
        barplot(eigen[1:7], xlab="Dimensions", ylab="Eigenvalues (%)", names.arg = 1:7, col = c(rep("black", 2), rep("grey", 5)), ylim=c(0,max(eigen)+5), cex.axis=1.2, cex.lab=1.4,cex.names=1.2)
      }
  }
  
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
683
684
685
  }
  
  ### PCA
686
  PCAPlot_meta <-function(input,dds, group_init, n = min(500, nrow(counts(dds))), type.trans = c("VST", "rlog"), 
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
687
688
                           col = c("lightblue", "orange", "MediumVioletRed", "SpringGreen"),plot="pca") 
  {
689
690
    ## Var of interest
    VarInt  = input$VarInt
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
691
    
692
693
    group = as.character(apply(group_init,1,paste, collapse = "-"))
    group_init = group
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
694
    
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
    ## Keep only some sample 
    val = c()
    for(i in 1:length(VarInt))
    { 
      Tinput = paste("input$","Mod",VarInt[i],sep="")
      expr=parse(text=Tinput)
      ## All the modalities for all the var of interest
      val = c(val,eval(expr))
    }
    if(length(VarInt)>1) Kval = apply(expand.grid(val,val),1,paste, collapse = "-")
    else Kval = val
    ind_kept = which(as.character(group)%in%Kval)
    
    ## Get the group corresponding to the modalities
    group = group[ind_kept]
    nb = length(unique((group)))
    group = as.factor(group)
    
    ## To select the colors
    indgrp =as.integer(as.factor(group_init))[ind_kept]
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
715
    
716
717
    
    if(nlevels(group)!=0)
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
718
    { 
719
      type.trans <- type.trans[1]
Stevenn Volant's avatar
Stevenn Volant committed
720
      
721
722
723
      if (type.trans == "VST") counts.trans <- assay(varianceStabilizingTransformation(dds))
      else counts.trans <- assay(rlogTransformation(dds))
      counts.trans = counts.trans[,ind_kept]
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
724
      
725
726
      rv = apply(counts.trans, 1, var, na.rm = TRUE)
      pca = prcomp(t(counts.trans[order(rv, decreasing = TRUE),][1:n, ]))
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
727
      
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
      if(plot=="pca")
      { 
        prp <- pca$sdev^2 * 100/sum(pca$sdev^2)
        prp <- round(prp, 2)
        ncol1 <- ncol(group) == 1
        
        abs = range(pca$x[, 1])
        abs = abs(abs[2] - abs[1])/25
        ord = range(pca$x[, 2])
        ord = abs(ord[2] - ord[1])/25
        
        par(mfrow = c(1, 2),cex=input$cexTitleDiag,mar=c(6,6,4,5))
        plot(pca$x[, 1], pca$x[, 2], las = 1, cex = input$cexTitleDiag, col = col[indgrp], 
             pch = 16,
             xlab = paste0("PC1 (", prp[1], "%)"),
             ylab = paste0("PC2 (", prp[2], "%)"), 
             main = "Principal Component Analysis"
              )
        abline(h = 0, v = 0, lty = 2, col = "lightgray")
        text(pca$x[, 1] - ifelse(pca$x[, 1] > 0, abs, -abs), pca$x[,2] - ifelse(pca$x[, 2] > 0, ord, -ord), colnames(counts.trans), col = col[indgrp],cex=input$cexLabelDiag)
        abs = range(pca$x[, 1])
        abs = abs(abs[2] - abs[1])/25
        ord = range(pca$x[, 3])
        ord = abs(ord[2] - ord[1])/25
        plot(pca$x[, 1], pca$x[, 3], las = 1, cex = input$cexTitleDiag, col = col[indgrp], 
             pch = 16,
             xlab = paste0("PC1 (", prp[1], "%)"), 
             ylab = paste0("PC3 (", prp[3], "%)"), 
             main = "Principal Component Analysis")
        abline(h = 0, v = 0, lty = 2, col = "lightgray")
        text(pca$x[, 1] - ifelse(pca$x[, 1] > 0, abs, -abs), pca$x[,3] - ifelse(pca$x[, 3] > 0, ord, -ord), colnames(counts.trans), col = col[indgrp],cex=input$cexLabelDiag)
      }
      if(plot=="eigen"){eigen = pca$sdev[1:min(7,ncol(counts.trans))]^2; barplot(eigen, xlab="Dimensions", ylab="Eigenvalues (%)", names.arg = 1:min(7,ncol(counts.trans)), col = c(rep("black", 3), rep("grey", 4)), ylim=c(0,max(eigen)+5), cex.axis=1.2, cex.lab=1.4,cex.names=1.2)}
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
761
762
763
764
765
    }
  }
  
  
  
stevenn's avatar
stevenn committed
766

stevenn's avatar
stevenn committed
767
768
769
770
771
  ############################################################
  ##
  ##              CREATE THE CONTRAST DATABASE
  ##
  ############################################################
stevenn's avatar
stevenn committed
772

stevenn's avatar
stevenn committed
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
  
  BaseContrast <- function(input,names,namesfile)
  {  

    v_tmp = c()
    filesize = file.info(namesfile)[,"size"]
    
    for(i in 1:length(names))
    {  
      Tinput = paste("input$",names[i],sep="")
      expr=parse(text=Tinput)
      val = eval(expr) 
      v_tmp[i] = as.numeric(val)
    }
    
    if(filesize!=0)
    { 
      oldContrast = read.table(namesfile,header=TRUE)
      colnamesTmp = c(colnames(oldContrast),input$ContrastName)
      mat = cbind(oldContrast,v_tmp)
    }
    else{ colnamesTmp = input$ContrastName; mat = v_tmp}
    
    write.table(mat,namesfile,row.names=FALSE,col.names = colnamesTmp)
  }
  
  
  ## Remove nul counts
  removeNulCounts <-function (counts) 
  {
    return(counts[rowSums(counts) > 0, ])
  }
stevenn's avatar
stevenn committed
805
806

  
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
807
808
809
810
811
  ############################################################
  ##
  ##              VISUALISATION PLOTS
  ##
  ############################################################
stevenn's avatar
stevenn committed
812
  
813
  GetDataToPlot <- function(input,resDiff,VarInt,ind_taxo,aggregate=TRUE)
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
814
815
  {
    dds = resDiff$dds
816
817
    val = c()
    list.val = list()
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
818
819
820
821
822
    counts = as.data.frame(round(counts(dds, normalized = TRUE)))
    target = resDiff$target
    counts_tmp_combined = NULL
    prop_tmp_combined = NULL
    targetInt = NULL
Stevenn Volant's avatar
Stevenn Volant committed
823
    namesCounts = NULL
824
825
    levelsMod = NULL
    prop_all=NULL
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
826
827
828
    ## Select a subset within the taxonomy level (default is the 12 most abundant)
    nbKept = length(ind_taxo)
    Taxonomy = rownames(counts)
stevenn's avatar
stevenn committed
829
    
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
830
831
    if (length(VarInt)>0 && nbKept>0)
    { 
832
833
      ## Get the modalities to keep
      for(i in 1:length(VarInt))
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
834
      { 
835
836
837
838
839
        Tinput = paste("input$","ModVisu",VarInt[i],sep="")
        expr=parse(text=Tinput)
        ## All the modalities for all the var of interest
        val = c(val,eval(expr))
        list.val[[i]] = eval(expr)
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
840
      }
841
842
      if (!is.null(val) && !is.null(list.val))
      {
svolant's avatar
svolant committed
843

844
845
846
847
848
849
850
851
852
        ## Create the variable to plot
        targetInt = as.data.frame(target[,VarInt])
        rownames(targetInt)=target[,1]  
        ## Combining the Varint
        if(length(VarInt)>1){targetInt$AllVar = apply(targetInt,1,paste, collapse = "-"); targetInt$AllVar = factor(targetInt$AllVar,levels =  expand.grid2.list(list.val))}
        if(length(VarInt)<=1){targetInt$AllVar = target[,VarInt]; targetInt$AllVar = factor(targetInt$AllVar,levels = val)}
        colnames(targetInt) = c(VarInt,"AllVar")
        
        ## Keep only the selected modalities
svolant's avatar
svolant committed
853
854
855
        ind_kept = which(!is.na(targetInt$AllVar))
        targetInt = targetInt[ind_kept,]
          
856
857
858
859
        levelsMod = levels(targetInt$AllVar)
        
        ## Create the counts matrix only for the selected subset
        counts_tmp = counts[Taxonomy%in%ind_taxo,]
svolant's avatar
svolant committed
860
        counts_tmp = counts_tmp[,colnames(counts_tmp)%in%rownames(targetInt)]
861
862
863
864
865
866
867
868
        
        ## Proportions over all the taxonomies
        prop_all = t(counts)/rowSums(t(counts))
        prop_all = as.data.frame(prop_all[,Taxonomy%in%ind_taxo])
        prop_all = as.matrix(prop_all[rownames(prop_all)%in%rownames(targetInt),])
        rownames(prop_all) = targetInt$AllVar
        
        ## Be careful transposition !
svolant's avatar
svolant committed
869
        if(aggregate && nrow(counts_tmp)>0 && nrow(targetInt)>0)
870
871
872
873
874
875
        { 
          counts_tmp_combined = aggregate(t(counts_tmp),by=list(targetInt$AllVar),sum)
          rownames(counts_tmp_combined) = counts_tmp_combined$Group.1
          namesCounts = counts_tmp_combined$Group.1
          counts_tmp_combined = as.matrix(counts_tmp_combined[,-1])
        }
svolant's avatar
svolant committed
876
        if(!aggregate && nrow(counts_tmp)>0 && nrow(targetInt)>0)
877
878
879
880
881
882
883
884
885
        {  
          counts_tmp_combined = t(counts_tmp)
          prop_tmp_combined = counts_tmp_combined/colSums(counts_tmp)
          rownames(counts_tmp_combined) = targetInt$AllVar
          namesCounts = targetInt$AllVar
          rownames(prop_tmp_combined) = targetInt$AllVar
        }
        
        ## Ordering the counts
svolant's avatar
svolant committed
886
887
888
889
890
891
892
893
        if(!is.null(counts_tmp_combined))
        {
          MeanCounts = apply(counts_tmp_combined,2,mean)
          ord = order(MeanCounts,decreasing=TRUE)
          counts_tmp_combined = as.matrix(counts_tmp_combined[,ord])
          if(!aggregate) prop_tmp_combined = as.matrix(prop_tmp_combined[,ord])
          prop_all = as.matrix(prop_all[,ord])
        }
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
894
895
896
      }
    }
    
897
      return(list(counts = counts_tmp_combined,targetInt=targetInt,prop=prop_tmp_combined,namesCounts=namesCounts,levelsMod=levelsMod,prop_all=prop_all))
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
898
899
900
    
    
  }
stevenn's avatar
stevenn committed
901
902
  
  
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
903
904
905
906
907
908
  
  ###########################
  ## Plots for visualisation
  ###########################
  
  Plot_Visu_Barplot <- function(input,resDiff)
stevenn's avatar
stevenn committed
909
  {
910
    
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
911
    ## Get Input for BarPlot
912
913
    VarInt = input$VisuVarInt
    ind_taxo = input$selectTaxoPlot
stevenn's avatar
stevenn committed
914
    
915
    tmp_combined = GetDataToPlot(input,resDiff,VarInt,ind_taxo)
Stevenn Volant's avatar
Stevenn Volant committed
916
    counts_tmp_combined = tmp_combined$counts
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
917
    nbKept = length(ind_taxo)
Stevenn Volant's avatar
Stevenn Volant committed
918
919
920
921
    SamplesNames = tmp_combined$namesCounts
    
    if(nbKept>1) namesTax = colnames(counts_tmp_combined)
    if(nbKept==1) namesTax = ind_taxo
stevenn's avatar
stevenn committed
922
    
923
924
925
926
927
    dataNull = data.frame(x=c(1,2),y=c(1,2))
    plotd3 = nvd3Plot(x ~ y , data = dataNull, type = "multiBarChart", id = 'barplotTaxoNyll',height = input$heightVisu,width=input$widthVisu)
    gg = NULL
    
    if(!is.null(counts_tmp_combined) && nrow(counts_tmp_combined)>0 && length(VarInt)>0)
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
928
    { 
Stevenn Volant's avatar
Stevenn Volant committed
929
      
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
930
931
932
933
      ## Create the data frame for the plot function
      dataBarPlot_mat = c()
      tmp_mat = matrix(0,ncol=3,nrow=nbKept)
      tmp_counts = c()
stevenn's avatar
stevenn committed
934
      
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
935
936
937
        for(i in 1:(nrow(counts_tmp_combined)))
        {
          ## Taxo
Stevenn Volant's avatar
Stevenn Volant committed
938
          tmp_mat[1:nbKept,1] = namesTax
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
939
940
941
942
943
944
945
946
947
948
949
950
          
          ## Counts
  
          tmpProp = counts_tmp_combined[i,]
          if(input$CountsOrProp=="prop")
          { 
            tmpProp = round(tmpProp/sum(tmpProp),3)
            tmpProp = as.numeric(tmpProp/sum(tmpProp) * 100)
          }
          tmp_counts = c(tmp_counts,tmpProp)      
          
          ## Meta data
Stevenn Volant's avatar
Stevenn Volant committed
951
          tmp_mat[1:nbKept,3] = as.character(rep(SamplesNames[i],nbKept))
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
952
953
954
955
956
957
958
959
960
961
962
          
          ## Conbined the sample
          dataBarPlot_mat = rbind(dataBarPlot_mat,tmp_mat)
        }
        
        
        ## Add numeric vector to the dataframe
        dataBarPlot_mat = as.data.frame(dataBarPlot_mat)
        
        colnames(dataBarPlot_mat) = c("Taxonomy","Proportions","AllVar")
        dataBarPlot_mat[,2] = tmp_counts
963
964
965
966
        if(input$SensPlotVisu == "Vertical") Sens = "multiBarChart"
        if(input$SensPlotVisu == "Horizontal") Sens = "multiBarHorizontalChart"
      
        plotd3 <- nvd3Plot(Proportions ~ AllVar | Taxonomy, data = dataBarPlot_mat, type = Sens, id = 'barplotTaxo',height = input$heightVisu,width=input$widthVisu)
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
967
        plotd3$chart(stacked = TRUE)
Stevenn Volant's avatar
Stevenn Volant committed
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
      
        ##################################
        ## Same plot in ggplot2 for export
        ##################################
      
        tax.colors=rep(c("#1f77b4","#aec7e8","#ff7f0e","#ffbb78", "#2ca02c","#98df8a","#d62728","#ff9896","#9467bd","#c5b0d5","#8c564b",
                         "#c49c94","#e377c2","#f7b6d2","#7f7f7f", "#c7c7c7","#bcbd22","#dbdb8d","#17becf","#9edae5"),ceiling(nbKept/20))
        
        dataBarPlot_mat$Taxonomy = factor(dataBarPlot_mat$Taxonomy,levels = namesTax)
      
        gg= ggplot(dataBarPlot_mat, aes(x=AllVar, y=Proportions, fill=Taxonomy)) 
        gg= gg + geom_bar(stat="identity")
        gg= gg + theme_bw()+ scale_fill_manual(values=tax.colors)
        gg = gg +theme(panel.grid.minor.x=element_blank(),panel.grid.major.x=element_blank()) 
        if(input$CountsOrProp=="prop") gg = gg+labs(y="Relative abundance (%)",x="")
        if(input$CountsOrProp=="counts") gg = gg+labs(y="Abundance",x="")
        if(input$SensPlotVisu == "Horizontal") gg = gg + coord_flip()
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
985
    } 
Stevenn Volant's avatar
Stevenn Volant committed
986
    return(list(plotd3=plotd3,gg=gg))
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
987
988
989
990
991
992
993
994
995
996
997
  }
  
  
  
######################################################
##
##            HEATMAP
##
######################################################
  
  
Stevenn Volant's avatar
Stevenn Volant committed
998
  Plot_Visu_Heatmap <- function(input,resDiff,export=FALSE){
Amine  GHOZLANE's avatar
Amine GHOZLANE committed
999
  
1000
  VarInt = input$VisuVarInt