---
description: Identify genes that are significantly over or
  under-expressed between conditions in specific cell populations.
subtitle:  Bioconductor Toolkit
title:  Differential gene expression
---

<div>

> **Note**
>
> Code chunks run R commands unless otherwise specified.

</div>

In this tutorial we will cover differential gene expression, which
comprises an extensive range of topics and methods. In single cell,
differential expresison can have multiple functionalities such as
identifying marker genes for cell populations, as well as identifying
differentially regulated genes across conditions (healthy vs control).
We will also cover controlling batch effect in your test.

We can first load the data from the clustering session. Moreover, we can
already decide which clustering resolution to use. First let's define
using the `louvain` clustering to identifying differentially expressed
genes.

``` {r}
suppressPackageStartupMessages({
    library(scater)
    library(scran)
    # library(venn)
    library(patchwork)
    library(ggplot2)
    library(pheatmap)
    library(igraph)
    library(dplyr)
})
```

``` {r}
# download pre-computed data if missing or long compute
fetch_data <- TRUE

# url for source and intermediate data
path_data <- "https://export.uppmax.uu.se/naiss2023-23-3/workshops/workshop-scrnaseq"
path_file <- "data/covid/results/bioc_covid_qc_dr_int_cl.rds"
if (!dir.exists(dirname(path_file))) dir.create(dirname(path_file), recursive = TRUE)
if (fetch_data && !file.exists(path_file)) download.file(url = file.path(path_data, "covid/results/bioc_covid_qc_dr_int_cl.rds"), destfile = path_file)
sce <- readRDS(path_file)
print(reducedDims(sce))
```

## Cell marker genes

Let us first compute a ranking for the highly differential genes in each
cluster. There are many different tests and parameters to be chosen that
can be used to refine your results. When looking for marker genes, we
want genes that are positively expressed in a cell type and possibly not
expressed in others.

``` {r}
# Compute differentiall expression
markers_genes <- scran::findMarkers(
    x = sce,
    groups = as.character(sce$louvain_SNNk15),
    lfc = .5,
    pval.type = "all",
    direction = "up"
)

# List of dataFrames with the results for each cluster
markers_genes

# Visualizing the expression of one
markers_genes[["1"]]
```

We can now select the top 25 overexpressed genes for plotting.

``` {r}
# Colect the top 25 genes for each cluster and put the into a single table
top25 <- lapply(names(markers_genes), function(x) {
    temp <- markers_genes[[x]][1:25, 1:2]
    temp$gene <- rownames(markers_genes[[x]])[1:25]
    temp$cluster <- x
    return(temp)
})
top25 <- as_tibble(do.call(rbind, top25))
top25$p.value[top25$p.value == 0] <- 1e-300
top25
```

``` {r}
#| fig-height: 6
#| fig-width: 7

par(mfrow = c(1, 5), mar = c(4, 6, 3, 1))
for (i in unique(top25$cluster)) {
    barplot(sort(setNames(-log10(top25$p.value), top25$gene)[top25$cluster == i], F),
        horiz = T, las = 1, main = paste0(i, " vs. rest"), border = "white", yaxs = "i", xlab = "-log10FC"
    )
    abline(v = c(0, -log10(0.05)), lty = c(1, 2))
}
```

We can visualize them as a heatmap. Here we are selecting the top 5.

``` {r}
#| fig-height: 6
#| fig-width: 8

as_tibble(top25) %>%
    group_by(cluster) %>%
    top_n(-5, p.value) -> top5

scater::plotHeatmap(sce[, order(sce$louvain_SNNk15)],
    features = unique(top5$gene),
    center = T, zlim = c(-3, 3),
    colour_columns_by = "louvain_SNNk15",
    show_colnames = F, cluster_cols = F,
    fontsize_row = 6,
    color = colorRampPalette(c("purple", "black", "yellow"))(90)
)
```

We can also plot a violin plot for each gene.

``` {r}
#| fig-height: 12
#| fig-width: 13

scater::plotExpression(sce, features = unique(top5$gene), x = "louvain_SNNk15", ncol = 5, colour_by = "louvain_SNNk15", scales = "free")
```

## DGE across conditions

The second way of computing differential expression is to answer which
genes are differentially expressed within a cluster. For example, in our
case we have libraries comming from patients and controls and we would
like to know which genes are influenced the most in a particular cell
type. For this end, we will first subset our data for the desired cell
cluster, then change the cell identities to the variable of comparison
(which now in our case is the **type**, e.g. Covid/Ctrl).

``` {r}
#| fig-height: 5
#| fig-width: 5

# Filter cells from that cluster
cell_selection <- sce[, sce$louvain_SNNk15 == 8]

# Compute differentiall expression
DGE_cell_selection <- findMarkers(
    x = cell_selection,
    groups = cell_selection@colData$type,
    lfc = .25,
    pval.type = "all",
    direction = "any"
)
top5_cell_selection <- lapply(names(DGE_cell_selection), function(x) {
    temp <- DGE_cell_selection[[x]][1:5, 1:2]
    temp$gene <- rownames(DGE_cell_selection[[x]])[1:5]
    temp$cluster <- x
    return(temp)
})
top5_cell_selection <- as_tibble(do.call(rbind, top5_cell_selection))
top5_cell_selection
```

We can now plot the expression across the **type**.

