diff --git a/panels/deseq-server.R b/panels/deseq-server.R
index 979b71729a2666ad1f2b73210d5c6956b1251c29..7522154081ff9d7698888d05cb1e471e8ae5ce36 100644
--- a/panels/deseq-server.R
+++ b/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
})
})
diff --git a/panels/deseq-ui.R b/panels/deseq-ui.R
index 941fe163456b51ad9c0eaaa7280c51ff0d947a68..5b3a6b71ea501745b86a8d530586ff476e057937 100644
--- a/panels/deseq-ui.R
+++ b/panels/deseq-ui.R
@@ -1,2 +1,3 @@
deseq <- fluidPage(outputCodeButton(withLoader(plotOutput("deseq", height = 700))),
- uiOutput("deseqUI"))
+ #table
+ uiOutput("deseqUI"))