A Shiny app for visualizing DESeq2 results

Implement from scratch

First we run DESeq2 analysis on the airway dataset:

library(airway)
data(airway)
se = airway
library(DESeq2)
dds = DESeqDataSet(se, design = ~ dex)
keep = rowSums(counts(dds)) >= 10
dds = dds[keep, ]

dds$dex = relevel(dds$dex, ref = "untrt")

dds = DESeq(dds)
res = results(dds)
res = as.data.frame(res)

Since we use InteractiveComplexHeatmap package, we start with the heatmap. Following are the thoughts of this Shiny app:

1. The app starts from the heatmap which shows the significantly differentially expressed genes.
2. When selecting from the heatmap, the selected genes are highlighted in a MA-plot and a volcano plot so it is easy to correspond to these genes’ base means, log2 fold changes and FDRs.
3. There should also be a table which contains the statistics from DESeq2 analysis for the selected genes.
4. The heatmap can be regenerated by different cutoffs for selecting the signficant genes.

We first define several functions:

1. make_heatmap(): make the heatmap for the differential genes under specific cutoffs.
2. make_maplot(): make the MA-plot and highlighted a subset of selected genes.
3. make_volcano(): make the volcano plot and highlighted a subset of selected genes.

make_maplot() and make_volcano() are basically scatter plot and they are defined very similarly.

The heatmap only contains the information of the differential genes, which means, all the indices captured with InteractiveComplexHeatmap only restricted to the matrix of the differential genes. In make_maplot() and make_volcano(), they expect the indices of the genes to be from the complete gene sets, thus, we create an environment variable env which saves the indices of the differential genes in the complete and share between functions.

library(InteractiveComplexHeatmap)
library(ComplexHeatmap)
library(circlize)
library(GetoptLong)

env = new.env()

make_heatmap = function(fdr = 0.01, base_mean = 0, log2fc = 0) {
    l = res$padj <= fdr & res$baseMean >= base_mean & 
                        abs(res$log2FoldChange) >= log2fc; l[is.na(l)] = FALSE

    if(sum(l) == 0) return(NULL)

    m = counts(dds, normalized = TRUE)
    m = m[l, ]

    env$row_index = which(l)

    ht = Heatmap(t(scale(t(m))), name = "z-score",
        top_annotation = HeatmapAnnotation(
            dex = colData(dds)$dex,
            sizeFactor = anno_points(colData(dds)$sizeFactor)
        ),
        show_row_names = FALSE, show_column_names = FALSE, row_km = 2,
        column_title = paste0(sum(l), " significant genes with FDR < ", fdr),
        show_row_dend = FALSE) + 
        Heatmap(log10(res$baseMean[l]+1), show_row_names = FALSE, width = unit(5, "mm"),
            name = "log10(baseMean+1)", show_column_names = FALSE) +
        Heatmap(res$log2FoldChange[l], show_row_names = FALSE, width = unit(5, "mm"),
            name = "log2FoldChange", show_column_names = FALSE,
            col = colorRamp2(c(-2, 0, 2), c("green", "white", "red")))
    ht = draw(ht, merge_legend = TRUE)
    ht
}

# make the MA-plot with some genes highlighted
make_maplot = function(res, highlight = NULL) {
    col = rep("#00000020", nrow(res))
    cex = rep(0.5, nrow(res))
    names(col) = rownames(res)
    names(cex) = rownames(res)
    if(!is.null(highlight)) {
        col[highlight] = "red"
        cex[highlight] = 1
    }
    x = res$baseMean
    y = res$log2FoldChange
    y[y > 2] = 2
    y[y < -2] = -2
    col[col == "red" & y < 0] = "darkgreen"
    par(mar = c(4, 4, 1, 1))

    suppressWarnings(
        plot(x, y, col = col, 
            pch = ifelse(res$log2FoldChange > 2 | res$log2FoldChange < -2, 1, 16), 
            cex = cex, log = "x",
            xlab = "baseMean", ylab = "log2 fold change")
    )
}

# make the volcano plot with some genes highlited
make_volcano = function(res, highlight = NULL) {
    col = rep("#00000020", nrow(res))
    cex = rep(0.5, nrow(res))
    names(col) = rownames(res)
    names(cex) = rownames(res)
    if(!is.null(highlight)) {
        col[highlight] = "red"
        cex[highlight] = 1
    }
    x = res$log2FoldChange
    y = -log10(res$padj)
    col[col == "red" & x < 0] = "darkgreen"
    par(mar = c(4, 4, 1, 1))

    suppressWarnings(
        plot(x, y, col = col, 
            pch = 16, 
            cex = cex,
            xlab = "log2 fold change", ylab = "-log10(FDR)")
    )
}

We use shinydashboard to organize the interactive heatmap components. It has a three-column layout:

• The first column: the original heatmap,
• The second column: the sub-heatmap and the default output,
• The thrid column: the self-defined output.

We separatedly specifiy the three interaction heatmap components by originalHeatmapOutput(), subHeatmapOutput() and HeatmapInfoOutput().

One additional htmlOutput("note") is only for printing some friendly message in the app.

