3.9 多样本间差异分析
刘小泽写于2020.7.10+7.15
scRNA-seq的一个强大之处在于:方便多个实验条件下的样本比较,比如可以检测药物处理后 (Richard et al. 2018) 或基因修饰后 (Scialdone et al. 2016) 的细胞类型变化,比起常规的单一实验因素能提供更多信息。
对于多样本的scRNA差异分析,主要分成两大类:差异表达(differential expression)和差异丰度(differential abundance),前者检测相同类型细胞在不同条件下表达的变化,后者检测不同条件下细胞类型(或状态等)组成的变化。下面就利用小鼠早期胚胎细胞数据(Pijuan-Sala et al. 2019) 来展示这两种分析。
数据来自小鼠E8.5时期的嵌合体胚胎,每个嵌合体胚胎都是将td-Tomato-positive的胚胎干细胞(ESCs)注射到野生型囊胚中得到的。这里注射的细胞和囊胚细胞除了td-Tomato报告基因,没有其他基因上的差异。利用野生型嵌合体进行研究的目的是为了确定注射过程本身是否引入了其他的差异。
在小鼠胚胎干细胞(ES 细胞)中进行DNA 同源重组,将ES 细胞重新注射到囊胚腔中,可以形成嵌合胚胎,并能在假孕小鼠体内发育,最终发育为嵌合体小鼠,是判定胚胎干细胞系是否具有种系分化能力的唯一方法
sce.chimera
数据总共包含6个样本,包含3个重复批次,每个批次2个处理,一个处理样本属于 td-Tomato阳性细胞,另一个处理样本属于阴性。样本是从6-7个嵌合胚胎中利用荧光活化分选得到的。每个样本利用10X Genomics方法建库(Zheng et al. 2017),得到2000-7000个细胞。其中merging这一步中指定出了批次信息,方便更好地进行数据整合
#--- loading ---#
library(MouseGastrulationData)
sce.chimera <- WTChimeraData(samples=5:10)
sce.chimera
#--- feature-annotation ---#
library(scater)
rownames(sce.chimera) <- uniquifyFeatureNames(
rowData(sce.chimera)$ENSEMBL, rowData(sce.chimera)$SYMBOL)
#--- quality-control ---#
drop <- sce.chimera$celltype.mapped %in% c("stripped", "Doublet")
sce.chimera <- sce.chimera[,!drop]
#--- normalization ---#
sce.chimera <- logNormCounts(sce.chimera)
#--- variance-modelling ---#
library(scran)
dec.chimera <- modelGeneVar(sce.chimera, block=sce.chimera$sample)
chosen.hvgs <- dec.chimera$bio > 0
#--- merging ---#
library(batchelor)
set.seed(01001001)
merged <- correctExperiments(sce.chimera,
batch=sce.chimera$sample,
subset.row=chosen.hvgs,
PARAM=FastMnnParam(
merge.order=list(
list(1,3,5), # WT (3 replicates)
list(2,4,6) # td-Tomato (3 replicates)
)
)
)
#--- clustering ---#
g <- buildSNNGraph(merged, use.dimred="corrected")
clusters <- igraph::cluster_louvain(g)
colLabels(merged) <- factor(clusters$membership)
#--- dimensionality-reduction ---#
merged <- runTSNE(merged, dimred="corrected", external_neighbors=TRUE)
merged <- runUMAP(merged, dimred="corrected", external_neighbors=TRUE)
拿到数据,先要有个初步的认识
> merged
class: SingleCellExperiment
dim: 14699 19426
metadata(2): merge.info pca.info
assays(3): reconstructed counts logcounts
rownames(14699): Xkr4 Rp1 ... Vmn2r122
CAAA01147332.1
rowData names(3): rotation ENSEMBL SYMBOL
colnames(19426): cell_9769 cell_9770 ...
cell_30701 cell_30702
colData names(13): batch cell ... sizeFactor
label
reducedDimNames(3): corrected TSNE UMAP
altExpNames(0):
# 有哪些样本信息
> names(colData(merged))
[1] "batch" "cell"
[3] "barcode" "sample"
[5] "stage" "tomato"
[7] "pool" "stage.mapped"
[9] "celltype.mapped" "closest.cell"
[11] "doub.density" "sizeFactor"
[13] "label"
# 看下三个批次
> table(merged$pool)
3 4 5
3324 5644 10458
# 看下两种处理
> table(merged$tomato)
FALSE TRUE
10331 9095
# 看下分群结果
> table(colLabels(merged))
1 2 3 4 5 6 7 8 9 10
947 112 868 680 606 507 2208 424 1259 252
11 12 13 14 15 16 17 18 19 20
104 727 58 423 1044 872 432 1264 454 1022
21 22 23 24 25 26
211 160 156 1640 860 2136
# 当然也可以把三个批次和分群结果放一起看
table(colLabels(merged), merged$pool)
# 或者把两个处理和分群结果放一起看
table(colLabels(merged), merged$tomato)
# 再直观一点,看看数据整合的效果如何,也就是批次效应处理如何
gridExtra::grid.arrange(
plotTSNE(merged, colour_by="tomato", text_by="label"),
plotTSNE(merged, colour_by=data.frame(pool=factor(merged$pool))),
ncol=2
)