``` {r}
#| fig-height: 4
#| fig-width: 6

scater::plotExpression(cell_selection, features = unique(top5_cell_selection$gene), x = "type", ncol = 5, colour_by = "type")
```

#DGE_ALL6.2:

``` {r}
#| fig-height: 8
#| fig-width: 13

plotlist <- list()
for (i in unique(top5_cell_selection$gene)) {
    plotlist[[i]] <- plotReducedDim(sce, dimred = "UMAP_on_MNN", colour_by = i, by_exprs_values = "logcounts") +
        ggtitle(label = i) + theme(plot.title = element_text(size = 20))
}
wrap_plots(plotlist, ncol = 3)
```

## Gene Set Analysis (GSA)

### Hypergeometric enrichment test

Having a defined list of differentially expressed genes, you can now
look for their combined function using hypergeometric test.

``` {r}
# Load additional packages
library(enrichR)

# Check available databases to perform enrichment (then choose one)
enrichR::listEnrichrDbs()

# Perform enrichment
top_DGE <- DGE_cell_selection$Covid[(DGE_cell_selection$Covid$p.value < 0.01) & (abs(DGE_cell_selection$Covid[, grep("logFC.C", colnames(DGE_cell_selection$Covid))]) > 0.25), ]

enrich_results <- enrichr(
    genes     = rownames(top_DGE),
    databases = "GO_Biological_Process_2017b"
)[[1]]
```

Some databases of interest:\
`GO_Biological_Process_2017b``KEGG_2019_Human``KEGG_2019_Mouse``WikiPathways_2019_Human``WikiPathways_2019_Mouse`\
You visualize your results using a simple barplot, for example:

``` {r}
{
    par(mfrow = c(1, 1), mar = c(3, 25, 2, 1))
    barplot(
        height = -log10(enrich_results$P.value)[10:1],
        names.arg = enrich_results$Term[10:1],
        horiz = TRUE,
        las = 1,
        border = FALSE,
        cex.names = .6
    )
    abline(v = c(-log10(0.05)), lty = 2)
    abline(v = 0, lty = 1)
}
```

## Gene Set Enrichment Analysis (GSEA)

Besides the enrichment using hypergeometric test, we can also perform
gene set enrichment analysis (GSEA), which scores ranked genes list
(usually based on fold changes) and computes permutation test to check
if a particular gene set is more present in the Up-regulated genes,
among the DOWN_regulated genes or not differentially regulated.

``` {r}
#| fig-height: 5
#| fig-width: 5

# Create a gene rank based on the gene expression fold change
gene_rank <- setNames(DGE_cell_selection$Covid[, grep("logFC.C", colnames(DGE_cell_selection$Covid))], casefold(rownames(DGE_cell_selection$Covid), upper = T))
```

Once our list of genes are sorted, we can proceed with the enrichment
itself. We can use the package to get gene set from the Molecular
Signature Database (MSigDB) and select KEGG pathways as an example.

``` {r}
#| fig-height: 5
#| fig-width: 5

library(msigdbr)

# Download gene sets
msigdbgmt <- msigdbr::msigdbr("Homo sapiens")
msigdbgmt <- as.data.frame(msigdbgmt)

# List available gene sets
unique(msigdbgmt$gs_subcat)

# Subset which gene set you want to use.
msigdbgmt_subset <- msigdbgmt[msigdbgmt$gs_subcat == "CP:WIKIPATHWAYS", ]
gmt <- lapply(unique(msigdbgmt_subset$gs_name), function(x) {
    msigdbgmt_subset[msigdbgmt_subset$gs_name == x, "gene_symbol"]
})
names(gmt) <- unique(paste0(msigdbgmt_subset$gs_name, "_", msigdbgmt_subset$gs_exact_source))
```

Next, we will run GSEA. This will result in a table containing
information for several pathways. We can then sort and filter those
pathways to visualize only the top ones. You can select/filter them by
either `p-value` or normalized enrichment score (`NES`).

``` {r}
#| fig-height: 5
#| fig-width: 12

library(fgsea)

# Perform enrichemnt analysis
fgseaRes <- fgsea(pathways = gmt, stats = gene_rank, minSize = 15, maxSize = 500, nperm = 10000)
fgseaRes <- fgseaRes[order(fgseaRes$NES, decreasing = T), ]

# Filter the results table to show only the top 10 UP or DOWN regulated processes (optional)
top10_UP <- fgseaRes$pathway[1:10]

# Nice summary table (shown as a plot)
plotGseaTable(gmt[top10_UP], gene_rank, fgseaRes, gseaParam = 0.5)
```

<div>

> **Discuss**
>
> Which KEGG pathways are upregulated in this cluster? Which KEGG
> pathways are dowregulated in this cluster? Change the pathway source
> to another gene set (e.g. CP:WIKIPATHWAYS or CP:REACTOME or
> CP:BIOCARTA or GO:BP) and check the if you get similar results?

</div>

Finally, let's save the integrated data for further analysis.

``` {r}
saveRDS(sce, "data/covid/results/bioc_covid_qc_dr_int_cl_dge.rds")
```

## Session info

```{=html}
<details>
```
```{=html}
<summary>
```
Click here
```{=html}
</summary>
```
``` {r}
sessionInfo()
```

```{=html}
</details>
```