This notebook implements an efficient version of pseudobulk nb-glm based differential expression analysis with DESeq2. Pseudobulk means that all reads from a single batch group (e.g. donor) get pooled into a single observation.
In general, pseudobulk is statistically preferable to but much slower than Wilcoxon, especially when you need to consider covariates. A more robust but considerably slower alternative to pseudobulk is including donors as random effects. Random effects are preferable for small cell count groups but likely give similar results to pseudobulk estimates for large groups.
This idea is not at all new. The earliest reference I know for is from Lun et al:
1) To find markers most upregulated in a cluster, I divide samples into those in and out of the cluster. An alternative is to let each out group remain an independent pseudobulk sample. This is in fact the recommended way from Mike Love: https://support.bioconductor.org/p/118090/. While this is certainly faster than re-estimate size factors for each cluster-specific analysis, I find it gives strange results. Namely, I get more inflated p-values and significant p-values for the wrong canonical marker genes (e.g. CD14 for B cells).
2) On my laptop, it takes ~20 seconds to run do ~3000 genes from 2700 cells, 3 donors, 2 batches, and 9 cell types.
%% Cell type:markdown id: tags:
# Load some data
%% Cell type:code id: tags:
``` R
suppressPackageStartupMessages({
library(tidyverse)
library(presto)
library(singlecellmethods)
library(SeuratData)
library(Seurat)
library(DESeq2)
})
```
%% Output
Warning message:
“package ‘DESeq2’ was built under R version 3.6.1”Warning message:
“package ‘S4Vectors’ was built under R version 3.6.1”Warning message:
“package ‘BiocGenerics’ was built under R version 3.6.1”Warning message:
“package ‘IRanges’ was built under R version 3.6.1”Warning message:
“package ‘GenomicRanges’ was built under R version 3.6.1”Warning message:
“package ‘GenomeInfoDb’ was built under R version 3.6.1”Warning message:
“package ‘SummarizedExperiment’ was built under R version 3.6.1”Warning message:
“package ‘Biobase’ was built under R version 3.6.1”Warning message:
“package ‘DelayedArray’ was built under R version 3.6.1”Warning message:
“package ‘BiocParallel’ was built under R version 3.6.1”
In this artificial example, donor and batch are fictitious, so DESeq2's GLM $\beta$ estimates should not be that different from the Wilcoxon estimates. Here, we'll compare $\beta$s to auROC, which is essentially equivalent to the Wilxocon statistic.
“Column `group` joining character vector and factor, coercing into character vector”
%% Cell type:markdown id: tags:
Most of the results agree, more or less. Interestingly, the Wilcoxon labels almost all genes as upregulated in DCs and CD16+ Monocytes and downregulated in Platelets. What's going on here? It turns out that DCs and CD16+ Monos are mRNA rich cells while platelets are mRNA poor cells overall. DESeq2 is able to account for this effect better than CP10K normalization.
