volcano-plot

panels/deseq-server.R

@@ -8,7 +8,10 @@ output$deseqContrastVarUI <- renderUI({
 })
 
 output$deseqContrastModUI <- renderUI({
-  validate(need(physeq(), ""), need(input$deseqContrastVar, ""))
+  validate(need(physeq(), ""), 
+           need(input$deseqContrastVar, ""),
+           need(class(get_variable(physeq(), input$deseqContrastVar)) != "numeric", "")
+           )
   checkboxGroupInput(
     "deseqContrastMod",
     label = "Contrast (exactly two required) : ",
@@ -38,10 +41,11 @@ output$deseqTitleUI <- renderUI({
 output$deseqPadjUI <- renderUI({
   validate(need(physeq(), ""))
   sliderInput("deseqPadj",
-              label = "Adjusted p-value threshold :",
+              label = "Adjusted p-value threshold (recommended 0.05 ):",
               min = 0,
               max = 1,
-              value = 0.05)
+              value = 0.05,
+              step = 0.01)
 })
 
 output$deseqUI <- renderUI({
@@ -60,59 +64,98 @@ output$deseqUI <- renderUI({
 output$deseq <- metaRender2(renderPlot, {
   validate(
     need(physeq(), "Requires an abundance dataset"),
-    need(length(input$deseqContrastMod) == 2, "Requires two conditions"))
+    need(class(get_variable(physeq(), input$deseqContrastVar)) != "numeric" || 
+           length(input$deseqContrastMod) == 2, "Requires a continuous design or a selection of two modalities for a discrete design.")
+    )
   data <- physeq()
+
+  design <- metaExpr({as.formula(..(paste("~", input$deseqContrastVar)))})
+  cds <- metaExpr({phyloseq_to_deseq2(data, design = design)})
+  dds <- metaExpr({DESeq2::DESeq(cds, sfType = "poscounts")})
   
-  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"
-        )
-      )
+  results <- if (class(get_variable(data, input$deseqContrastVar)) == "numeric") {
+    # First case: regression against a continuous variable
+    metaExpr({
+      DESeq2::results(object = dds,
+                      name = ..(input$deseqContrastVar),
+                      tidy = TRUE) %>%
+        as_tibble() %>%
+        rename(OTU = row) %>%
+        inner_join(tax_table(data) %>% as("matrix") %>% as_tibble(rownames = "OTU"), by = "OTU")
+      })
+    } else {
+      if (length(levels(get_variable(data, input$deseqContrastVar))) == 2) {
+        # Second case: regression against a binary variable
+        metaExpr({
+          DESeq2::results(object = dds,
+                          name = ..(DESeq2::resultsNames(dds)[-1]),
+                          tidy = TRUE) %>%
+            as_tibble() %>% rename(OTU = row) %>%
+            inner_join(tax_table(data) %>% as("matrix") %>% as_tibble(rownames = "OTU"), by = "OTU")
+          })
+      } else {
+        # Third case: regression against a qualiative variable with three or more levels
+        metaExpr({
+          DESeq2::results(object = dds,
+                          contrast = ..(c(input$deseqContrastVar, input$deseqContrastMod[1], input$deseqContrastMod[2])),
+                          tidy = TRUE) %>%
+            as_tibble() %>% rename(OTU = row) %>%
+            inner_join(tax_table(data) %>% as("matrix") %>% as_tibble(rownames = "OTU"), by = "OTU")
+        })
+      }
+    }
+  
+  detail <- if (class(get_variable(data, input$deseqContrastVar)) == "numeric") {
+    # First case
+    metaExpr({
+      ..(paste0("You compare low and high values of the continuous variable ", input$deseqContrastVar, ".\nA positive log2FoldChange means more abundant for high values of ", input$deseqContrastVar, "."))
+    })
+  } else {
+    if (length(levels(get_variable(data, input$deseqContrastVar))) == 2) {
+      # Second case
+      metaExpr({
+        ..(paste0("You compare ", input$deseqContrastMod[1], " to ", input$deseqContrastMod[2], " for the variable ", input$deseqContrastVar, ".\nA positive log2FoldChange means more abundant in ", input$deseqContrastMod[2], " than in ", input$deseqContrastMod[1], "."))
+      })
+    } else {
+      # Third case
+      metaExpr({
+        ..(paste0("You choose to compare ", input$deseqContrastMod[1], " to ", input$deseqContrastMod[2], " for the variable", input$deseqContrastVar, ".\nA positive log2FoldChange means more abundant in ", input$deseqContrastMod[1], " than in ", input$deseqContrastMod[2], "."))
+      })
+    }
+  }
+  
+  deseqTable <- metaExpr({
     
-    # 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() +
+  })
+  
+  deseqPlot <- metaExpr({
+    ggplot(results %>% mutate(evidence = -log10(padj), 
+                              evolution = case_when(
+                                padj <= ..(input$deseqPadj) & log2FoldChange < 0 ~ "Down", 
+                                padj <= ..(input$deseqPadj) & log2FoldChange > 0 ~ "Up", 
+                                TRUE                              ~ "Not DA"
+                              )),
+           aes(x = log2FoldChange, y = evidence)) + 
+      geom_point(aes(color = evolution), size = 1.75, alpha = 0.8, na.rm = T) +     # base layer 
       theme_bw(base_size = 16) +                                                    # clean up theme
-      theme(legend.position = "none") +                                             # remove legend 
-      ggtitle(label = ..(input$deseqTitle)) +                                       # add title
+      theme(legend.position = "none",                                               # remove legend 
+            plot.subtitle = element_text(size = 12)) +                              # add subtitle
+      ggtitle(label = ..(input$deseqTitle), subtitle = detail) +                    # add informative 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")
+      geom_vline(xintercept = 0, colour = "grey80", linetype = 2) +                 # add line at 0
+      geom_hline(yintercept = -log10(..(input$deseqPadj)), colour = "grey80", linetype = 2) +
+      scale_color_manual(values = c("Down" = "red", "Not DA" = "grey20", "Up" = "green")) # change colors
+  })
+  
+  metaExpr({
+    design <- ..(design)
+    cds <- ..(cds)
+    dds <- ..(dds)
+    results <- ..(results)
+    detail <- ..(detail)
+    p <- ..(deseqPlot)
+    p
   })
 })
 

panels/deseq-ui.R

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