4.4 实战四 | 10X | 过滤后的PBMC

刘小泽写于2020.7.19

1 前言

这次使用的数据是来自Zheng et al. 2017 的三个PBMC数据,而且这个数据是经过前期过滤的

准备数据

library(TENxPBMCData)
all.sce <- list(
pbmc3k=TENxPBMCData('pbmc3k'),
pbmc4k=TENxPBMCData('pbmc4k'),
pbmc8k=TENxPBMCData('pbmc8k')
)
all.sce
# $pbmc3k
# class: SingleCellExperiment
# dim: 32738 2700
# metadata(0):
# assays(1): counts
# rownames(32738): ENSG00000243485 ENSG00000237613 ...
# ENSG00000215616 ENSG00000215611
# rowData names(3): ENSEMBL_ID Symbol_TENx Symbol
# colnames: NULL
# colData names(11): Sample Barcode ... Individual
# Date_published
# reducedDimNames(0):
# altExpNames(0):
#
# $pbmc4k
# class: SingleCellExperiment
# dim: 33694 4340
# metadata(0):
# assays(1): counts
# rownames(33694): ENSG00000243485 ENSG00000237613 ...
# ENSG00000277475 ENSG00000268674
# rowData names(3): ENSEMBL_ID Symbol_TENx Symbol
# colnames: NULL
# colData names(11): Sample Barcode ... Individual
# Date_published
# reducedDimNames(0):
# altExpNames(0):
#
# $pbmc8k
# class: SingleCellExperiment
# dim: 33694 8381
# metadata(0):
# assays(1): counts
# rownames(33694): ENSG00000243485 ENSG00000237613 ...
# ENSG00000277475 ENSG00000268674
# rowData names(3): ENSEMBL_ID Symbol_TENx Symbol
# colnames: NULL
# colData names(11): Sample Barcode ... Individual
# Date_published
# reducedDimNames(0):
# altExpNames(0):

多个数据集的批量处理,重点就是列表list和for循环的熟练使用,还有相关的apply家族函数。而且每个结果数据也要对应放在一个新列表中,比如下面质控使用的stats <- high.mito <- list() ,就是新建了两个空列表,然后把结果放进去

2 批量质控

数据备份

把unfiltered数据主要用在质控的探索上

unfiltered <- all.sce

还是先通过线粒体含量计算质控结果,然后根据这个结果进行过滤,一个for循环搞定

library(scater)
stats <- high.mito <- list()
for (n in names(all.sce)) {
current <- all.sce[[n]]
is.mito <- grep("MT", rowData(current)$Symbol_TENx)
stats[[n]] <- perCellQCMetrics(current, subsets=list(Mito=is.mito))
high.mito[[n]] <- isOutlier(stats[[n]]$subsets_Mito_percent, type="higher")
all.sce[[n]] <- current[,!high.mito[[n]]]
}

看一下根据线粒体过滤的结果

> lapply(high.mito, summary)
$pbmc3k
Mode FALSE TRUE
logical 2609 91
$pbmc4k
Mode FALSE TRUE
logical 4182 158
$pbmc8k
Mode FALSE TRUE
logical 8157 224

批量作图(也是把作图结果放进list,方便后期批量导出)

qcplots <- list()
for (n in names(all.sce)) {
current <- unfiltered[[n]]
colData(current) <- cbind(colData(current), stats[[n]])
current$discard <- high.mito[[n]]
qcplots[[n]] <- plotColData(current, x="sum", y="subsets_Mito_percent",
colour_by="discard") + scale_x_log10()
}
# do.call也是list处理中的常用函数
do.call(gridExtra::grid.arrange, c(qcplots, ncol=3))

3 批量归一化

这里使用的是最简单的方法logNormCounts(),就是用某个细胞中每个基因或spike-in转录本的表达量除以这个细胞计算的size factor,最后还进行一个log转换,得到一个新的矩阵:logcounts 【不过这个名字并不是很贴切,只是因为拼写简单,真实应该是:log-transformed normalized expression values。而不是单纯字面意思取个log】

all.sce <- lapply(all.sce, logNormCounts)

看到这里,可能会想,为什么没计算size factor就直接进行了logNormCounts?

前面提到的操作,一般都是:

# 常规方法
lib.sf.zeisel <- librarySizeFactors(sce.zeisel)
lib.sf.zeisel <- logNormCounts(lib.sf.zeisel)
# 去卷积方法
clusters <- quickCluster(sce.pbmc)
sce.pbmc <- computeSumFactors(sce.pbmc, cluster=clusters)
sce.pbmc <- logNormCounts(sce.pbmc)

