internal.R

## 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)
}


## 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,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))
  }
  
  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
    counts = biom_data(dataBIOM)
    counts = as.matrix(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))
  }
  
  
  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))
  }
  
  GetInteraction2 <- function(target)
  { 
    res=c()
    namesTarget = colnames(target)[2:ncol(target)]
    k=1
    for(i in 1:(length(namesTarget)-1))
    { 
      for(j in (i+1):length(namesTarget))
      { 
        res[k] = paste(namesTarget[i],":",namesTarget[j],sep="")
        k = k+1
      }
    }
    
    return(res)
  }
  


  ## Print the contrasts
  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)
      out = paste(out,paste("<b>",contnames[i], "</b> <br/>", paste(contrast[contrast != 0], coefs[contrast != 0], collapse = " ", sep = " ")),"<br/>")
    }
    return(out)
    
  }

  ## 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)) normFactorsOrder[i] = normFactors[ind]
    }
    colnames(CountsOrder) = labels
    return(list(CountsOrder=CountsOrder,normFactorsOrder = normFactorsOrder))
  }
  
  
  ## Get the dds object of DESeq2
  Get_dds_object <- function(input,counts,target,design,normFactorsOTU,CT_noNorm,CT_Norm)
  {
    dds <- DESeqDataSetFromMatrix(countData=counts, colData=target, design=design)
    sizeFactors(dds) = normFactorsOTU
    dds <- estimateDispersions(dds, fitType=input$fitType)
    if(as.numeric(R.Version()$major)+as.numeric(R.Version()$minor) >= 4.3){
      dds <- nbinomWaldTest(dds)
    }else{
      dds <- nbinomWaldTest(dds,modelMatrixType = "expanded")
    }
    countsNorm = counts(dds, normalized = TRUE)
    
    save(dds,file="dds.RData")
    return(list(dds = dds,raw_counts=counts,countsNorm=countsNorm,target=target,design=design,normFactors = normFactorsOTU,CT_noNorm=CT_noNorm,CT_Norm=CT_Norm))
  }

  ## 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)
  {
    
    VarInt = input$VarInt
    dds = resDiff$dds
    counts = resDiff$raw_counts
    if(input$CountsType=="Normalized") counts = resDiff$countsNorm
    target = resDiff$target
    normFactors = resDiff$normFactors
    
    ## Counts at the OTU level
    CT = resDiff$CT_noNorm
    if(input$CountsType=="Normalized") CT = resDiff$CT_Norm
    
    group = as.data.frame(target[,VarInt])
    rownames(group) = rownames(target)
    res = NULL
    
    if(ncol(group)>0 && nrow(counts)>0)
    { 
      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))
      
      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)
      if(input$DiagPlot=="boxplotNorm") res = boxplotNorm(input,CT,group = group, col=colors)
      if(input$DiagPlot=="DispPlot") res = plotDispEsts(dds)
      if(input$DiagPlot=="MajTax") res = majTaxPlot(input,counts, group = group, col=colors)
      if(input$DiagPlot=="SfactorsVStot") res = diagSFactors(input,normFactors,resDiff$raw_counts) 
      if(input$DiagPlot=="pcaPlot") res = PCAPlot_meta(input,dds, group,  type.trans = input$TransType, col = colors)
      if(input$DiagPlot=="pcoaPlot") res = PCoAPlot_meta(input,dds, group, col = colors) 
      if(input$DiagPlot=="clustPlot") res = HCPlot(input,dds,group,type.trans=input$TransType,counts,col=colors)
    }
    
    return(res)
  }

  
  HCPlot <- function (input,dds,group,type.trans,counts,col = c("lightblue", "orange", "MediumVioletRed", "SpringGreen")) 
  {
    
    res = NULL
    
    ## Get the counts
    if (input$DistClust == "euclidean" && type.trans == "VST") counts <- assay(varianceStabilizingTransformation(dds))
    if (input$DistClust == "euclidean" && type.trans == "rlog") counts <- assay(rlogTransformation(dds))
    
    ## Get the group of leaf
    group = apply(group,1,paste, collapse = "-")    
    nb = length(unique((group)))
    
    ## Get the dendrogram
    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")
    
    dend = as.dendrogram(hc)
    
    ## Get the type of dendrogram
    type <- input$typeHculst
    
    dend <- set(dend, "labels_cex", input$cexLabelDiag)
    if(input$colorHC) labels_colors(dend)<-col[as.integer(as.factor(group))][order.dendrogram(dend)]
    if(type=="hori") 
    { 
      par(cex=input$cexTitleDiag,mar=c(6,6,4,5))
      res = plot(dend, main = "Cluster dendrogram",xlab = paste(input$DistClust,"distance, Ward criterion",sep=" "),cex=input$cexLabelDiag)
    }  
    if(type!="hori")
    {
      par(cex=input$cexTitleDiag,mar=c(6,6,4,5))
      res = circlize_dendrogram(dend, labels_track_height = 0.2, dend_track_height = .3, main = "Cluster dendrogram",xlab = paste(input$DistClust,"distance, Ward criterion",sep=" "))
    }
    return(res)
  }
  
  
  ## Color for the horizontal dendro
  colLabdendo <- function(n,group) {
    
    group = apply(group,1,paste, collapse = "-")
    
    nb = length(unique((group)))
    namesGrp = names(group)

    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") 
  }
  
  
  my.boxplot <- function(x, pol.col = 1, pol.density = NULL, pol.angle = 45,
                         bxp.pars = list(boxwex = 0.8, staplewex = 0.5, outwex = 0.5), ...){
    res <- boxplot(x, pars = bxp.pars, ...) # que boxplot se d�merde avec ses arguments
    n <- ncol(res$stats) # nombre de boxplots
    density <- if(is.null(pol.density)){NULL}else{rep(pol.density, length = n)}
    angle <- if(is.null(pol.angle)){NULL}else{rep(pol.angle, length = n)}
    col <- if(is.null(pol.col)){NULL}else{rep(pol.col, length = n)}
    # Ajout des textures
    ex <- bxp.pars$boxwex/2 # j'ai juste besoin de la largeur des bo�tes pass�e � bxp
    for(i in 1:n){
      polygon(c(i - ex, i - ex, i + ex, i + ex),
              c(res$stats[2, i], res$stats[4, i], res$stats[4, i], res$stats[2, i]),
              density = density[i], angle = angle[i], col = col[i])
      segments(i-ex,res$stats[3,i],i+ex,res$stats[3,i],lwd=3,col="black",lend=1)
    }
  }
  
  
  
  ## Boxplot for the counts normalized/no normalized
  boxplotNorm <- function(input,CT, group, col = c("lightblue","orange", "MediumVioletRed", "SpringGreen")) 
  {
    
    ncol1 <- ncol(group) == 1
    par(cex=input$cexTitleDiag,mar=c(12,6,4,5))
    if(input$RemoveNullValue) CT[CT==0] = NA
    
    ### Boxplots of the counts
    my.boxplot(log2(CT+1), las = 2, pol.col = col[as.integer(group[,1])],
               pol.density = if (ncol1) {NULL}
               else {15}, 
               pol.angle = if (ncol1) {NULL}
            else {seq(0,160,length.out =nlevels(group[, 2]))[as.integer(group[, 2])]},
            main = paste(input$CountsType, "counts distribution"), ylab = expression(log[2] ~ ( count + 1)))
    legend("topright", levels(group[, 1]), fill = col[1:nlevels(group[,1])], bty = "n")
    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")
    
     }
  

  

  ## barplot total
  barplotTot <- function(input,counts, group, cex.names = 1, col = c("lightblue","orange", "MediumVioletRed", "SpringGreen")) 
  {

    ncol1 <- ncol(group) == 1
    par(cex=input$cexTitleDiag,mar=c(12,6,4,5))
    barplot(colSums(counts), cex.names = cex.names, main = "Total mapped read count per sample", ylab = "Total mapped read count", 
            ylim = c(0, max(colSums(counts)) * 1.2), density = if (ncol1) {NULL}
            else {15}, 
            angle = if (ncol1) {NULL}
            else {seq(0,160,length.out =nlevels(group[, 2]))[as.integer(group[, 2])]}, col = col[as.integer(group[, 1])], las = 2)
    legend("topright", levels(group[, 1]), fill = col[1:nlevels(group[,1])], bty = "n")
    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")
  
  }


  ## barplot Nul 
  barPlotNul <-function (input,counts, group, cex.names = 1, col = c("lightblue","orange", "MediumVioletRed", "SpringGreen")) 
  {
    
    percentage <- apply(counts, 2, function(x) {sum(x == 0)}) * 100/nrow(counts)
    percentage.allNull <- (nrow(counts) - nrow(removeNulCounts(counts))) * 100/nrow(counts)
    ncol1 <- ncol(group) == 1
    
    par(cex=input$cexTitleDiag,mar=c(12,6,4,5))

    barplot(percentage, las = 2, col = col[as.integer(group[,1])], 
            density = if (ncol1) {NULL}
            else {15}, 
            angle = if (ncol1) {NULL}
            else {seq(0,160,length.out =nlevels(group[, 2]))[as.integer(group[, 2])]},
            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")
    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")
  }


  ## Plot density
  densityPlotTot <-function (input,counts, group, col = c("lightblue", "orange", "MediumVioletRed", "SpringGreen")) 
  {
    
    counts <- removeNulCounts(counts)
    ncol1 <- ncol(group) == 1
    par(cex=input$cexTitleDiag,mar=c(12,5,4,5))
    plot(density(log2(counts[, 1] + 1)), las = 1, lwd = 2, main = "Density of counts distribution", 
         xlab = expression(log[2] ~ (raw ~ count + 1)), 
         ylim = c(0, max(apply(counts, 2, function(x) {max(density(log2(x + 1))$y)})) * 1.05), 
         lty = if (ncol1) {1}
         else{rep(seq(1:6),ceiling(nlevels(group[, 2])/6))[as.integer(group[, 2])[1]]}, 
         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}
            else{rep(seq(1:6),ceiling(nlevels(group[, 2])/6))[as.integer(group[, 2])[i]]})
    }
    legend("topright", levels(group[, 1]), lty = 1, col = col[1:nlevels(group[,1])], lwd = 2, bty = "n")
    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")
    
  }


  ## Table of maj. taxo
  majTab <- function(input,counts,n)
  {
    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
  majTaxPlot <-function (input,counts, n = 3, group, cex.names = 1, col = c("lightblue",  "orange", "MediumVioletRed", "SpringGreen")) 
  {
    p = majTab(input,counts,n)
    maj <- apply(p, 2, max)
    seqname <- rownames(p)[apply(p, 2, which.max)]
    ncol1 <- ncol(group) == 1
    par(cex=input$cexTitleDiag,mar=c(12,6,4,5))
    x <- barplot(maj, col = col[as.integer(group[, 1])], main = "Proportion of mapped reads from\nmost expressed sequence",
                 ylim = c(0, max(maj) * 1.2), cex.main = 1, 
                 cex.names = cex.names, las = 2, ylab = "Proportion of mapped reads", 
                 density = if (ncol1) {NULL}
                 else {15}, 
                 angle = if (ncol1) {NULL}
                 else {seq(0,160,length.out =nlevels(group[, 2]))[as.integer(group[, 2])]})
    
    legend("topright", levels(group[, 1]), fill = col[1:nlevels(group[,1])], bty = "n")
    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")
    
    for (i in 1:length(seqname)) text(x[i], maj[i]/2, seqname[i], cex=input$cexLabelDiag, srt = 90, adj = 0)
  }
  

  ## plot SERE Coefs
#   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")
#     
#   }
  
  
  ## Get the SERE COEF
  SEREcoef<-function(counts)
  {
    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)) {
        sere[i, j] <- sigfun_Pearson_meta(counts[, c(i, j)])
      }
    }
    sere=sere+t(sere)
    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)
  }


  ## Plots of size factors
  diagSFactors<-function (input,normFactors,counts) 
  {
    geomeans <- exp(rowMeans(log(counts)))
    samples <- colnames(counts)
      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),
           ylab = "Total number of reads", xlab = "Size factors", 
           main = "Diagnostic: size factors vs total number of reads")
      if(input$addLabelSFact) text(normFactors,colSums(counts),labels = samples,cex=input$cexLabelDiag)
      abline(lm(colSums(counts) ~ normFactors + 0), lty = 2, col = "grey")
  }

  
  ### PCoA
  PCoAPlot_meta <-function (input, dds, group_init, col = c("SpringGreen","dodgerblue","black","firebrick1"), plot = "pcoa") 
  {
    cval=c()
    time_set = 0
    # Set of shape
    shape=c(19,17,15,18)
    
    ## Var of interest
    VarInt  = input$VarInt
    
    ## 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)
    
    ## Get the group corresponding to the modalities
    group = group[ind_kept]
    nb = length(unique((group)))
    group = as.factor(group)
    
    if(nlevels(group)!=0)
    { 
      ## Get the norm data
      counts.norm = as.data.frame(round(counts(dds)))
      if(input$CountsType=="Normalized") counts.norm = as.data.frame(round(counts(dds, normalized = TRUE)))
      # was removed
      counts.norm = counts.norm[,ind_kept]
  
      ## Get the distance
      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))
      
      ## 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))
      }
      #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"){
        par(cex=input$cexTitleDiag,mar=c(6,6,4,5))
        ## 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),"%"),
             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 ')
        # 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),
                        add.plot = TRUE, cpoint = input$cexpoint, cell=input$cexcircle, clabel=input$cexLabelDiag,  cstar = input$cexstar)
        }  
        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)
      }
  }
  
  }
  
  ### PCA
  PCAPlot_meta <-function(input,dds, group_init, n = min(500, nrow(counts(dds))), type.trans = c("VST", "rlog"), 
                           col = c("lightblue", "orange", "MediumVioletRed", "SpringGreen"),plot="pca") 
  {
    ## Var of interest
    VarInt  = input$VarInt
    
    group = as.character(apply(group_init,1,paste, collapse = "-"))
    group_init = group
    
    ## 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]
    
    
    if(nlevels(group)!=0)
    { 
      type.trans <- type.trans[1]
      
      if (type.trans == "VST") counts.trans <- assay(varianceStabilizingTransformation(dds))
      else counts.trans <- assay(rlogTransformation(dds))
      counts.trans = counts.trans[,ind_kept]
      
      rv = apply(counts.trans, 1, var, na.rm = TRUE)
      pca = prcomp(t(counts.trans[order(rv, decreasing = TRUE),][1:n, ]))
      
      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)}
    }
  }
  
  
  

  ############################################################
  ##
  ##              CREATE THE CONTRAST DATABASE
  ##
  ############################################################

  
  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)
  }
  
  
  BaseContrastEasy <- function(input,names,namesfile,target)
  {  
    
    v_tmp = rep(0,length(names))
    print(names)
    filesize = file.info(namesfile)[,"size"]
    F1 = NULL
    nameContrast = ""
    
    ## Get the selected modalities
    M1 = input$Select1_contrast
    M2 = input$Select2_contrast
    print(M1)
    print(M2)
    
    if(length(input$Interaction2)>0) F1 = input$Select3_contrast
    ## Get the name of the parameter corresponding to the modalities
    InterVar = input$InterestVar
    Sel_Var = InterVar[which(unlist(lapply(target[,InterVar],FUN = function(x){M1%in%x})))]
    names1dds = paste(Sel_Var,M1,sep="")
    names2dds = paste(Sel_Var,M2,sep="")
    
    ## fill the vector
    ind1 = which(names%in%names1dds)
    ind2 = which(names%in%names2dds)
    print(ind1)
    print(ind2)
    if(length(ind1)>0) v_tmp[ind1] = 1
    if(length(ind2)>0) v_tmp[ind2] = -1
    print(v_tmp)
    
    nameContrast = paste(M1,"_vs_",M2,sep="")
    
    if(F1!="All" && !is.null(F1)){
      Sel_Var_For = InterVar[which(unlist(lapply(target[,InterVar],FUN = function(x){F1%in%x})))]
      ## Depends on the interation
      namesfor1 = paste(Sel_Var,M1,".",Sel_Var_For,F1,sep="")
      namesfor1.1 = paste(Sel_Var_For,F1,".",Sel_Var,M1,sep="")
      namesfor2 = paste(Sel_Var,M2,".",Sel_Var_For,F1,sep="")
      namesfor2.1 = paste(Sel_Var_For,F1,".",Sel_Var,M2,sep="")
      ind1.for = which(names%in%c(namesfor1,namesfor1.1))
      ind2.for = which(names%in%c(namesfor2,namesfor2.1))
      if(length(ind1.for)>0) v_tmp[ind1.for] = 1
      if(length(ind2.for)>0) v_tmp[ind2.for] = -1
      nameContrast = paste(nameContrast,"_for_",F1,sep="")
    }
    print(v_tmp)
    
    if(filesize!=0)
    { 
      oldContrast = read.table(namesfile,header=TRUE)
      print(oldContrast)
      colnamesTmp = c(colnames(oldContrast),nameContrast)
      mat = cbind(oldContrast,v_tmp)
    }
    else{ colnamesTmp = nameContrast; mat = v_tmp}
    
    write.table(mat,namesfile,row.names=FALSE,col.names = colnamesTmp)
  }
  
  
  ## Remove nul counts
  removeNulCounts <-function (counts) 
  {
    return(counts[rowSums(counts) > 0, ])
  }

  
  ############################################################
  ##
  ##              VISUALISATION PLOTS
  ##
  ############################################################