其实现在有了分群结果,下一步应该进行marker基因检测以及细胞注释了,但这里的目的不是这些,所以这些步骤略去了,直接使用(Pijuan-Sala et al. 2019) 小鼠早期胚胎发育数据集做的的细胞类型。这个信息放在了
merged$celltype.mapped
中。不过为了看一下使用他人注释数据用在我们自己的数据质量如何,可以用热图来看二者重叠程度by.label <- table(colLabels(merged), merged$celltype.mapped)
pheatmap::pheatmap(log2(by.label+1), cluster_cols=FALSE, cluster_rows=FALSE,
color=viridis::viridis(101))
通过这个图可以看到,我们这里存在多对多的关系。也反映出:针对发育分化的细胞,细胞类型注释是个难题。

这也是最最常见的分析思路,一般拿到转录组表达量就会先想着把差异分析做一做 针对单细胞的差异分析方法有多种,一般的综述都会提及。这里使用的方法是“模拟bulk转录组表达量”的模式(Tung et al. 2017)
它的思路是:根据不同细胞类型和不同分群,把对应的细胞表达量加起来。比如根据上面得到的细胞分群信息和注释信息
> colData(merged)[,c("celltype.mapped", "sample")]
DataFrame with 19426 rows and 2 columns
celltype.mapped sample
<character> <integer>
cell_9769 Mesenchyme 5
cell_9770 Endothelium 5
cell_9771 Allantois 5
cell_9772 Erythroid3 5
cell_9773 Erythroid1 5
... ... ...
cell_30698 Erythroid2 10
cell_30699 Erythroid3 10
cell_30700 Surface ectoderm 10
cell_30701 Forebrain/Midbrain/Hindbrain 10
cell_30702 Erythroid3 10
然后根据这两个信息作为分组,将各自的组内细胞表达量分别加起来。
summed <- aggregateAcrossCells(merged,
id=colData(merged)[,c("celltype.mapped", "sample")])
> dim(merged);dim(summed)
[1] 14699 19426
[1] 14699 186
可见大大降低了样本数量,原来是一个细胞当做一个样本对待;现在是整合后的一组细胞当做一个样本。这样做的原因有以下几个:
- bulk转录组软件的标准差异分析流程需要比较大的表达量,而单细胞中存在很多0表达量
- 将多个细胞汇成一组,也不是随便挑的几个细胞,就像上面它是根据两个条件选的同组细胞。就像是生物体内的复制过程,虽然这些细胞之间也存在着一些差异,但相比于组外的细胞,它们之间更相似。而如果直接把每个细胞的表达量传给差异分析工具,它会认为每个细胞就是一个独立的生物样本,这也是和真实情况不符的
差异分析基于edgeR的 quasi-likelihood (QL)方法,使用负二项广义线性模型(NB GLM) 来处理过度分散的表达量数据。
先来回忆一下bulk转录组的edegR分析怎么做吧:
rm(list = ls())
options(stringsAsFactors = F)
load('airway.expreSet.Rdata')
library(edgeR)
## 第一步:创建edgeR对象
# 这里需要表达矩阵和分组信息
e <- DGEList(counts=expr,group=factor(grp))
## 第二步:预处理
# 进行一下过滤
keep <- rowSums(cpm(e)>1) >= 2
e <- e[keep, , keep.lib.sizes=FALSE]
# 进行一下校正
e$samples$lib.size <- colSums(e$counts)
e <- calcNormFactors(e)
e$samples
DEG=e
## 第三步:创建模型并分析
# 这里得到分组矩阵(实验设计矩阵)
design <- model.matrix(~0+factor(grp))
rownames(design)<-colnames(DEG)
colnames(design)<-levels(factor(grp))
DEG <- estimateGLMCommonDisp(DEG,design)
DEG <- estimateGLMTrendedDisp(DEG, design)
DEG <- estimateGLMTagwiseDisp(DEG, design)
fit <- glmFit(DEG, design)
# edgeR User guide (Page. 30 => "GLM Approach")
lrt <- glmLRT(fit, contrast=c(-1,1)) # accoding to design to modify
# or we can use another way (glmQLFTest):
# CasevsCtrl <- makeContrasts(Case-Ctrl=trt-untrt, levels=design)
# lrt <- glmQLFTest(fit,contrast=CasevsCtrl)
nrDEG=topTags(lrt, n=nrow(DEG))
nrDEG=as.data.frame(nrDEG)
由于我们前面使用了细胞类型+批次信息两个分组条件,因此这里需要先锁定一个再分析另一个。意思就是我们可以先指定一个细胞类型,然后去进行这3个批次(各2个重复)共6个样本的差异分析,比如可以指定细胞类型是:间质细胞
label <- "Mesenchyme"
current <- summed[,label==summed$celltype.mapped]
# 于是从原来的186组样本又再次过滤为了6组样本
> dim(current)
[1] 14699 6
第一步:创建edgeR对象
library(edgeR)
y <- DGEList(counts(current), samples=colData(current))
y