看一下logNormCounts的帮助文档就能明白了,逻辑很清楚:

  • 函数默认的参数是:size_factors=NULL,如果没有计算size factor更新给函数,那么函数会执行normalizeCounts 的操作

  • 再来看normalizeCounts会有什么操作:如果没有提供size factor,它会根据数据类型去自己计算size factor

    • 对于count矩阵和SummarizedExperiment数据类型,会通过librarySizeFactors计算

    • 对于SingleCellExperiment这种数据类型,它会首先在数据中寻找size factor是否存在,如果找不到,也是会使用librarySizeFactors

也就是说,如果我们不提前计算,就会自动帮我们用最简单的librarySizeFactors计算,并添加到我们的数据中。正是因为我们只需要最简单的方法,所以才可以不提供。如果要使用去卷积方法,还是要自己先计算好

最后看下结果

> lapply(all.sce, function(x) summary(sizeFactors(x)))
$pbmc3k
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.2338 0.7478 0.9262 1.0000 1.1571 6.6042
$pbmc4k
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.3155 0.7109 0.8903 1.0000 1.1272 11.0267
$pbmc8k
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.2963 0.7043 0.8772 1.0000 1.1177 6.7942

4 批量找表达量高变化基因

同样也是使用了最简单的计算方法

library(scran)
all.dec <- lapply(all.sce, modelGeneVar)
all.hvgs <- lapply(all.dec, getTopHVGs, prop=0.1)

作图

par(mfrow=c(1,3))
for (n in names(all.dec)) {
curdec <- all.dec[[n]]
plot(curdec$mean, curdec$total, pch=16, cex=0.5, main=n,
xlab="Mean of log-expression", ylab="Variance of log-expression")
curfit <- metadata(curdec)
# 可视化一条线(下图的蓝线),这条线指所有的基因都会存在的一种偏差
curve(curfit$trend(x), col='dodgerblue', add=TRUE, lwd=2)
}
  • 每个点表示一个基因

  • 图中蓝线指的是:技术因素导致的偏差

  • 纵坐标表示总偏差:它等于技术偏差+生物因素偏差

因此,要衡量一个基因的生物因素偏差大小,就看对应的纵坐标减去对应的蓝线的值

5 批量降维

这里将每个PBMC数据单独进行降维,而并没有把它们混合起来再分析

关于降维方法,这里选择是与PCA近似的SVD算法(singular value decomposition,奇异值分解),scater或scran都可以直接通过函数计算SVD,利用参数BSPARAM=传递一个BiocSingularParam对象到runPCA

  • SVD是一种矩阵分解方法,相当于因式分解,目的纯粹就是将一个矩阵拆分成多个矩阵相乘的形式

  • 对于稀疏矩阵来说,SVD算法更适用,这样对于大数据来说节省了很大空间

library(BiocSingular)
set.seed(10000)
# 这里的runPCA是一个降维的函数名字,它可以用本身的PCA算法,还可以用SVD算法
all.sce <- mapply(FUN=runPCA, x=all.sce, subset_row=all.hvgs,
MoreArgs=list(ncomponents=25, BSPARAM=RandomParam()),
SIMPLIFY=FALSE)
set.seed(100000)
all.sce <- lapply(all.sce, runTSNE, dimred="PCA")
set.seed(1000000)
all.sce <- lapply(all.sce, runUMAP, dimred="PCA")

这里使用SVD其实还是为了帮助更好地进行PCA,支持4种方式:

  • ExactParam: exact SVD with runExactSVD.

  • IrlbaParam: approximate SVD with irlba via runIrlbaSVD.

  • RandomParam: approximate SVD with rsvd via runRandomSVD.

  • FastAutoParam: fast approximate SVD, chosen based on the matrix representation.

6 批量聚类

使用基于图形的聚类,最基础的想法是:我们首先构建一个图,其中每个节点都是一个细胞,它与高维空间中最邻近的细胞相连。连线基于细胞之间的相似性计算权重,权重越高,表示细胞间关系更密切。如果一群细胞之间的权重高于另一群细胞,那么这一群细胞就被当做一个群体 “communiity”。

