---
title: "A Shiny app for visualizing DESeq2 results"
author: "Zuguang Gu ( z.gu@dkfz.de )"
date: "`r Sys.Date()`"
output: 
    rmarkdown::html_vignette:
        width: 8
        fig_width: 5
vignette: >
  %\VignetteIndexEntry{6. A Shiny app for visualizing DESeq2 results}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

<style>
p {
    margin: 1em 0;
}
img {
    background-color: #FFFFFF;
    padding: 2px;
    border: 1px solid #DDDDDD;
    border-radius: 3px;
    border: 1px solid #CCCCCC;
    margin: 0 5px;
}
</style>

In this vignette, I demonstrate how to implement a Shiny app for visualizing DESeq2 results with **InteractiveComplexHeatmap** package.

## Implement from scratch

First we run DESeq2 analysis on the **airway** dataset:

```{r, eval = FALSE}
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.


```{r, eval = FALSE}
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.

```{r, eval = FALSE}
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.

```{r, eval = FALSE}
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.

```{r, eval = FALSE}
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.

```{r, eval = FALSE}
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()`:

```{r, eval = FALSE}
shinyApp(ui, server)
```

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

<script>
document.write('<img src="https://user-images.githubusercontent.com/449218/111832925-b16ae280-88f1-11eb-8530-290374f9f2c2.gif" width="100%" />');
</script>


## 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:

```{r, eval = FALSE}
library(DESeq2)
library(airway)
data(airway)

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

interactivate(dds)
```

<br>
<br>
<br>
<br>
<br>
<br>