第二步:预处理
首先是去除低质量样本(担心影响后面的归一化):不过这里的文库大小不是像原来一样计算每个样本的表达量加和。因为这里已经将细胞分组,看到的6不是6个细胞,而是6组细胞。因此这里的计算标准是:每组内含有不少于多少的细胞。 例如(Crowell et al. 2019)就定义了如果一组包含少于20个细胞,就是非常低的文库
discarded <- current$ncells < 20
y <- y[,!discarded]
summary(discarded)
## Mode FALSE
## logical 6
然后去除低质量基因(目的是增加后面模型构建的准确度,减少多重矫正的压力):先得到一个logCPM阈值,然后过滤。可以看到下面过滤掉了不少基因
keep <- filterByExpr(y, group=current$tomato)
y <- y[keep,]
summary(keep)
## Mode FALSE TRUE
## logical 9011 5688
最后,归一化矫正:使用TMM(trimmed mean of M-values)方法,也算得上是
calcNormFactors()
函数的标志了y <- calcNormFactors(y)

第三步:模型并分析
这部分的重点就是这个分组矩阵了,下面锁定样本中批次(pool)的差异,而保留注射导致的差异(tomato)
design <- model.matrix(~factor(pool) + factor(tomato), y$samples)
design
## (Intercept) factor(pool)4 factor(pool)5 factor(tomato)TRUE
## Sample1 1 0 0 1
## Sample2 1 0 0 0
## Sample3 1 1 0 1
## Sample4 1 1 0 0
## Sample5 1 0 1 1
## Sample6 1 0 1 0
## attr(,"assign")
## [1] 0 1 1 2
## attr(,"contrasts")
## attr(,"contrasts")$`factor(pool)`
## [1] "contr.treatment"
##
## attr(,"contrasts")$`factor(tomato)`
## [1] "contr.treatment"
接下来就是新增部分了,使用
estimateDisp()
来估计负二项分布negative binomial (NB)y <- estimateDisp(y, design)
summary(y$trended.dispersion)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0103 0.0167 0.0213 0.0202 0.0235 0.0266
接着用
glmQLFit()
估计 quasi-likelihood分布【具体的统计知识可以日后再理解,这里不是重点】fit <- glmQLFit(y, design, robust=TRUE)
用
glmQLFTest()
看一下表达量差异,差异基因是logFC不为0(logFC=1
也就是treat/control=1
)并且FDR小于5%的res <- glmQLFTest(fit, coef=ncol(design))
summary(decideTests(res))
## factor(tomato)TRUE
## Down 8
## NotSig 5672
## Up 8