library(shiny)
library(shinydashboard)
body = dashboardBody(
    fluidRow(
        column(width = 4,
            box(title = "Differential heatmap", width = NULL, solidHeader = TRUE, status = "primary",
                originalHeatmapOutput("ht", height = 800, containment = TRUE)
            )
        ),
        column(width = 4,
            id = "column2",
            box(title = "Sub-heatmap", width = NULL, solidHeader = TRUE, status = "primary",
                subHeatmapOutput("ht", title = NULL, containment = TRUE)
            ),
            box(title = "Output", width = NULL, solidHeader = TRUE, status = "primary",
                HeatmapInfoOutput("ht", title = NULL)
            ),
            box(title = "Note", width = NULL, solidHeader = TRUE, status = "primary",
                htmlOutput("note")
            ),
        ),
        column(width = 4,
            box(title = "MA-plot", width = NULL, solidHeader = TRUE, status = "primary",
                plotOutput("ma_plot")
            ),
            box(title = "Volcanno plot", width = NULL, solidHeader = TRUE, status = "primary",
                plotOutput("volcanno_plot")
            ),
            box(title = "Result table of the selected genes", width = NULL, solidHeader = TRUE, status = "primary",
                DTOutput("res_table")
            )
        ),
        tags$style("
            .content-wrapper, .right-side {
                overflow-x: auto;
            }
            .content {
                min-width:1500px;
            }
        ")
    )
)

In previous code, we add self-defined CSS code to set the minimal width of the area that contains the boxes.

We will have several self-defined actions to respond brushing event on heatmap. We put all these actions into one brush_action() call. Note here we use env$row_index[row_index] to get the indices of the selected genes that correspond to the complete gene set.

library(DT)
library(GetoptLong) # for qq() function
brush_action = function(df, input, output, session) {
    
    row_index = unique(unlist(df$row_index))
    selected = env$row_index[row_index]
        
    output[["ma_plot"]] = renderPlot({
        make_maplot(res, selected)
    })

    output[["volcanno_plot"]] = renderPlot({
        make_volcano(res, selected)
    })

    output[["res_table"]] = renderDT(
        formatRound(datatable(res[selected, c("baseMean", "log2FoldChange", "padj")], rownames = TRUE), columns = 1:3, digits = 3)
    )

    output[["note"]] = renderUI({
        if(!is.null(df)) {
            HTML(qq("<p>Row indices captured in <b>Output</b> only correspond to the matrix of the differential genes. To get the row indices in the original matrix, you need to perform:</p>
<pre>
l = res$padj <= @{input$fdr} & 
    res$baseMean >= @{input$base_mean} & 
    abs(res$log2FoldChange) >= @{input$log2fc}
l[is.na(l)] = FALSE
which(l)[row_index]
</pre>
<p>where <code>res</code> is the complete data frame from DESeq2 analysis and <code>row_index</code> is the <code>row_index</code> column captured from the code in <b>Output</b>.</p>"))
        }
    })
}

Side bar contains settings for cutoffs to select significant genes, i.e., FDR, base mean and log2 fold change. We also add an action button to trigger the generation of heatmap.

ui = dashboardPage(
    dashboardHeader(title = "DESeq2 results"),
    dashboardSidebar(
        selectInput("fdr", label = "Cutoff for FDRs:", c("0.001" = 0.001, "0.01" = 0.01, "0.05" = 0.05)),
        numericInput("base_mean", label = "Minimal base mean:", value = 0),
        numericInput("log2fc", label = "Minimal abs(log2 fold change):", value = 0),
        actionButton("filter", label = "Generate heatmap")
    ),
    body
)

On the server side, we only respond to the action button, because heatmap generation does not happen instantly (1 to 2 seconds lag). We want to generate the heatmap only when all the cutoffs are configured by users.

ignoreNULL = FALSE is set to ensure that the heatmap will be generated right after the app is loaded.

server = function(input, output, session) {
    observeEvent(input$filter, {
        ht = make_heatmap(fdr = as.numeric(input$fdr), base_mean = input$base_mean, log2fc = input$log2fc)
        if(!is.null(ht)) {
            makeInteractiveComplexHeatmap(input, output, session, ht, "ht",
                brush_action = brush_action)
        } else {
            # The ID for the heatmap plot is encoded as @{heatmap_id}_heatmap, thus, it is ht_heatmap here.
            output$ht_heatmap = renderPlot({
                grid.newpage()
                grid.text("No row exists after filtering.")
            })
        }
    }, ignoreNULL = FALSE)
}

Everything is ready, and we generate the app with shinyApp():

shinyApp(ui, server)

This app can be directly generated with htShinyExample(10.5). A demonstration of the app is in following figure:

The function interactivate()

InteractiveComplexHeatmap has a generic function interactivate() which aims to provide an API to generate Shiny apps for objects that contain results for specific analysis. Currently, it only has an implementation for the DESeqDataSet object, which is from DESeq2 analysis. It can be simply used as:

library(DESeq2)
library(airway)
data(airway)

dds = DESeqDataSet(airway, design = ~ dex)
dds = DESeq(dds)

interactivate(dds)