road to volcano

panels/deseq-server.R

+output$deseqContrastVarUI <- renderUI({
+  validate(need(physeq(), ""))
+  selectInput(
+    "deseqContrastVar",
+    label = "Experimental design : ",
+    choices = c(sample_variables(physeq()))
+  )
+})
+
+output$deseqContrastModUI <- renderUI({
+  validate(need(physeq(), ""), need(input$deseqContrastVar, ""))
+  checkboxGroupInput(
+    "deseqContrastMod",
+    label = "Contrast (exactly two required) : ",
+    choices = NULL,
+    inline = TRUE
+  )
+})
+
+observe({
+  validate(need(physeq(), ""), need(input$deseqContrastVar, ""))
+  var <- levels(as.factor(get_variable(physeq(), input$deseqContrastVar)))
+  updateCheckboxGroupInput(session,
+                           inputId = "deseqContrastMod",
+                           choices = var,
+                           selected = var[c(1, 2)],
+                           inline = TRUE
+  )
+})
+
+output$deseqTitleUI <- renderUI({
+  validate(need(physeq(), ""))
+  textInput("deseqTitle",
+            label = "Title : ",
+            value = "Volcano Plot")
+})
+
+output$deseqPadjUI <- renderUI({
+  validate(need(physeq(), ""))
+  sliderInput("deseqPadj",
+              label = "Adjusted p-value threshold :",
+              min = 0,
+              max = 1,
+              value = 0.05)
+})
+
+output$deseqUI <- renderUI({
+  validate(need(physeq(), ""))
+  box(
+    title = "Setting : " ,
+    width = NULL,
+    status = "primary",
+    uiOutput("deseqContrastVarUI"),
+    uiOutput("deseqContrastModUI"),
+    uiOutput("deseqTitleUI"),
+    uiOutput("deseqPadjUI")
+  )
+})
+
+output$deseq <- metaRender2(renderPlot, {
+  validate(
+    need(physeq(), "Requires an abundance dataset"),
+    need(length(input$deseqContrastMod) == 2, "Requires two conditions"))
+  data <- physeq()
+  
+  metaExpr({
+    deseq_data <- phyloseq_to_deseq2(data, ..(as.formula(paste("~", input$deseqContrastVar))))
+    dds <- DESeq2::DESeq(deseq_data, sfType = "poscounts")
+    
+    results <- DESeq2::results(dds, 
+                               tidy = TRUE, 
+                               contrast = ..(c(input$deseqContrastVar, input$deseqContrastMod[1], input$deseqContrastMod[2]))) %>% 
+      rename(OTU = row)
+    
+    da_volcano <- data.frame(
+      #otu      = row.names(results),
+      otu      = results$OTU,
+      evidence = -log10(results$padj), 
+      lfc      = results$log2FoldChange) %>% 
+      na.omit()
+    # add a threshol line
+    y_axix_volcano_line <- -log10(..(input$deseqPadj))
+    # Modify dataset to add new coloumn of colors
+    da_volcano <- da_volcano %>%
+      mutate(
+        color = case_when(
+          lfc > 0 & evidence > y_axix_volcano_line ~ "More", 
+          lfc < 0 & evidence > y_axix_volcano_line ~ "Less", 
+          TRUE                                     ~ "Equal"
+        )
+      )
+    
+    # Color corresponds to fold change directionality
+    volcano_plot <- ggplot(da_volcano, 
+                           aes(x = lfc, y = evidence, label = otu)) + 
+      geom_point(aes(color = factor(color)), size = 1.75, alpha = 0.8, na.rm = T) + # add gene points
+      geom_text() +
+      theme_bw(base_size = 16) +                                                    # clean up theme
+      theme(legend.position = "none") +                                             # remove legend 
+      ggtitle(label = ..(input$deseqTitle)) +                                       # add title
+      xlab(expression(log[2]("FoldChange"))) +                                      # x-axis label
+      ylab(expression(-log[10]("adjusted p-value"))) +                              # y-axis label
+      geom_vline(xintercept = 0, colour = "grey80", linetype = 2) +                                # add line at 0
+      geom_hline(aes(yintercept = y_axix_volcano_line), yintercept = y_axix_volcano_line, colour = "grey80", linetype = 2) +
+      annotate(geom = "text", 
+               label = paste("padj =", ..(input$deseqPadj)), 
+               x = min(da_volcano$lfc), 
+               y = y_axix_volcano_line + 0.25, 
+               size = 4,
+               colour = "black",
+               vjust = 0,
+               hjust = 0) + # add pvalue threshold           
+      scale_color_manual(values = c("More" = "red", "Less" = "chartreuse", "Equal" = "darkgray")) # change colors
+    # Plot figure
+    volcano_plot + scale_y_continuous(trans = "log1p")
+  })
+})
+
+observeEvent(input$deseq_output_code,
+             {
+               displayCodeModal(
+                 expandChain(
+                   quote(library(phyloseq)),
+                   quote(library(phyloseq.extended)),
+                   quote(library(DESeq2)),
+                   quote(library(ggplot2)),
+                   quote(library(magrittr)),
+                   quote(library(dplyr)),
+                   
+                   "# Replace `data` with you own data.",
+                   output$deseq()
+                 )
+               )
+             }
+)

panels/deseq-ui.R

+deseq <- fluidPage(outputCodeButton(withLoader(plotOutput("deseq", height = 700))),
+                     uiOutput("deseqUI"))

server.R

@@ -23,6 +23,7 @@ shinyServer
   source("panels/richnessA-server.R", local = TRUE)
   source("panels/richnessB-server.R", local = TRUE)
   source("panels/pca-server.R", local = TRUE)
+  source("panels/deseq-server.R", local = TRUE)
   source("panels/tree-server.R", local = TRUE)
 
   physeq <- reactiveVal()

ui.R

@@ -10,6 +10,7 @@ source("panels/rarefactionCurve-ui.R", local = TRUE)
 source("panels/richnessA-ui.R", local = TRUE)
 source("panels/richnessB-ui.R", local = TRUE)
 source("panels/pca-ui.R", local = TRUE)
+source("panels/deseq-ui.R", local = TRUE)
 source("panels/tree-ui.R", local = TRUE)
 source("panels/Help-ui.R", local = TRUE)
 
@@ -75,6 +76,7 @@ dashboardHeader(title = "Easy16S"),
                menuSubItem("Table", tabName = "betaTable")
       ),
       menuItem("PCA", tabName = "pca", icon = icon("bullseye")),
+      menuItem("Differential expression analysis", tabName = "deseq", icon = icon("balance-scale-left")),
       menuItem("Phylogenetic tree", tabName = "tree", icon = icon("tree")),
       menuItem("Help", tabName = "Help", icon = icon("info-circle"))
   )),
@@ -97,6 +99,7 @@ dashboardHeader(title = "Easy16S"),
       tabItem(tabName = "betaHeatmap", betaHeatmap),
       tabItem(tabName = "betaTable", betaTable),
       tabItem(tabName = "pca", pca),
+      tabItem(tabName = "deseq", deseq),
       tabItem(tabName = "tree", tree),
       tabItem(tabName = "Help", Help)
       )