Ajout css, choix axes pca/pcoa

LoadPackages.R

@@ -131,6 +131,19 @@ if (!require(scales)) {
   library(scales)
 }
 
+
+if(!require(biomformat)){
+  devtools::install_github("biomformat", "aghozlane")
+}
+
+
 # Allow to upload 50M files
 options(shiny.maxRequestSize=50*1024^2) 
-source("internal.R")
+source("Rfunctions/Data_Management.R")
+source("Rfunctions/Stat_Model.R")
+source("Rfunctions/DiagPlot.R")
+source("Rfunctions/VisuPlot.R")
+source("Rfunctions/CompPlot.R")
+source("Rfunctions/DiffTable.R")
+
+

Rfunctions/CompPlot.R

+#@ This file contains all the functions for the 
+#@ comparison plots of SHAMAN
+
+
+##############################
+##       HEATMAP
+##############################
+Plot_Visu_Heatmap_FC <- function(input,BaseContrast,resDiff,export=FALSE){
+  
+  res = NULL
+  SelContrast = input$ContrastList_table_FC
+  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
+    ind = rownames(log2FC)%in%ind_taxo
+    log2FC = as.matrix(log2FC[ind,])
+    
+    if(input$SortHeatComp =="Selection") tmp_ord = match(ind_taxo, rownames(log2FC))
+    if(input$SortHeatComp =="Names") tmp_ord = order(rownames(log2FC))
+    if(input$SortHeatComp =="Values") tmp_ord = order(log2FC[,1])
+    
+    if(input$SortHeatComp !="Auto") log2FC = log2FC[tmp_ord,]
+    
+    col1 <- c(colorRampPalette(c("royalblue4","royalblue3","royalblue2","royalblue1","white"))(n = 100),colorRampPalette(c("white",  "firebrick1", "firebrick2", "firebrick3", "firebrick4"))(n = 100))
+    breaks <- c(seq(min(log2FC,-0.01,na.rm = TRUE), 0,length=100),seq(0.01,max(log2FC,0.02,na.rm = TRUE),length=100))
+    colorFunc <- col_bin(col1, bins = rescale(breaks))
+    ## Transpose matrix if Horizontal
+    if(input$SensPlotVisuComp=="Horizontal") log2FC = t(as.matrix(log2FC))
+    
+    if(!export && nrow(log2FC)>0) res = d3heatmap(log2FC, dendrogram = "none", Rowv = (input$SortHeatComp =="Auto"), Colv = FALSE, na.rm = TRUE, height = input$heightVisuComp, show_grid = FALSE, colors = colorFunc, scale = input$scaleHeatmapComp,cexRow = input$LabelSizeHeatmapComp,cexCol =input$LabelSizeHeatmapComp, offsetCol=input$LabelColOffsetHeatmapComp,offsetRow=input$LabelRowOffsetHeatmapComp)
+    if(export && nrow(log2FC)>0) heatmap.2(log2FC, dendrogram = "none", Rowv = (input$SortHeatComp =="Auto"), Colv = FALSE, na.rm = TRUE, margins=c(input$lowerMarginComp,input$rightMarginComp), density.info="none", trace="none", col = col1, scale = input$scaleHeatmapComp,cexRow = input$LabelSizeHeatmapComp,cexCol =input$LabelSizeHeatmapComp, 
+                                           offsetCol=input$LabelColOffsetHeatmapComp,offsetRow=input$LabelRowOffsetHeatmapComp,symm=FALSE,symkey=TRUE,symbreaks=TRUE)
+  }
+  return(res)
+}
+
+
+##############################
+##      VENN DIAGRAM
+##############################
+Plot_Visu_Venn <- function(input,BaseContrast,resDiff,export=FALSE){
+  
+  res = NULL
+  SelContrast = input$ContrastList_table_FC
+  data = GetData_venn(input,SelContrast,BaseContrast,resDiff)$res
+  res = venn_tooltip(d3vennR(data=data))
+  
+  return(res)
+}
+
+
+
+
+##############################################################
+##
+##          Useful functions
+##
+##############################################################
+
+## Get the non NA data at the top of the dataframe
+Go_data_top <- function(vect)
+{
+  n = length(vect)
+  tmp = rep(NA,n)
+  ind = which(!is.na(vect))
+  n1 = length(ind)
+  if(n1>0)    tmp[1:n1] =  vect[ind]
+  return(tmp)
+}
+
+
+## Get the data to plot the heatmap with the FC
+Get_log2FC_padj <-function(input,BaseContrast,resDiff, info = NULL)
+{
+  log2FC = NULL
+  padj = NULL
+  dds = resDiff$dds
+  counts = resDiff$counts
+  target = resDiff$target
+  SelContrast = colnames(BaseContrast)
+  nbCont = length(SelContrast)
+  result = list()
+  alpha = as.numeric(input$AlphaVal)
+  cooksCutoff = ifelse(input$CooksCutOff!='Auto',ifelse(input$CooksCutOff!=Inf,input$CutOffVal,Inf),TRUE)
+  
+  if(nbCont>=1)
+  {
+    for(i in 1:nbCont)
+    { 
+      cont = as.character(SelContrast[i])
+      result[[cont]] <- results(dds,contrast=BaseContrast[,cont],pAdjustMethod=input$AdjMeth,
+                                cooksCutoff=cooksCutoff,
+                                independentFiltering=input$IndFiltering,alpha=alpha)
+    }
+    log2FC = as.matrix(round(result[[SelContrast[1]]][, "log2FoldChange"], 3))
+    padj = as.matrix(round(result[[SelContrast[1]]][, "padj"], 3))
+    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))
+      }
+      colnames(log2FC) = names(result)
+      colnames(padj) = names(result)
+    }
+    
+    rownames(log2FC) = rownames(result[[SelContrast[1]]])
+    rownames(padj) = rownames(result[[SelContrast[1]]])
+  }
+  return(list(log2FC=as.data.frame(log2FC),padj=padj))
+}
+
+
+## Add tooltips on venn digramm
+venn_tooltip <- function(venn){
+  venn$x$tasks[length(venn$x$tasks)+1] <- list(
+    htmlwidgets::JS('
+                    function(){
+                    var div = d3.select(this);
+                    
+                    // add a tooltip
+                    var tooltip = d3.select("body").append("div")
+                    .attr("class", "venntooltip")
+                    .style("position", "absolute")
+                    .style("text-align", "center")
+                    .style("width", 128)
+                    .style("height", 16)
+                    .style("background", "#333")
+                    .style("color","#ddd")
+                    .style("padding","2px")
+                    .style("border","0px")
+                    .style("border-radius","8px")
+                    .style("opacity",0);
+                    
+                    div.selectAll("path")
+                    .style("stroke-opacity", 0)
+                    .style("stroke", "#fff")
+                    .style("stroke-width", 0)
+                    
+                    // add listeners to all the groups to display tooltip on mousover
+                    div.selectAll("g")
+                    .on("mouseover", function(d, i) {
+                    
+                    // sort all the areas relative to the current item
+                    venn.sortAreas(div, d);
+                    
+                    // Display a tooltip with the current size
+                    tooltip.transition().duration(400).style("opacity", .9);
+                    tooltip.text(d.size);
+                    
+                    // highlight the current path
+                    var selection = d3.select(this).transition("tooltip").duration(400);
+                    selection.select("path")
+                    .style("stroke-width", 3)
+                    .style("fill-opacity", d.sets.length == 1 ? .4 : .1)
+                    .style("stroke-opacity", 1);
+                    })
+                    
+                    .on("mousemove", function() {
+                    tooltip.style("left", (d3.event.pageX) + "px")
+                    .style("top", (d3.event.pageY - 28) + "px");
+                    })
+                    
+                    .on("mouseout", function(d, i) {
+                    tooltip.transition().duration(50).style("opacity", 0);
+                    var selection = d3.select(this).transition("tooltip").duration(400);
+                    selection.select("path")
+                    .style("stroke-width", 0)
+                    .style("fill-opacity", d.sets.length == 1 ? .25 : .0)
+                    .style("stroke-opacity", 0);
+                    });
+                    }
+                    ')
+    )
+  return(venn)
+  }
+
+
+## Transform the data for the venn diagram
+GetData_venn <-function(input,SelContrast,BaseContrast,resDiff)
+{
+  res = list()
+  df.tot = NULL
+  VarInt = input$VarInt
+  dds = resDiff$dds
+  counts = resDiff$counts
+  target = resDiff$target
+  nbCont = length(SelContrast)
+  result = list()
+  alpha = as.numeric(input$AlphaVal)
+  cooksCutoff = ifelse(input$CooksCutOff!='Auto',ifelse(input$CooksCutOff!=Inf,input$CutOffVal,Inf),TRUE)
+  
+  if(nbCont>=2)
+  {
+    for(i in 1:nbCont)
+    { 
+      cont = as.character(SelContrast[i])
+      result[[cont]] <- results(dds,contrast=BaseContrast[,cont],pAdjustMethod=input$AdjMeth,
+                                cooksCutoff=cooksCutoff,
+                                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(nbCont>1)
+    {
+      for(i in 2:nbCont)
+      {
+        padj = round(result[[SelContrast[i]]][, "padj"], 3)
+        df.tmp = as.matrix(rownames(result[[i]]))
+        if(length(which(padj>alpha))>0) df.tmp[which(padj>alpha),]=NA 
+        if(any(is.na(padj))) df.tmp[which(is.na(padj)),]=NA  
+        df = cbind(df,df.tmp)
+      }
+      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)
+    cmp=1
+    if(ncont>1 && !is.null(ncont))
+    {
+      for(i in 1:(ncont))
+      {
+        for(j in i:ncont)
+        {
+          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
+        }
+      }
+    }
+    
+  }
+  return(list(res=res,df.tot=df.tot))
+}
+

Rfunctions/Data_Management.R

+#@ This file contains all the functions needed to
+#@ to load, check and transform the data
+
+## Add news to the home page
+addNews <- function(date ="",title="",text="")
+{
+  res=list()
+  res$r1 = paste("<b><font size='+1'>",date,"</font></b>", " - ", "<b><font size='+1'>",title,"</font></b><br/>")
+  res$r2 = paste("<p><font color='grey'>",text,"</font></p><hr/>")
+  
+  return(HTML(unlist(res)))
+}
+
+
+
+## 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)
+}
+
+
+
+
+
+## Check the format of the counts table
+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))
+}
+
+## Check the format of the taxonomy table
+CheckTaxoTable <- function(taxo,counts)
+{
+  Error = NULL
+  Warning = NULL
+  if(ncol(taxo)<1){Error = "The number of columns of the taxonomy table must be at least 1" }
+  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"}
+  }
+  
+  ## Annotated features without counts
+  if(!any(rownames(taxo)%in%rownames(counts))){ Error = "Some annotated features are not in the count table"}
+  
+  return(list(Error=Error,Warning=Warning))
+}
+
+## Get the percentage of annotated OTU
+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))
+}
+
+
+## Get the counts and the taxo tables from the BIOM format file.
+GetDataFromBIOM <-function(dataBIOM)
+{
+  ## Counts table
+  counts = biom_data(dataBIOM)
+  counts = as.matrix(counts)
+  ## Change of - to . is risky
+  colnames(counts) = gsub("-",".",colnames(counts))
+  counts = as.data.frame(counts)
+  CheckCounts = CheckCountsTable(counts)
+  counts = CheckCounts$counts
+  
+  ## Taxonomy table
+  taxo = as.data.frame(observation_metadata(dataBIOM))
+  OTUnames = rownames(taxo)
+  ## Modif taxo table (remove p__,... and change the colnames)
+  taxo = as.data.frame(sapply(taxo,gsub,pattern="^.*__",replacement=""))
+  colnames(taxo) = c("Kingdom", "Phylum","Class","Order","Family","Genus","Species")
+  rownames(taxo) = OTUnames
+  ## Remove empty row
+  taxo[taxo==""] = NA
+  taxo[taxo=="Unassigned"] = NA
+  taxo=taxo[rowSums(is.na(taxo))!=dim(taxo)[2], ]
+  
+  CheckTaxo = CheckTaxoTable(taxo,counts)
+  
+  ## Pourcentage of annotation
+  tmp = PercentAnnot(counts,taxo)
+  
+  return(list(counts=counts,taxo=taxo,CheckCounts=CheckCounts,CheckTaxo=CheckTaxo,Percent=tmp$Percent,CheckPercent=tmp$Error))
+}
+
+## Check the data
+GetDataFromCT <-function(dataC,dataT)
+{
+  
+  ## Counts table
+  counts = dataC
+  CheckCounts = CheckCountsTable(counts)
+  counts = CheckCounts$counts
+  ## Taxonomy table
+  taxo = as.data.frame(dataT)
+  CheckTaxo = CheckTaxoTable(taxo,counts)
+  
+  ## Pourcentage of annotation
+  tmp = PercentAnnot(counts,taxo)
+  
+  return(list(counts=counts,taxo=taxo,CheckCounts=CheckCounts,CheckTaxo=CheckTaxo,Percent=tmp$Percent,CheckPercent=tmp$Error))
+}
+
+
+
+## Get the counts for the selected taxonomy
+GetCountsMerge <- function(input,dataInput,taxoSelect,target,design)
+{
+  ## Init
+  counts= NULL
+  CheckTarget = FALSE
+  CT_noNorm = NULL
+  normFactors = NULL
+  FeatureSize = NULL
+  
+  ## Counts and taxo tables
+  CT = dataInput$counts
+  taxo = dataInput$taxo
+  
+  ## Select cols in the target
+  labels = target[,1]
+  ind = which(colnames(CT)%in%labels)
+  
+  ## Get the normalization variable (normalization can be done according to this variable)
+  VarNorm = input$VarNorm
+  
+  if(length(ind)==length(labels))
+  { 
+    if(input$TypeTable == "MGS"){
+      ## Get the feature size for the normalisation
+      Size_indcol = which(toupper(colnames(CT))%in%"SIZE")
+      if(length(Size_indcol)==1) FeatureSize = CT[,Size_indcol]
+      else print("Size parameter is missing in the count matrix")
+      # Consider only counts
+      CT = CT[,ind]
+      # Divide by gene length
+      CT = CT / FeatureSize * 1000
+      # Convert matrix as integer
+      CT_int=t(apply(CT,1,as.integer))
+      rownames(CT_int)=rownames(CT)
+      colnames(CT_int)=colnames(CT)
+      CT=CT_int
+    } else CT = CT[,ind]
+    
+    ## Order CT according to the target
+    CT = OrderCounts(counts=CT,labels=labels)$CountsOrder
+    CT_noNorm = CT
+    RowProd = sum(apply(CT_noNorm,1,prod))
+    
+    ## Create the dds object
+    dds <- DESeqDataSetFromMatrix(countData=CT, colData=target, design=design)
+    
+    if(is.null(VarNorm)){
+      ## Counts normalisation
+      ## 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)))
+      normFactors = sizeFactors(dds)
+      
+    } else{
+      group = as.data.frame(target[,VarNorm])
+      group = apply(group,1,paste, collapse = "-")
+      normFactors = c()
+      mod = unique(group)
+      ## At least 2 samples are needed for the normalization
+      if(min(table(group))>1){
+        for(i in unique(group))
+        {
+          indgrp = which(group==i) 
+          CT_tmp = CT[,indgrp]
+          CT_tmp = removeNulCounts(CT_tmp) 
+          target_tmp = data.frame(labels = rownames(target)[indgrp])
+          dds_tmp <- DESeqDataSetFromMatrix(countData=CT_tmp, colData=target_tmp, design=~labels)
+          if(input$AccountForNA) dds_tmp = estimateSizeFactors(dds_tmp,locfunc=eval(as.name(input$locfunc)),geoMeans=GeoMeansCT(CT_tmp))
+          if(!input$AccountForNA) dds_tmp = estimateSizeFactors(dds_tmp,locfunc=eval(as.name(input$locfunc)))
+          normFactors[indgrp] = sizeFactors(dds_tmp)
+        }
+      } else{
+        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)))
+        normFactors = sizeFactors(dds)
+      }
+      
+      sizeFactors(dds) = normFactors
+    }
+    
+    ## Keep normalized OTU table
+    CT_Norm = counts(dds, normalized=TRUE)
+    
+    # Only interesting OTU
+    merged_table = merge(CT, taxo[order(rownames(CT)),], by="row.names")
+    CT = merged_table[,2: (dim(CT)[2]+1)]
+    taxo = merged_table[,(dim(CT)[2]+2):dim(merged_table)[2]]
+    rownames(CT) = merged_table[,1]
+    rownames(taxo) = merged_table[,1]
+    #ordOTU = order(rownames(taxo))
+    counts_annot = CT
+    #       ordOTU = order(rownames(taxo))
+    #       indOTU_annot = which(rownames(CT)%in%rownames(taxo))
+    #       counts_annot = CT[indOTU_annot[ordOTU],]
+    ## Aggregate matrix
+    if(taxoSelect=="OTU/Gene") counts = counts_annot
+    else{
+      if(input$TypeTable == "MGS"){
+        taxoS = taxo[,input$TypeTable]
+        counts = aggregate(counts_annot,by=list(Taxonomy = taxoS),mean)
+        rownames(counts)=counts[,1]
+        counts=counts[,-1]
+        counts_int=t(apply(counts,1,as.integer))
+        rownames(counts_int)=rownames(counts)
+        colnames(counts_int)=colnames(counts)
+        counts=counts_int
+      }
+      if(taxoSelect != "MGS"){
+        #taxoS = taxo[ordOTU,taxoSelect]
+        taxoS = taxo[,taxoSelect]
+        counts = aggregate(counts_annot,by=list(Taxonomy = taxoS),sum)
+        rownames(counts)=counts[,1];counts=counts[,-1]
+      }
+    }
+    
+    ## Ordering the counts table according to the target labels 
+    tmpOrder = OrderCounts(counts,normFactors,labels)
+    counts = tmpOrder$CountsOrder
+    normFactors = tmpOrder$normFactorsOrder
+    CheckTarget = TRUE
+  }
+  return(list(counts=counts,CheckTarget=CheckTarget,normFactors=normFactors, CT_noNorm=CT_noNorm, CT_Norm =CT_Norm))
+  #return(list(counts=counts,target=target[ind,],labeled=labeled,normFactors=normFactors, CT_noNorm=CT_noNorm))
+}
+
+
+## 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)
+}
+
+
+## Order the counts 
+OrderCounts <- function(counts,normFactors=NULL,labels)
+{
+  n = length(labels)
+  CountsOrder = counts
+  normFactorsOrder = normFactors
+  for(i in 1:n)
+  {
+    ind = which(labels[i]==colnames(counts))
+    CountsOrder[,i] = counts[,ind]
+    if(!is.null(normFactors))