Add plots to analyse the differential abundances

LoadPackages.R

@@ -65,11 +65,11 @@ if(!require(treeWeightD3)){
    devtools::install_github('pierreLec/treeWeightD3')
    library(treeWeightD3)
 }
-if (!require(BiocInstaller)){
-  source("https://bioconductor.org/biocLite.R")
-  biocLite("BiocInstaller")
-  library(BiocInstaller)
-}
+# if (!require(BiocInstaller)){
+#   source("https://bioconductor.org/biocLite.R")
+#   biocLite("BiocInstaller")
+#   library(BiocInstaller)
+# }
 
 if (!require(d3heatmap)) {
   devtools::install_github('aghozlane/d3heatmap')
@@ -223,4 +223,27 @@ if (!require("htmltools")){
 #   library(plotly)  
 # }
 
+if (!require("rAmCharts")){
+  install.packages("rAmCharts")
+  library(rAmCharts)
+}
+
+if(!require("colourpicker")){
+  install.packages("colourpicker")
+  library(coulourpicker)
+}  
+
+if(!require("data.table")){
+  install.packages("data.table")
+  library(data.table)
+}
+
+if(!require("UpSetR")){
+  install.packages("UpSetR")
+  library(UpSetR)
+}
 
+if(!require("ggrepel")){
+  install.packages("ggrepel")
+  library(ggrepel)
+}
\ No newline at end of file

Rfunctions/CompPlot.R

@@ -5,20 +5,22 @@
 ##############################
 ##       HEATMAP
 ##############################
-Plot_Visu_Heatmap_FC <- function(input,BaseContrast,resDiff,export=FALSE){
+Plot_Visu_Heatmap_FC <- function(input, BaseContrast, resDiff, ContrastListDebounce, SelectTaxoPlotCompDebounce, export=FALSE){
   
   res = NULL
-  SelContrast = input$ContrastList_table_FC
+  #SelContrast = input$ContrastList_table_FC
+  SelContrast = ContrastListDebounce()
   selcontrast_matrix = as.matrix(BaseContrast[,SelContrast])
   colnames(selcontrast_matrix) = SelContrast
   log2FC = Get_log2FC_padj(input,selcontrast_matrix,resDiff, info = NULL)$log2FC
   
-  
   if(!is.null(log2FC) && length(SelContrast)>=2)
   { 
     cont = which(colnames(log2FC)%in%SelContrast)
-    log2FC = log2FC[,SelContrast] 
-    ind_taxo = input$selectTaxoPlotComp
+    log2FC = log2FC[,SelContrast]
+    ###
+    ind_taxo = SelectTaxoPlotCompDebounce()
+    #ind_taxo = input$selectTaxoPlotComp
     ind = rownames(log2FC)%in%ind_taxo
     log2FC = as.matrix(log2FC[ind,])
     
@@ -41,22 +43,235 @@ Plot_Visu_Heatmap_FC <- function(input,BaseContrast,resDiff,export=FALSE){
   return(res)
 }
 
+##############################
+##      P VALUE DENSITY PLOT
+##############################
+Plot_pValue_Density <- function(input, BaseContrast, resDiff, ContrastListDebounce){
+  res = NULL
+  #SelContrast = input$ContrastList_table_FC
+  SelContrast = ContrastListDebounce()
+  if(length(SelContrast)>=1){
+  selcontrast_matrix = as.matrix(BaseContrast[,SelContrast])
+  colnames(selcontrast_matrix) = SelContrast
+  
+  padj = data.frame(na.omit(Get_log2FC_padj(input,selcontrast_matrix,resDiff)$padj))
+  
+  data <- data.frame()
+  for (cont in SelContrast){
+    data_cont <- data.frame(padj[,cont])
+    data_cont$contrast <- cont
+    colnames(data_cont) <- c("padj","contrast")
+    data <- rbind(data, data_cont)
+  }
+  data$contrast <- factor(data$contrast, levels = SelContrast)
+  p <- ggplot(data, aes(x = padj, color = contrast, fill = contrast)) + geom_density(alpha = input$fillOpacity, size = input$lineWidth)
+  p <- p + theme(axis.title = element_text(size = input$FontSizepValueDensity), 
+                 axis.text = element_text(size = input$FontSizepValueDensity),
+                 legend.title = element_text(size = input$FontSizepValueDensity),
+                 legend.text = element_text(size = input$FontSizepValueDensity))
+  res = p
+  }
+  return(res)
+}
+
+##############################
+##      LOGIT PLOT
+##############################
+Plot_Comp_Logit <- function(input, BaseContrast, resDiff, SelectTaxoPlotCompDebounce, export = FALSE)
+  {
+  plot = NULL
+  if(input$Contrast1 != input$Contrast2 && input$Contrast1 != "..." && input$Contrast2 != "..."){
+    SelContrast = c(input$Contrast1, input$Contrast2)
+    selcontrast_matrix = as.matrix(BaseContrast[,SelContrast])
+    colnames(selcontrast_matrix) = SelContrast
+    padj = na.omit(Get_log2FC_padj(input,selcontrast_matrix,resDiff)$padj)
+   
+    #data to plot
+    x <- padj[,input$Contrast1]
+    logit_x <- log(x/(1-x))
+    y <- padj[,input$Contrast2]
+    logit_y <- log(y/(1-y))
+    
+    #classification for colours
+    alpha_logit <- log(as.numeric(input$AlphaVal)/(1-as.numeric(input$AlphaVal)))
+    x_signif <- logit_x < alpha_logit
+    y_signif <- logit_y < alpha_logit
+    matrix <- cbind(x_signif, y_signif)
+    color_var <-
+      apply(matrix, 1, function(point)
+        if (point[[1]]) {
+          if (point[[2]]) {"Significant for both contrasts"}
+          else {"Significant for contrast 1"}
+        } 
+        else{
+          if (point[[2]]) {"Significant for contrast 2"}
+          else{"Not significant"}
+          }
+      )
+    
+    #labels
+    names <- row.names(padj)
+    ###
+    points_to_label <- SelectTaxoPlotCompDebounce()
+    #points_to_label <- input$selectTaxoPlotComp
+    labels <-
+      sapply(names, function(name)
+        if (is.element(name, points_to_label)) {
+          points_to_label[match(name, points_to_label)]
+        } else{
+          ""
+        })
+    
+    if(!export){
+      plot <- scatterD3(x = logit_x, 
+                        y = logit_y,
+                        xlab = paste("logit p value", input$Contrast1),
+                        ylab = paste("logit p value", input$Contrast2),
+                        hover_opacity = 1,
+                        tooltip_text = row.names(padj),
+                        # xlim = xlimits,
+                        # ylim = ylimits,
+                        lines = if (input$showSignifThresholdsLogitPlot) {
+                          if(input$showDiagonal){
+                            data.frame(
+                              slope = c(Inf, 0, 1),
+                              intercept = c(alpha_logit,alpha_logit,0),
+                              stroke = "black",
+                              stroke_width = input$linesWidthLogitPlot,
+                              stroke_dasharray = 5
+                            )
+                          }
+                          else{
+                          data.frame(
+                            slope = c(Inf, 0),
+                            intercept = c(alpha_logit,alpha_logit),
+                            stroke = "black",
+                            stroke_width = input$linesWidthLogitPlot,
+                            stroke_dasharray = 5
+                          )}
+                        }
+                        else{if (input$showDiagonal){
+                          data.frame(
+                            slope = c(1),
+                            intercept = c(0),
+                            stroke = "black",
+                            stroke_width = input$linesWidthLogitPlot,
+                            stroke_dasharray = 5
+                          )}
+                        else{NULL}},
+                        fixed = input$fixed11,
+                        col_var = color_var,
+                        col_lab = if(input$showLegendLogit){"Legend"}else{NA}
+                        ,
+                        colors = c(
+                          "Not significant" = input$colour01,
+                          "Significant for contrast 1" = input$colour02,
+                          "Significant for contrast 2" = input$colour03,
+                          "Significant for both contrasts" = input$colour04
+                        ),
+                        lab = labels,
+                        point_opacity = input$pointOpacityLogit,
+                        point_size = input$pointSizeLogit,
+                        axes_font_size = paste(input$axisFontSizeLogit, "0%", sep = ""),
+                        labels_size = input$labelsSizeLogit,
+                        legend_font_size = paste(input$legendFontSizeLogit, "0%", sep = ""), 
+                        legend_width = input$legendWidthLogit,
+                        # dom_id_reset_zoom = "scatterD3-reset-zoomLogit",
+                        # dom_id_svg_export = "scatterD3-svg-exportLogit",
+                        # menu = FALSE,
+                        # height = input$heightVolcanoPlot,
+                        # width=if(input$modifwidthVolcano){input$widthVolcanoPlot},
+                        # disable_wheel = TRUE
+                        )
+      }
+    
+    if(export){
+      data <- data.frame(logit_x, logit_y, color_var)
+      plot <- ggplot(data, aes(x = logit_x, y = logit_y))
+      plot <- plot + geom_point(aes(color = color_var), size = (input$pointSizeLogit %/% 30), alpha = input$pointOpacityLogit)
+      plot <- plot + scale_x_continuous(paste("logit p value", input$Contrast1))
+      plot <- plot + scale_y_continuous(paste("logit p value", input$Contrast2))
+      plot <- plot + theme(axis.title = element_text(size = input$axisFontSizeLogit), 
+                           axis.text = element_text(size = input$axisFontSizeLogit),
+                           legend.title = element_text(size = input$legendFontSizeLogit),
+                           legend.text = element_text(size = input$legendFontSizeLogit))
+      plot <- plot + geom_text(aes(label = labels, color = color_var), size = input$labelsSizeLogit %/% 3, hjust = 1, vjust = 1, show.legend = FALSE)
+      plot <- plot + scale_color_manual(name = "Legend title", values = c("Not significant" = input$colour01, "Significant for contrast 1" = input$colour02,
+                                                                          "Significant for contrast 2" = input$colour03, "Significant for both contrasts" = input$colour04))
+      if (input$fixed11) {
+        plot <- plot + coord_fixed()
+      }
+      if (input$showSignifThresholdsLogitPlot) {
+        plot <- plot + geom_hline(yintercept = alpha_logit, linetype = "dashed", size = input$linesWidthLogitPlot/2)
+        plot <- plot + geom_vline(xintercept = alpha_logit, linetype = "dashed", size = input$linesWidthLogitPlot/2)
+      }
+      if(input$showDiagonal){
+        plot <- plot + geom_abline(slope = 1, intercept = 0, linetype = "dashed", size = input$linesWidthLogitPlot/2)
+      }
+  }
+  return(plot)
+}}
+
 
 ##############################
 ##      VENN DIAGRAM
 ##############################
-Plot_Visu_Venn <- function(input,BaseContrast,resDiff,export=FALSE){
-  
+Plot_Visu_Venn <- function(input,BaseContrast,resDiff, ContrastListDebounce, export=FALSE){
   res = NULL
-  SelContrast = input$ContrastList_table_FC
+  #SelContrast = input$ContrastList_table_FC
+  SelContrast = ContrastListDebounce()
   data = GetData_venn(input,SelContrast,BaseContrast,resDiff)$res
   res = venn_tooltip(d3vennR(data=data))
-  
   return(res)
 }
 
+##############################
+##      CONTRASTS COMPAIR
+##############################
+Plot_UpSet <- function(input,BaseContrast, resDiff, ContrastListDebounce, export=FALSE){
+  plot = NULL
+  #SelContrast = input$ContrastList_table_FC
+  SelContrast = ContrastListDebounce()
+  if(length(SelContrast)>=2){
+    data = GetData_venn(input,SelContrast,BaseContrast,resDiff)$df.tot
+    listInput <- list()
+    n <- ncol(data)
+    for(i in 1:n){
+      new <- list(na.omit(data[,i]))
+      listInput <- c(listInput, new)
+    }
+    print(input$pointSizeContrastsCompair)
+    print(class(input$pointSizeContrastsCompair))
+    names(listInput) <- colnames(data)
+    plot <- upset(fromList(listInput), nsets = n, nintersects = NA, set_size.show = input$showNumbers, 
+                  show.numbers = if(input$showNumbers){"yes"}else{"no"}, sets = SelContrast, keep.order = !(input$orderBySize),
+                  point.size = input$pointSizeContrastsCompair, line.size = 1.5, order.by = if(input$orderBySize){"freq"}else{"degree"}, 
+                  mainbar.y.label = "Number of differential features in common", sets.x.label = "Number of differential features", 
+                  text.scale = c(2, 1.6, 1.6, 1.6, 2, 1.5))
+  }
+  return(plot)
+}
 
 
+##############################
+##   MULTIPLE VENN COMPAIR
+##############################
+Plot_MultipleVenn <- function(input,BaseContrast, resDiff, ContrastListDebounce){
+  plot = NULL
+  SelContrast = ContrastListDebounce()
+  if(length(SelContrast)>=2){
+  data = GetData_venn(input,SelContrast,BaseContrast,resDiff)$res_multiple_venn
+
+  plot <- ggplot(data,aes(x=x,y=y)) + geom_point() + theme_bw() + 
+          geom_label_repel(aes(label=name), size = input$labelSizemultipleVenn) + xlim(c(0,1)) + ylim(c(0,1)) + 
+          xlab(bquote(Contrast1 *intersect(Contrast2) ~ "/ Contrast1")) + 
+          ylab(bquote(Contrast1 *intersect(Contrast2) ~ "/ Contrast2"))
+  plot <- plot + theme(axis.title = element_text(size = input$FontSizeMultipleVenn), 
+                 axis.text = element_text(size = input$FontSizeMultipleVenn - 2))
+  }
+  return(plot)
+}
+
 
 ##############################################################
 ##
@@ -77,6 +292,7 @@ Go_data_top <- function(vect)
 
 
 ## Get the data to plot the heatmap with the FC
+## and for the logit plot
 Get_log2FC_padj <-function(input,BaseContrast,resDiff, info = NULL)
 {
   log2FC = NULL
@@ -100,17 +316,20 @@ Get_log2FC_padj <-function(input,BaseContrast,resDiff, info = NULL)
                                 independentFiltering=input$IndFiltering,alpha=alpha)
     }
     log2FC = as.matrix(round(result[[SelContrast[1]]][, "log2FoldChange"], 3))
-    padj = as.matrix(round(result[[SelContrast[1]]][, "padj"], 3))
+    #padj = as.matrix(round(result[[SelContrast[1]]][, "padj"], 3))
+    padj = as.matrix((result[[SelContrast[1]]][, "padj"]))
     if(nbCont>1)
     {
       for(i in 2:nbCont)
       {
         log2FC = cbind(log2FC,round(result[[SelContrast[i]]][, "log2FoldChange"], 3))
-        padj = cbind(padj,round(result[[SelContrast[i]]][, "padj"], 7))
+        #padj = cbind(padj,round(result[[SelContrast[i]]][, "padj"], 7))
+        padj = cbind(padj,(result[[SelContrast[i]]][, "padj"]))
       }
+    }#  
       colnames(log2FC) = names(result)
       colnames(padj) = names(result)
-    }
+    #}
     
     rownames(log2FC) = rownames(result[[SelContrast[1]]])
     rownames(padj) = rownames(result[[SelContrast[1]]])
@@ -188,6 +407,7 @@ venn_tooltip <- function(venn){
 GetData_venn <-function(input,SelContrast,BaseContrast,resDiff)
 {
   res = list()
+  res_multiple_venn = NULL
   df.tot = NULL
   VarInt = input$VarInt
   dds = resDiff$dds
@@ -208,11 +428,12 @@ GetData_venn <-function(input,SelContrast,BaseContrast,resDiff)
                                 independentFiltering=input$IndFiltering,alpha=alpha)
     }
     padj = round(result[[SelContrast[1]]][, "padj"], 3)
+    
     # save(result,padj,SelContrast,file = "test1.RData")
     df = as.matrix(rownames(result[[1]]))
-    if(length(which(padj>alpha))>0) df[which(padj>alpha),]=NA 
-    if(any(is.na(padj))) df[which(is.na(padj)),]=NA 
-    
+    if(length(which(padj>alpha))>0) df[which(padj>alpha),]=NA
+    if(any(is.na(padj))) df[which(is.na(padj)),]=NA
+
     if(nbCont>1)
     {
       for(i in 2:nbCont)
@@ -226,17 +447,16 @@ GetData_venn <-function(input,SelContrast,BaseContrast,resDiff)
       colnames(df) = SelContrast
       df = as.data.frame(df)
     }
-    
     ## Keep the entire dataframe
     df.tot = as.data.frame(apply(df,2,Go_data_top))
     maxRow = max(apply(as.data.frame(apply(df,2,Go_data_top)),2,FUN=function(x) length(which(!is.na(x)))))
-    
     df.tot = df.tot[1:max(maxRow,1),]
     ## Remove col with only NA
     df = df[,which(apply(!is.na(df),2,any))]
     
     ncont = ncol(as.data.frame(df))
     names.df = names(df)
+    rownames_multiple_venn = c()
     cmp=1
     if(ncont>1 && !is.null(ncont))
     {
@@ -247,11 +467,21 @@ GetData_venn <-function(input,SelContrast,BaseContrast,resDiff)
           if(i!=j) res[[cmp]] = list(sets=list(names.df[i],names.df[j]),size= length(which(!is.na(intersect(df[,i],df[,j])))))
           if(i==j) res[[cmp]] = list(sets=list(names.df[i]),size= length(which(!is.na(intersect(df[,i],df[,i])))))
           cmp=cmp+1
+          if(i!=j) {
+            if (is.null(res_multiple_venn))
+              {res_multiple_venn <- data.frame(name = c(rownames_multiple_venn, paste(names.df[i],names.df[j], sep = " vs ")), 
+                                                   x = c(length(which(!is.na(intersect(df[,i],df[,j]))))/length(which(!is.na(intersect(df[,i],df[,i]))))), 
+                                                   y = c(length(which(!is.na(intersect(df[,i],df[,j]))))/length(which(!is.na(intersect(df[,j],df[,j]))))),
+                                                   stringsAsFactors=FALSE)}
+            else{new_row <- data.frame(name = c(rownames_multiple_venn, paste(names.df[i],names.df[j], sep = " vs ")),
+                         x = c(length(which(!is.na(intersect(df[,i],df[,j]))))/length(which(!is.na(intersect(df[,i],df[,i]))))),
+                         y = c(length(which(!is.na(intersect(df[,i],df[,j]))))/length(which(!is.na(intersect(df[,j],df[,j]))))))
+                res_multiple_venn <- rbind(res_multiple_venn, new_row)}
+          }
         }
       }
     }
-    
   }
-  return(list(res=res,df.tot=df.tot))
+  return(list(res=res,df.tot=df.tot,res_multiple_venn=res_multiple_venn))
 }
 

