Differential gene expression

In this tutorial we will cover about Differetial gene expression, which comprises an extensive range of topics and methods. In single cell, differential expresison can have multiple functionalities such as of identifying marker genes for cell populations, as well as differentially regulated genes across conditions (healthy vs control). We will also exercise on how to account the batch information 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.

suppressPackageStartupMessages({
    library(scater)
    library(scran)
    # library(venn)
    library(cowplot)
    library(ggplot2)
    # library(rafalib)
    library(pheatmap)
    library(igraph)
    library(dplyr)
})

sce <- readRDS("data/results/covid_qc_dr_int_cl.rds")

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 positivelly expressed in a cell type and possibly not expressed in the others.

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

# List of dataFrames with the results for each cluster
markers_genes

## List of length 8
## names(8): 1 2 3 4 5 6 7 8

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

## DataFrame with 18121 rows and 10 columns
##                 p.value          FDR summary.logFC     logFC.2     logFC.3
##               <numeric>    <numeric>     <numeric>   <numeric>   <numeric>
## CD79A      3.12219e-307 5.65771e-303       3.33918     3.50769     3.49569
## MS4A1      2.70271e-242 2.44879e-238       2.34911     2.39361     2.37237
## IGHM       1.34591e-215 8.12973e-212       3.66201     3.82441     3.84624
## TNFRSF13C  6.93770e-201 3.14295e-197       2.05748     2.10440     2.03702
## LINC00926  3.06863e-161 1.11213e-157       1.87617     1.91202     1.91947
## ...                 ...          ...           ...         ...         ...
## AC011043.1            1            1    0.00431686  0.00431686  0.00431686
## AC007325.4            1            1   -0.00900931 -0.00900931 -0.00362024
## AL354822.1            1            1    0.01468456  0.01468456  0.01091020
## AC233755.1            1            1    0.00407473  0.00407473  0.00407473
## AC240274.1            1            1    0.00771093  0.01016628  0.00771093
##                 logFC.4     logFC.5     logFC.6     logFC.7     logFC.8
##               <numeric>   <numeric>   <numeric>   <numeric>   <numeric>
## CD79A           3.33918     3.42574     3.48286     3.45316     3.48489
## MS4A1           2.37006     2.38049     2.34911     2.36439     2.38240
## IGHM            3.66201     3.81737     3.84129     3.81403     3.87463
## TNFRSF13C       2.05748     2.08148     2.06205     2.08998     2.08595
## LINC00926       1.87617     1.91044     1.89659     1.89392     1.91356
## ...                 ...         ...         ...         ...         ...
## AC011043.1  0.002419979  0.00431686  0.00197413  0.00320700  0.00431686
## AC007325.4  0.000964593 -0.01778832 -0.00330748 -0.00525116 -0.00195822
## AL354822.1  0.015844207  0.01290013  0.01070242  0.01376169  0.00073809
## AC233755.1  0.004074732  0.00407473  0.00407473  0.00407473  0.00407473
## AC240274.1 -0.001286558  0.00684239  0.00875890  0.00663751  0.01270911

We can now select the top 25 up regulated genes for plotting.

# 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

We can now select the top 25 up regulated genes for plotting.

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))
}