for (n in names(all.sce)) {
g <- buildSNNGraph(all.sce[[n]], k=10, use.dimred='PCA')
clust <- igraph::cluster_walktrap(g)$membership
colLabels(all.sce[[n]]) <- factor(clust)
}
# 看看各自分了多少群
lapply(all.sce, function(x) table(colLabels(x)))
## $pbmc3k
##
## 1 2 3 4 5 6 7 8 9 10
## 487 154 603 514 31 150 179 333 147 11
##
## $pbmc4k
##
## 1 2 3 4 5 6 7 8 9 10 11 12 13
## 497 185 569 786 373 232 44 1023 77 218 88 54 36
##
## $pbmc8k
##
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
## 1004 759 1073 1543 367 150 201 2067 59 154 244 67 76 285 20 15
## 17 18
## 64 9

作图

# 还是先新建一个空列表,方便保存数据
all.tsne <- list()
for (n in names(all.sce)) {
all.tsne[[n]] <- plotTSNE(all.sce[[n]], colour_by="label") + ggtitle(n)
}
do.call(gridExtra::grid.arrange, c(all.tsne, list(ncol=3)))

7 数据整合

前面只是批量进行了各个数据集的分析,现在要把它们整合起来再分析一下

找共有基因

# 先看看各自多少基因
> lapply(all.sce, function(x) length(rownames(x)))
$pbmc3k
[1] 32738
$pbmc4k
[1] 33694
$pbmc8k
[1] 33694
# 找共同基因
universe <- Reduce(intersect, lapply(all.sce, rownames))
> length(universe)
[1] 31232

对每个数据批量取子集

# 这个操作可以记下来,lapply批量取子集
all.sce2 <- lapply(all.sce, "[", i=universe,)
> lapply(all.sce2, dim)
$pbmc3k
[1] 31232 2609
$pbmc4k
[1] 31232 4182
$pbmc8k
[1] 31232 8157
# 同样的,对找高变异基因的结果取子集
all.dec2 <- lapply(all.dec, "[", i=universe,)

把三个数据当做一个数据的三个批次,重新进行归一化

library(batchelor)
normed.sce <- do.call(multiBatchNorm, all.sce2)

根据重新归一化的结果,再次找HVGs

这次是把3个批次放在一起再找的

combined.dec <- do.call(combineVar, all.dec2)
combined.hvg <- getTopHVGs(combined.dec, n=5000)

对一个大数据进行降维

结合:单细胞交响乐10-数据集整合后的批次矫正 中的【第4部分 MNN矫正】

fastMNN先进行PCA降维,以加速下面的聚类环节

set.seed(1000101)
merged.pbmc <- do.call(fastMNN, c(normed.sce,
list(subset.row=combined.hvg, BSPARAM=RandomParam())))
merged.pbmc
# class: SingleCellExperiment
# dim: 5000 14948
# metadata(2): merge.info pca.info
# assays(1): reconstructed
# rownames(5000): ENSG00000090382 ENSG00000163220 ...
# ENSG00000122068 ENSG00000011132
# rowData names(1): rotation
# colnames: NULL
# colData names(2): batch label
# reducedDimNames(1): corrected
# altExpNames(0):
# 这时就出现了批次信息
> table(merged.pbmc$batch)
pbmc3k pbmc4k pbmc8k
2609 4182 8157

检查一下结果,使用lost.var ,值越大表示丢失的真实生物异质性越多

  • It contains a matrix of the variance lost in each batch (column) at each merge step (row).

  • Large proportions of lost variance (>10%) suggest that correction is removing genuine biological heterogeneity.

metadata(merged.pbmc)$merge.info$lost.var
## pbmc3k pbmc4k pbmc8k
## [1,] 7.003e-03 3.126e-03 0.000000
## [2,] 7.137e-05 5.125e-05 0.003003

对一个大数据进行聚类

g <- buildSNNGraph(merged.pbmc, use.dimred="corrected")
colLabels(merged.pbmc) <- factor(igraph::cluster_louvain(g)$membership)
table(colLabels(merged.pbmc), merged.pbmc$batch)
##
## pbmc3k pbmc4k pbmc8k
## 1 113 387 825
## 2 507 395 806
## 3 175 344 581
## 4 295 539 1018
## 5 346 638 1210
## 6 11 3 9
## 7 17 27 111
## 8 33 113 185
## 9 423 754 1546
## 10 4 36 67
## 11 197 124 221
## 12 150 180 293
## 13 327 588 1125
## 14 11 54 160

可视化

set.seed(10101010)
merged.pbmc <- runTSNE(merged.pbmc, dimred="corrected")
# 查看3个数据混合后的聚类结果,以及有没有批次效应
gridExtra::grid.arrange(
plotTSNE(merged.pbmc, colour_by="label", text_by="label", text_colour="red"),
plotTSNE(merged.pbmc, colour_by="batch")
)