Rfunctions/TablesPlot.R

+#@ This file contains all the functions for the
+#@ plots in tab "Tables" of SHAMAN
+
+
+###########################
+##        Bar chart
+###########################
+Bar_Chart_Tables <- function(input, data, export = FALSE) {
+  data_table <- data.table::data.table(data$complete)
+  data_table <-
+    data_table[order(eval(parse(text = input$ColumnOrder)), decreasing = input$Decreasing)]
+  
+  data_table <-
+    data_table[data_table$Id %in% input$selectTaxoPlotBarChart,]
+  
+  #set colors
+  x <- data_table[["log2FoldChange"]]
+  y <- data_table[["pvalue_adjusted"]]
+  y_log <- -1 * log10(y)
+  y_signif <- -log10(as.numeric(input$AlphaVal))
+  x_signif <-
+    if (!input$setXsignifThreshold) {
+      0
+    } else{
+      input$XsignifThreshold
+    }
+  x_up <- x > x_signif
+  x_down <- x < -x_signif
+  y_sig <- y_log > y_signif
+  matrix <- cbind(x_up, x_down, y_sig)
+  Group <- 
+    apply(matrix, 1, function(point)
+    if ((anyNA(point)) |
+        (!point[[3]])) {
+      "Not significant"
+    } else{
+      if (point[[1]]) {
+        "Up"
+      } else{
+        if (point[[2]]) {
+          "Down"
+        } else{
+          "Not significant"
+        }
+      }
+    })
+  colors <-
+    apply(matrix, 1, function(point)
+      if ((anyNA(point)) |
+          (!point[[3]])) {
+        input$colour2
+      } else{
+        if (point[[1]]) {
+          input$colour3
+        } else{
+          if (point[[2]]) {
+            input$colour1
+          } else{
+            input$colour2
+          }
+        }
+      })
+  data_table[["color"]] <- colors
+  data_table[["Group"]] <- Group
+  
+  if(!export){
+    plot <- amBarplot(
+      x = "Id",
+      y = "log2FoldChange",
+      data = data_table,
+      horiz = TRUE,
+      ylab = "log 2 Fold Change",
+      creditsPosition = "bottom-right"
+    )
+    plot <- setProperties(plot, fontSize = input$fontSize)
+    return(plot)
+  }
+  
+  if(export){
+    print(input$fontSize)
+    data <- data_table
+    data$Id <- factor(data$Id, levels = rev(data$Id))
+    plot <- ggplot(data, aes(x = Id, y = log2FoldChange))
+    plot <- plot + geom_bar(stat = "identity", aes(color = Group, fill = Group))
+    plot <- plot + coord_flip()
+    plot <- plot + scale_color_manual(values = c("Down" = input$colour1,"Not significant" = input$colour2,"Up" = input$colour3))
+    plot <- plot + scale_fill_manual(values = c("Down" = input$colour1,"Not significant" = input$colour2,"Up" = input$colour3))
+    plot <- plot + labs(x = "", y = "log 2 Fold Change")
+    plot <- plot + theme(axis.title = element_text(size = input$fontSize),
+                         axis.text = element_text(size = input$fontSize),
+                         legend.title = element_text(size = input$fontSize),
+                         legend.text = element_text(size = input$fontSize))
+    if(input$removeLabelsBarChart){print("ici")
+                                         plot <- plot + theme(axis.ticks.y = element_blank(),
+                                                              axis.text.y = element_blank())}
+    return(plot)
+  }
+}
+
+###########################
+##        Volcano plot
+###########################
+Volcano_Plot <- function(input, data, export = FALSE) {
+
+  x <- data$complete[["log2FoldChange"]]
+  y <- data$complete[["pvalue_adjusted"]]
+  y_log <- -1 * log10(y)
+  
+  #set axes limits, horizontally centering the chart
+  xmin <- -1 * max(-1 * min(x[!is.na(x)]), max(x[!is.na(x)]))
+  xmax <- -xmin
+  deltaX <- xmax - xmin
+  ymin <- min(y_log[!is.na(y_log)])
+  ymax <- max(y_log[!is.na(y_log)])
+  deltaY <- ymax - ymin
+  margin <- 0.1
+  xlimits <- c(xmin - margin * deltaX, xmax + margin * deltaX)
+  ylimits <- c(ymin - margin * deltaY / 10, ymax + margin * deltaY)
+  
+  #significance thresholds
+  y_signif <- -log10(as.numeric(input$AlphaVal))
+  x_signif <-
+    if (!input$setXsignifThreshold) {
+      0
+    } else{
+      input$XsignifThreshold
+    }
+  
+  #comparison to thresholds
+  x_up <- x > x_signif
+  x_down <- x < -x_signif
+  y_sig <- y_log > y_signif
+  matrix <- cbind(x_up, x_down, y_sig)
+  color_var <-
+    apply(matrix, 1, function(point)
+      if ((anyNA(point)) |
+          (!point[[3]])) {
+        "Not significant"
+      } else{
+        if (point[[1]]) {
+          "Up"
+        } else{
+          if (point[[2]]) {
+            "Down"
+          } else{
+            "Not significant"
+          }
+        }
+      })
+  
+  names <- data$complete[["Id"]]
+  points_to_label <- input$selectTaxoLabelVolcano
+  labels <-
+    sapply(names, function(name)
+      if (is.element(name, points_to_label)) {
+        points_to_label[match(name, points_to_label)]
+      } else{
+        ""
+      })
+  # labels <-
+  #   sapply(names, function(name)
+  #     if (is.element(name, points_to_label)) {
+  #       name
+  #     } else{
+  #       ""
+  #     })
+  if(!export){
+  plot <-
+    scatterD3(
+      x = x,
+      y = y_log,
+      xlab = "Log 2 Fold Change",