server.R

library(shinydashboard)

shinyServer
(function(input, output, session)
{
  checkNull <- function(x) {
    if (!exists(as.character(substitute(x)))) {
      return(NULL)
    } else if (is.null(x)) {
      return(NULL)
    } else if (length(x) > 1) {
      return(x)
    }
    else if (x %in% c(0, "", NA, "NULL")) {
      return(NULL)
    } else {
      return(x)
    }
  }
  
  beautifulTable <- function(data)  {
    DT::datatable(
      data = data,
      rownames = FALSE,
      filter = "top",
      extensions = c("Buttons", "ColReorder", "FixedColumns"),
      options = list(
        dom = "lBtip",
        pageLength = 10,
        lengthMenu = list(c(10, 25, 50, 100, -1), list('10', '25', '50', '100', 'All')),
        buttons = list(
          'colvis',
          list(
            extend = 'collection',
            buttons = c('copy', 'csv', 'excel', 'pdf'),
            text = 'Download'
          )
        ),
        colReorder = TRUE,
        scrollX = TRUE,
        fixedColumns = list(leftColumns = 1, rightColumns = 0)
      ),
      width = "auto",
      height = "auto"
    )
  }
  
  source(
    "https://raw.githubusercontent.com/mahendra-mariadassou/phyloseq-extended/master/R/load-extra-functions.R"
  )
  
  data16S <- reactive({
    if (input$demo != "input")
    {
      load("demo/demo.RData")
      return(get(input$demo))
    }
    else {
      if (is.null(input$fileBiom))
      {
        return()
      }
      if (input$biomFormat == "std")
      {
        d <- import_biom(
          BIOMfilename = input$fileBiom$datapath,
          treefilename = input$fileTree$datapath# ,
          # refseqfilename = input$fileSeq$datapath
        )
      } else if (input$biomFormat == "frogs") {
        d <- import_frogs(
          biom = input$fileBiom$datapath,
          treefilename = input$fileTree$datapath# ,
          # refseqfilename = input$fileSeq$datapath
        )
      }
      
      colnames(tax_table(d)) <-
        c("Kingdom",
          "Phylum",
          "Class",
          "Order",
          "Family",
          "Genus",
          "Species",
          "Strain")[1:length(rank_names(d))]
      tax_table(d)[grep("unknown ", tax_table(d))] <- NA
      #tax_table(d)[grep("Unclassified", tax_table(d))] <- NA
      if (!is.null(input$fileMeta)) {
        if (input$CSVsep == "excel") {
          sample_data(d) <-
            RcmdrMisc::readXL(input$fileMeta$datapath,
                              rownames = TRUE,
                              header = TRUE)
        } else {
          sample_data(d) <- read.csv(
            input$fileMeta$datapath,
            header = TRUE,
            sep = input$CSVsep,
            row.names = 1,
            na.strings = NA
          )
        }
      } else {
        n <- data.frame(sample_names(d) , row.names = sample_names(d))
        names(n) <- "sample_names"
        sample_data(d) <- n
      }
      if (input$rareData) {
        d <- rarefy_even_depth(
          d,
          replace = FALSE,
          rngseed = as.integer(Sys.time()),
          verbose = FALSE
        )
      }
      return(d)
    }
  })
  
  output$rarefactionMin <- renderText({
    if (!is.null(input$fileBiom)) {
      paste("(min sample =", format(min(sample_sums(data16S(
        
      ))), big.mark = " "), "reads)")
    } else {
      paste("(min sample =", 0, "reads)")
    }
  })
  
  output$phyloseqPrint <- renderPrint({
    validate(
      need(
        data16S(),
        "Merci de commencer par importer un fichier d'abondance au format BIOM. Celui-ci peut etre obtenu a l'issue du workflow FROGS avec l'operation 'FROGS BIOM to std BIOM'.\nVous pouvez egalement choisir un dataset de demo"
      )
    )
    data16S()
  })
  
  output$summaryTable <- renderUI({
    validate(need(data16S(), ""))
    box(
      title = "Tables",
      width = NULL,
      status = "primary",
      tabsetPanel(
        tabPanel("otu_table",
                 beautifulTable(
                   data.frame(OTU = taxa_names(data16S()), otu_table(data16S()))
                 )),
        tabPanel("tax_table",
                 beautifulTable(
                   data.frame(OTU = taxa_names(data16S()), tax_table(data16S()))
                 )),
        tabPanel("sample_data",
                 beautifulTable(
                   data.frame(SAMPLE = sample_names(data16S()), sample_data(data16S()))
                 ))
      )
    )
  })
  
  output$histUI <- renderUI({
    validate(need(data16S(), ""))
    box(
      title = "Paramètres",
      width = NULL,
      status = "primary",
      selectInput(
        "barFill",
        label = "Niveau taxo :",
        choices = rank_names(data16S())
      ),
      selectInput(
        "barGrid",
        label = "Regroupement :",
        choices = c("..." = 0, sample_variables(data16S()))
      )
      ,
      selectInput(
        "barX",
        label = "X :",
        choices = c("..." = 0, sample_variables(data16S()))
      )
    )
  })
  
  output$histo <- renderPlot({
    validate(need(data16S(),
                  "Merci d'importer un fichier d'abondance"))
    p <- plot_bar(
      physeq = data16S(),
      fill = input$barFill,
      x = ifelse(is.null(checkNull(input$barX)), "Sample", input$barX)
    )
    if (!is.null(checkNull(input$barGrid))) {
      p <-
        p + facet_grid(paste(".", "~", input$barGrid), scales = "free_x")
    }
    return(p)
  })
  
  output$histFocusUIfocusRank <- renderUI({
    validate(need(data16S(), ""))
    radioButtons(
      "focusRank",
      label = "Niveau taxo :",
      choices = rank_names(data16S())[-length(rank_names(data16S()))],
      inline = TRUE
    )
  })
  
  output$histFocusUIfocusTaxa <- renderUI({
    validate(need(data16S(), ""))
    selectInput(
      "focusTaxa",
      label = "Taxa :",
      choices = unique(as.vector(tax_table(data16S(
        
      ))[, input$focusRank])),
      selected = TRUE
    )
  })
  
  output$histFocusUIfocusNbTaxa <- renderUI({
    validate(need(data16S(), ""))
    sliderInput(
      "focusNbTaxa",
      label = "Nombre de sous-taxons :",
      min = 0,
      #max = sum(tax_table(tax_glom(data16S(), rank_names(data16S())[1+as.integer(input$focusRank)]))[, as.integer(input$focusRank)]==input$focusTaxa)
      max = 30
      ,
      value = 10
    )
  })
  
  output$histFocusUIfocusGrid <- renderUI({
    validate(need(data16S(), ""))
    selectInput("focusGrid",
                label = "Regroupement :",
                choices = c("..." = 0, sample_variables(data16S())))
  })
  
  output$histFocusUIfocusX <- renderUI({
    validate(need(data16S(), ""))
    selectInput("focusX",
                label = "X :",
                choices = c("..." = 0, sample_variables(data16S())))
  })
  
  output$histoFocus <- renderPlot({
    validate(need(data16S(),
                  "Merci d'importer un fichier d'abondance"))
    p <- plot_composition(
      physeq = data16S(),
      taxaRank1 = input$focusRank,
      taxaSet1 = input$focusTaxa,
      taxaRank2 = rank_names(data16S())[which(rank_names(data16S()) == input$focusRank) + 1],
      numberOfTaxa = input$focusNbTaxa,
      fill = rank_names(data16S())[which(rank_names(data16S()) == input$focusRank) + 1],
      x = ifelse(is.null(checkNull(input$focusX)), "Sample", input$focusX)
    )
    if (!is.null(checkNull(input$focusGrid))) {
      p <-
        p + facet_grid(paste(".", "~", input$focusGrid), scales = "free_x")
    }
    return(p)
  })
  
  output$clustUI <- renderUI({
    validate(need(data16S(), ""))
    box(
      title = "Paramètres",
      width = NULL,
      status = "primary",
      selectInput(
        "clustDist",
        label = "Distance :",
        choices = list(
          "bray",
          "jaccard",
          "unifrac",
          "wunifrac",
          "dpcoa",
          "jsd",
          "euclidean"
        )
      ),
      selectInput(
        "clustMethod",
        label = "Methode :",
        choices = list(
          "ward.D2",
          "ward.D",
          "single",
          "complete",
          "average",
          "mcquitty",
          "median",
          "centroid"
        )
      ),
      selectInput(
        "clustCol",
        label = "Couleur :",
        choices = c("..." = 0, sample_variables(data16S()))
      )
    )
  })
  
  output$clust <- renderPlot({
    validate(need(data16S(),
                  "Merci d'importer un fichier d'abondance"))
    plot_clust(
      physeq = data16S(),
      dist = input$clustDist,
      method = input$clustMethod,
      color = checkNull(input$clustCol)
    )
  })
  
  output$richnessAUI <- renderUI({
    validate(need(data16S(), ""))
    box(
      title = "Paramètres",
      width = NULL,
      status = "primary",
      checkboxGroupInput(
        "richnessMeasures",
        label = "Mesures :",
        choices = c(
          "Observed",
          "Chao1",
          "ACE",
          "Shannon",
          "Simpson",
          "InvSimpson",
          "Fisher"
        ),
        selected = c(
          "Observed",
          "Chao1",
          "ACE",
          "Shannon",
          "Simpson",
          "InvSimpson",
          "Fisher"
        ),
        inline = TRUE
      ),
      selectInput(
        "richnessX",
        label = "X :",
        choices = c("..." = 0, sample_variables(data16S()))
      ),
      selectInput(
        "richnessColor",
        label = "Couleur :",
        choices = c("..." = 0, sample_variables(data16S()))
      ),
      selectInput(
        "richnessShape",
        label = "Forme :",
        choices = c("..." = 0, sample_variables(data16S()))
      ),
      radioButtons(
        "richnessBoxplot",
        label = "Representation :",
        choices = list(
          "Points seuls" = 1,
          "Boxplot et points" = 2,
          "Boxplot seul" = 3
        ),
        selected = 2,
        inline = TRUE
      )
    )
  })
  
  output$richnessA <- renderPlot({
    validate(need(data16S(),
                  "Merci d'importer un fichier d'abondance"))
    p <- plot_richness(
      physeq = data16S(),
      x = ifelse(is.null(checkNull(
        input$richnessX
      )), "samples", input$richnessX),
      color = checkNull(input$richnessColor),
      shape = checkNull(input$richnessShape),
      measures = checkNull(input$richnessMeasures)
    )
    if (input$richnessBoxplot >= 2) {
      p <- p + geom_boxplot()
    }
    if (input$richnessBoxplot <= 2) {
      p <- p + geom_point()
    }
    return(p)
  })
  
  output$richnessATable <- renderUI({
    validate(need(data16S(),
                  "Merci d'importer un fichier d'abondance"))
    p(beautifulTable(data.frame(
      SAMPLE = sample_names(data16S()), round(estimate_richness(data16S()), digits = 2)
    )))
  })
  
  output$richnessBUI <- renderUI({
    box(
      title = "Paramètres",
      width = NULL,
      status = "primary",
      selectInput(
        "richnessOrder",
        label = "Ordre de tri des echantillons :",
        choices = c("..." = 0, sample_variables(data16S()))
      )
    )
  })
  
  output$richnessB <- renderPlot({
    validate(need(data16S(),
                  "Merci d'importer un fichier d'abondance"))
    beta <-
      melt(as(distance(data16S(), method = input$richnessBDist), "matrix"))
    colnames(beta) <- c("x", "y", "distance")
    if (!is.null(checkNull(input$richnessOrder)))
    {
      new_factor = as.factor(get_variable(data16S(), input$richnessOrder))
      variable_sort <-
        as.factor(get_variable(data16S(), input$richnessOrder)[order(new_factor)])
      L = levels(reorder(sample_names(data16S()), as.numeric(new_factor)))
      beta$x <- factor(beta$x, levels = L)
      beta$y <- factor(beta$y, levels = L)
      palette <- hue_pal()(length(levels(new_factor)))
      tipColor <-
        col_factor(palette, levels = levels(new_factor))(variable_sort)
    }
    p <-
      ggplot(beta, aes(x = x, y = y, fill = distance)) + geom_tile()
    p <- p + theme(
      axis.text.x = element_text(
        angle = 90,
        hjust = 1,
        color = checkNull(tipColor)
      ),
      axis.text.y = element_text(color = checkNull(tipColor)),
      axis.title.x = element_blank(),
      axis.title.y = element_blank()
    )
    return(p + scale_fill_gradient2())
  })
  
  output$networkBUI <- renderUI({
    validate(need(data16S(), ""))
    box(
      title = "Paramètres",
      width = NULL,
      status = "primary",
      sliderInput(
        "netwMax",
        label = "Cutoff :",
        min = 0,
        max = 1,
        value = 0.7
      ),
      checkboxInput("netwOrphan",
                    label = "Garder les points orphelins",
                    value = TRUE),
      selectInput(
        "netwCol",
        label = "Couleur :",
        choices = c("..." = 0, sample_variables(data16S()))
      ),
      selectInput(
        "netwShape",
        label = "Forme :",
        choices = c("..." = 0, sample_variables(data16S()))
      ),
      selectInput(
        "netwLabel",
        label = "Label :",
        choices = c(
          "..." = 0,
          "Sample name" = "value",
          sample_variables(data16S())
        )
      )
    )
  })
  
  output$networkB <- renderPlot({
    validate(need(data16S(),
                  "Merci d'importer un fichier d'abondance"))
    g <- make_network(
      data16S(),
      distance = input$richnessBDist,
      max.dist = input$netwMax,
      keep.isolates = input$netwOrphan
    )
    p <- plot_network(
      g,
      physeq = data16S(),
      color = checkNull(input$netwCol),
      shape = checkNull(input$netwShape),
      label = checkNull(input$netwLabel),
      hjust = 2,
      title = NULL
    )
    return(p)
  })
  
  output$richnessBTable <- renderUI({
    validate(need(data16S(),
                  "Merci d'importer un fichier d'abondance"))
    p(beautifulTable(data.frame(
      SAMPLE = sample_names(data16S()), round(as.matrix(
        distance(data16S(), method = input$richnessBDist)
      ), digits = 2)
    )))
  })
  
  output$rarefactionCurve <- renderPlot({
    validate(need(data16S(),
                  "Merci d'importer un fichier d'abondance"))
    p <- ggrare(
      physeq = data16S(),
      step = 100,
      #step = input$rarefactionStep,
      color = checkNull(input$rarefactionColor),
      se = FALSE
    )
    if (!is.null(checkNull(input$rarefactionGrid))) {
      p <- p + facet_grid(paste(".", "~", input$rarefactionGrid))
    }
    
    if (input$rarefactionMin) {
      p <-
        p + geom_vline(xintercept = min(sample_sums(data16S())),
                       color = "gray60")
    }
    return(p)
    
  })
  
  output$rarefactionCurveUI <- renderUI({
    validate(need(data16S(), ""))
    box(
      title = "Paramètres",
      width = NULL,
      status = "primary",
      # sliderInput(
      #   "rarefactionStep",
      #   label = "Etapes de calcul :",
      #   min = 1,
      #   max = 1000,
      #   value = 100
      # ),
      checkboxInput("rarefactionMin", label = "Afficher le seuil de l'echantillon minimal", value = FALSE),
      selectInput(
        "rarefactionColor",
        label = "Couleur :",
        choices = c("..." = 0, sample_variables(data16S()))
      ),
      selectInput(
        "rarefactionGrid",
        label = "Regroupement :",
        choices = c("..." = 0, sample_variables(data16S()))
      )
    )
  })
  
  output$HeatmapUI <- renderUI({
    validate(need(data16S(), ""))
    box(
      title = "Paramètres",
      width = NULL,
      status = "primary",
      selectInput(
        "heatmapGrid",
        label = "Regroupement :",
        choices = c("..." = 0, sample_variables(data16S()))
      ),
      selectInput(
        "heatmapX",
        label = "X :",
        choices = c("..." = 0, sample_variables(data16S()))
      ),
      sliderInput(
        "heatmapTopOtu",
        label = "Selection des n OTU les plus abondant :",
        min = 1,
        max = ntaxa(data16S()),
        value = 250
      ),
      selectInput(
        "heatmapDist",
        label = "Distance :",
        selected = "bray",
        choices = list(
          "bray",
          "jaccard",
          "unifrac",
          "wunifrac",
          "dpcoa",
          "jsd",
          "euclidean"
        )
      ),
      selectInput(
        "heatmapMethod",
        label = "Methode :",
        selected = "NMDS",
        choices = list(
          "NMDS",
          "ward.D2",
          "ward.D",
          "single",
          "complete",
          "average",
          "mcquitty",
          "median",
          "centroid"
        )
      )
    )
  })
  
  output$Heatmap <- renderPlot({
    validate(need(data16S(),
                  "Merci d'importer un fichier d'abondance"))
    p <- plot_heatmap(
      physeq = prune_taxa(names(sort(
        taxa_sums(data16S()), decreasing = TRUE
      )[1:input$heatmapTopOtu]), data16S()),
      distance = input$heatmapDist,
      method = input$heatmapMethod,
      sample.order = checkNull(input$heatmapX),
      low = "yellow",
      high = "red",
      na.value = "white"
    )
    if (!is.null(checkNull(input$heatmapGrid))) {
      p <-
        p + facet_grid(paste(".", "~", input$heatmapGrid), scales = "free_x")
    }
    return(p)
  })
  
  output$treeUI <- renderUI({
    validate(need(data16S(), ""))
    box(
      title = "Paramètres",
      width = NULL,
      status = "primary",
      radioButtons(
        "treeRank",
        label = "Niveau taxonomique légendé :",
        choices = c(aucun = "",
                    rank_names(data16S()),
                    OTU = "taxa_names"),
        inline = TRUE
      ),
      sliderInput(
        "treeTopOtu",
        label = "Selection des n OTU les plus abondant :",
        min = 1,
        max = ntaxa(data16S()),
        value = 20
      ),
      checkboxInput("treeRadial", label = "Arbre radial", value = FALSE),
      checkboxInput("treeSample", label = "Show samples", value = TRUE),
      selectInput(
        "treeCol",
        label = "Couleur :",
        choices = c("..." = 0, sample_variables(data16S()))
      ),
      selectInput(
        "treeShape",
        label = "Forme :",
        choices = c("..." = 0, sample_variables(data16S()))
      )
    )
  })
  
  output$tree <- renderPlot({
    validate(
      need(data16S(), "Merci d'importer un fichier d'abondance"),
      need(
        phy_tree(data16S(), errorIfNULL = FALSE),
        "Arbre phylo invalide"
      )
    )
    p <- plot_tree(
      physeq = prune_taxa(names(sort(
        taxa_sums(data16S()), decreasing = TRUE
      )[1:input$treeTopOtu]), data16S()),
      method = ifelse(input$treeSample, "sampledodge", "treeonly"),
      color = checkNull(input$treeCol),
      shape = checkNull(input$treeShape),
      size = "abundance",
      label.tips = checkNull(input$treeRank),
      sizebase = 5,
      ladderize = "left",
      plot.margin = 0
    )
    if (checkNull(input$treeRadial)) {
      return(p + coord_polar(theta = "y"))
    } else {
      return(p)
    }
  })
  
  output$acpUI <- renderUI({
    validate(need(data16S(), ""))
    box(
      title = "Paramètres",
      width = NULL,
      status = "primary",
      selectInput(
        "acpDist",
        label = "Distance :",
        selected = "bray",
        choices = list(
          "bray",
          "jaccard",
          "unifrac",
          "wunifrac",
          "dpcoa",
          "jsd",
          "euclidean"
        )
      ),
      selectInput(
        "acpMethod",
        label = "Methode :",
        selected = "MDS",
        choices = list("DCA", "CCA", "RDA", "CAP", "DPCoA", "NMDS", "MDS", "PCoA")
      ),
      selectInput(
        "acpCol",
        label = "Couleur :",
        choices = c("..." = 0, sample_variables(data16S()))
      ),
      selectInput(
        "acpShape",
        label = "Forme :",
        choices = c("..." = 0, sample_variables(data16S()))
      ),
      selectInput(
        "acpEllipse",
        label = "Ellipses :",
        choices = c("..." = 0, sample_variables(data16S()))
      ),
      selectInput(
        "acpRep",
        label = "Barycentre :",
        choices = c("..." = 0, sample_variables(data16S()))
      )
    )
  })
  
  output$acp <- renderPlot({
    validate(need(data16S(),
                  "Merci d'importer un fichier d'abondance"))
    p <- plot_samples(
      data16S(),
      ordination = ordinate(
        data16S(),
        method = input$acpMethod,
        distance = input$acpDist
      ),
      axes = c(1, 2),
      color = checkNull(input$acpCol),
      replicate = checkNull(input$acpRep),
      shape = checkNull(input$acpShape)
    )
    if (!is.null(checkNull(input$acpEllipse))) {
      p <- p + stat_ellipse(aes_string(group = input$acpEllipse))
    }
    return(p + theme_bw())
  })
})