Differential expression testing

Load in the data

This vignette highlights some example workflows for performing differential expression in Seurat. For demonstration purposes, we will be using the 2,700 PBMC object that is created in the first guided tutorial. You can download the pre-computed object here.

library(Seurat)
pbmc <- readRDS(file = "~/Projects/datasets/pbmc3k_final.rds")

Perform default differential expression tests

The bulk of Seurat’s differential expression features can be accessed through the FindMarkers function. As a default, Seurat performs differential expression based on the non-parameteric Wilcoxon rank sum test. This replaces the previous default test (‘bimod’). To test for differential expression between two specific groups of cells, specify the ident.1 and ident.2 parameters.

# list options for groups to perform differential expression on
levels(pbmc@ident)

## [1] "B cells"           "CD4 T cells"       "CD8 T cells"      
## [4] "CD14+ Monocytes"   "Dendritic cells"   "FCGR3A+ Monocytes"
## [7] "Megakaryocytes"    "NK cells"

# Find differentially expressed genes between CD14+ and FCGR3A+ Monocytes
monocyte_de_genes <- FindMarkers(pbmc, ident.1 = "CD14+ Monocytes", ident.2 = "FCGR3A+ Monocytes")

# view results
head(monocyte_de_genes)

##               p_val avg_logFC pct.1 pct.2    p_val_adj
## FCGR3A 2.655706e-96 -2.569489 0.134 0.962 3.642036e-92
## LYZ    7.356623e-71  1.738480 1.000 0.987 1.008887e-66
## RHOC   5.664625e-65 -1.602845 0.165 0.854 7.768466e-61
## S100A8 5.369175e-64  2.656454 0.973 0.503 7.363287e-60
## IFITM2 2.099722e-62 -1.432368 0.676 1.000 2.879558e-58
## S100A9 1.540888e-61  2.250824 0.996 0.879 2.113174e-57

The results data frame has the following columns :

p_val : p_val (unadjusted)
avg_logFC : log fold-chage of the average expression between the two groups. Positive values indicate that the gene is more highly expressed in the first group.
pct.1 : The percentage of cells where the gene is detected in the first group
pct.2 : The percentage of cells where the gene is detected in the second group
p_val_adj : Adjusted p-value, based on bonferroni correction using all genes in the dataset.

If the ident.2 parameter is omitted or set to NULL, FindMarkers will test for differentially expressed genes between the group specified by ident.1 and all other cells.

# Find differentially expressed genes between CD14+ Monocytes and all other
# cells, only search for positive markers
monocyte_de_genes <- FindMarkers(pbmc, ident.1 = "CD14+ Monocytes", ident.2 = NULL, 
    only.pos = TRUE)

# view results
head(monocyte_de_genes)

##                p_val avg_logFC pct.1 pct.2     p_val_adj
## S100A9  0.000000e+00  3.830954 0.996 0.216  0.000000e+00
## S100A8  0.000000e+00  3.786963 0.973 0.122  0.000000e+00
## LGALS2  0.000000e+00  2.639056 0.908 0.060  0.000000e+00
## FCN1    0.000000e+00  2.367232 0.954 0.151  0.000000e+00
## CD14   5.052830e-289  1.947447 0.662 0.029 6.929451e-285
## TYROBP 1.281069e-283  2.111357 0.994 0.265 1.756858e-279

Prefilter genes or cells to increase the speed of DE testing

To increase the speed of marker discovery, particularly for large datasets, Seurat allows for pre-filtering of genes or cells. For example, genes that are very infrequently detected in either group of cells, or genes that are expressed at similar average levels, are unlikely to be differentially expressed. Example use cases of the min.pct, logfc.threshold, min.diff.pct, and max.cells.per.ident parameters are demonstrated below.

# Pre-filter genes that are detected at <50% frequency in either CD14+
# Monocytes or FCGR3A+ Monocytes
monocyte_de_genes <- FindMarkers(pbmc, ident.1 = "CD14+ Monocytes", ident.2 = "FCGR3A+ Monocytes", 
    min.pct = 0.5)

# Pre-filter genes that have less than a two-fold change between the average
# expression of CD14+ Monocytes vs FCGR3A+ Monocytes
monocyte_de_genes <- FindMarkers(pbmc, ident.1 = "CD14+ Monocytes", ident.2 = "FCGR3A+ Monocytes", 
    logfc.threshold = log(2))

# Pre-filter genes whose detection percentages across the two groups are
# similar (within 0.25)
monocyte_de_genes <- FindMarkers(pbmc, ident.1 = "CD14+ Monocytes", ident.2 = "FCGR3A+ Monocytes", 
    min.diff.pct = 0.25)

# Increasing min.pct, logfc.threshold, and min.diff.pct, will increase the
# speed of DE testing, but could also miss genes that are prefiltered

# Subsample each group to a maximum of 200 cells. Can be very useful for
# large clusters, or computationally-intensive DE tests
monocyte_de_genes <- FindMarkers(pbmc, ident.1 = "CD14+ Monocytes", ident.2 = "FCGR3A+ Monocytes", 
    max.cells.per.ident = 200)

Perform DE analysis using alternative tests

The following differential expression tests are currently supported:

“wilcox” : Wilcoxon rank sum test (default)
“bimod” : Likelihood-ratio test for single cell gene expression, (McDavid et al., Bioinformatics, 2013)
“roc” : Standard AUC classifier
“t” : Student’s t-test
“tobit” : Tobit-test for differential gene expression (Trapnell et al., Nature Biotech, 2014)
“poisson” : Likelihood ratio test assuming an underlying negative binomial distribution. Use only for UMI-based datasets
“negbinom” : Likelihood ratio test assuming an underlying negative binomial distribution. Use only for UMI-based datasets
“MAST” : GLM-framework that treates cellular detection rate as a covariate (Finak et al, Genome Biology, 2015) INSTALL
“DESeq2” : DE based on a model using the negative binomial distribution (Love et al, Genome Biology, 2014) INSTALL

For MAST and DESeq2 please ensure that these packages are installed separately in order to use them as part of Seurat. Once installed, use the test.use parameter can be used to specify which DE test to use.

# Test for DE genes using the MAST package
monocyte_de_genes <- FindMarkers(pbmc, ident.1 = "CD14+ Monocytes", ident.2 = "FCGR3A+ Monocytes", 
    test.use = "MAST")

# Test for DE genes using the DESeq2 package. Throws an error if DESeq2 has
# not already been installed Note that the DESeq2 workflows can be
# computationally intensive for large datasets, but are incompatible with
# some gene pre-filtering options We therefore suggest initially limiting
# the number of cells used for testing
monocyte_de_genes <- FindMarkers(pbmc, ident.1 = "CD14+ Monocytes", ident.2 = "FCGR3A+ Monocytes", 
    test.use = "DESeq2", max.cells.per.ident = 50)

Acknowledgements

We thank the authors of the MAST and DESeq2 packages for their kind assistance and advice. We also point users to the following study by Charlotte Soneson and Mark Robinson, which performs careful and extensive evaluation of methods for single cell differential expression testing.