4.4 实战四 | 10X | 过滤后的PBMC
刘小泽写于2020.7.19

1 前言

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

准备数据

1
library(TENxPBMCData)
2
all.sce <- list(
3
pbmc3k=TENxPBMCData('pbmc3k'),
4
pbmc4k=TENxPBMCData('pbmc4k'),
5
pbmc8k=TENxPBMCData('pbmc8k')
6
)
7
8
all.sce
9
# $pbmc3k
10
# class: SingleCellExperiment
11
# dim: 32738 2700
12
# metadata(0):
13
# assays(1): counts
14
# rownames(32738): ENSG00000243485 ENSG00000237613 ...
15
# ENSG00000215616 ENSG00000215611
16
# rowData names(3): ENSEMBL_ID Symbol_TENx Symbol
17
# colnames: NULL
18
# colData names(11): Sample Barcode ... Individual
19
# Date_published
20
# reducedDimNames(0):
21
# altExpNames(0):
22
#
23
# $pbmc4k
24
# class: SingleCellExperiment
25
# dim: 33694 4340
26
# metadata(0):
27
# assays(1): counts
28
# rownames(33694): ENSG00000243485 ENSG00000237613 ...
29
# ENSG00000277475 ENSG00000268674
30
# rowData names(3): ENSEMBL_ID Symbol_TENx Symbol
31
# colnames: NULL
32
# colData names(11): Sample Barcode ... Individual
33
# Date_published
34
# reducedDimNames(0):
35
# altExpNames(0):
36
#
37
# $pbmc8k
38
# class: SingleCellExperiment
39
# dim: 33694 8381
40
# metadata(0):
41
# assays(1): counts
42
# rownames(33694): ENSG00000243485 ENSG00000237613 ...
43
# ENSG00000277475 ENSG00000268674
44
# rowData names(3): ENSEMBL_ID Symbol_TENx Symbol
45
# colnames: NULL
46
# colData names(11): Sample Barcode ... Individual
47
# Date_published
48
# reducedDimNames(0):
49
# altExpNames(0):
Copied!
多个数据集的批量处理,重点就是列表list和for循环的熟练使用,还有相关的apply家族函数。而且每个结果数据也要对应放在一个新列表中,比如下面质控使用的stats <- high.mito <- list() ,就是新建了两个空列表,然后把结果放进去

2 批量质控

数据备份

把unfiltered数据主要用在质控的探索上
1
unfiltered <- all.sce
Copied!
还是先通过线粒体含量计算质控结果,然后根据这个结果进行过滤,一个for循环搞定
1
library(scater)
2
stats <- high.mito <- list()
3
4
for (n in names(all.sce)) {
5
current <- all.sce[[n]]
6
is.mito <- grep("MT", rowData(current)$Symbol_TENx)
7
stats[[n]] <- perCellQCMetrics(current, subsets=list(Mito=is.mito))
8
high.mito[[n]] <- isOutlier(stats[[n]]$subsets_Mito_percent, type="higher")
9
all.sce[[n]] <- current[,!high.mito[[n]]]
10
}
Copied!

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

1
> lapply(high.mito, summary)
2
$pbmc3k
3
Mode FALSE TRUE
4
logical 2609 91
5
6
$pbmc4k
7
Mode FALSE TRUE
8
logical 4182 158
9
10
$pbmc8k
11
Mode FALSE TRUE
12
logical 8157 224
Copied!

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

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

3 批量归一化

这里使用的是最简单的方法logNormCounts(),就是用某个细胞中每个基因或spike-in转录本的表达量除以这个细胞计算的size factor,最后还进行一个log转换,得到一个新的矩阵:logcounts 【不过这个名字并不是很贴切,只是因为拼写简单,真实应该是:log-transformed normalized expression values。而不是单纯字面意思取个log】
1
all.sce <- lapply(all.sce, logNormCounts)
Copied!
看到这里,可能会想,为什么没计算size factor就直接进行了logNormCounts?
前面提到的操作,一般都是:
1
# 常规方法
2
lib.sf.zeisel <- librarySizeFactors(sce.zeisel)
3
lib.sf.zeisel <- logNormCounts(lib.sf.zeisel)
4
# 去卷积方法
5
clusters <- quickCluster(sce.pbmc)
6
sce.pbmc <- computeSumFactors(sce.pbmc, cluster=clusters)
7
sce.pbmc <- logNormCounts(sce.pbmc)
Copied!
看一下logNormCounts的帮助文档就能明白了,逻辑很清楚:
    函数默认的参数是:size_factors=NULL,如果没有计算size factor更新给函数,那么函数会执行normalizeCounts 的操作
    再来看normalizeCounts会有什么操作:如果没有提供size factor,它会根据数据类型去自己计算size factor
      对于count矩阵和SummarizedExperiment数据类型,会通过librarySizeFactors计算
      对于SingleCellExperiment这种数据类型,它会首先在数据中寻找size factor是否存在,如果找不到,也是会使用librarySizeFactors
也就是说,如果我们不提前计算,就会自动帮我们用最简单的librarySizeFactors计算,并添加到我们的数据中。正是因为我们只需要最简单的方法,所以才可以不提供。如果要使用去卷积方法,还是要自己先计算好
最后看下结果
1
> lapply(all.sce, function(x) summary(sizeFactors(x)))
2
$pbmc3k
3
Min. 1st Qu. Median Mean 3rd Qu. Max.
4
0.2338 0.7478 0.9262 1.0000 1.1571 6.6042
5
6
$pbmc4k
7
Min. 1st Qu. Median Mean 3rd Qu. Max.
8
0.3155 0.7109 0.8903 1.0000 1.1272 11.0267
9
10
$pbmc8k
11
Min. 1st Qu. Median Mean 3rd Qu. Max.
12
0.2963 0.7043 0.8772 1.0000 1.1177 6.7942
Copied!

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

同样也是使用了最简单的计算方法
1
library(scran)
2
all.dec <- lapply(all.sce, modelGeneVar)
3
all.hvgs <- lapply(all.dec, getTopHVGs, prop=0.1)
Copied!
作图
1
par(mfrow=c(1,3))
2
for (n in names(all.dec)) {
3
curdec <- all.dec[[n]]
4
plot(curdec$mean, curdec$total, pch=16, cex=0.5, main=n,
5
xlab="Mean of log-expression", ylab="Variance of log-expression")
6
curfit <- metadata(curdec)
7
# 可视化一条线(下图的蓝线),这条线指所有的基因都会存在的一种偏差
8
curve(curfit$trend(x), col='dodgerblue', add=TRUE, lwd=2)
9
}
Copied!
    每个点表示一个基因
    图中蓝线指的是:技术因素导致的偏差
    纵坐标表示总偏差:它等于技术偏差+生物因素偏差
因此,要衡量一个基因的生物因素偏差大小,就看对应的纵坐标减去对应的蓝线的值

5 批量降维

这里将每个PBMC数据单独进行降维,而并没有把它们混合起来再分析
关于降维方法,这里选择是与PCA近似的SVD算法(singular value decomposition,奇异值分解),scater或scran都可以直接通过函数计算SVD,利用参数BSPARAM=传递一个BiocSingularParam对象到runPCA
    SVD是一种矩阵分解方法,相当于因式分解,目的纯粹就是将一个矩阵拆分成多个矩阵相乘的形式
    对于稀疏矩阵来说,SVD算法更适用,这样对于大数据来说节省了很大空间
1
library(BiocSingular)
2
set.seed(10000)
3
# 这里的runPCA是一个降维的函数名字,它可以用本身的PCA算法,还可以用SVD算法
4
all.sce <- mapply(FUN=runPCA, x=all.sce, subset_row=all.hvgs,
5
MoreArgs=list(ncomponents=25, BSPARAM=RandomParam()),
6
SIMPLIFY=FALSE)
7
8
set.seed(100000)
9
all.sce <- lapply(all.sce, runTSNE, dimred="PCA")
10
11
set.seed(1000000)
12
all.sce <- lapply(all.sce, runUMAP, dimred="PCA")
Copied!
这里使用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”。
1
for (n in names(all.sce)) {
2
g <- buildSNNGraph(all.sce[[n]], k=10, use.dimred='PCA')
3
clust <- igraph::cluster_walktrap(g)$membership
4
colLabels(all.sce[[n]]) <- factor(clust)
5
}
6
# 看看各自分了多少群
7
lapply(all.sce, function(x) table(colLabels(x)))
8
## $pbmc3k
9
##
10
## 1 2 3 4 5 6 7 8 9 10
11
## 487 154 603 514 31 150 179 333 147 11
12
##
13
## $pbmc4k
14
##
15
## 1 2 3 4 5 6 7 8 9 10 11 12 13
16
## 497 185 569 786 373 232 44 1023 77 218 88 54 36
17
##
18
## $pbmc8k
19
##
20
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
21
## 1004 759 1073 1543 367 150 201 2067 59 154 244 67 76 285 20 15
22
## 17 18
23
## 64 9
Copied!
作图
1
# 还是先新建一个空列表,方便保存数据
2
all.tsne <- list()
3
4
for (n in names(all.sce)) {
5
all.tsne[[n]] <- plotTSNE(all.sce[[n]], colour_by="label") + ggtitle(n)
6
}
7
do.call(gridExtra::grid.arrange, c(all.tsne, list(ncol=3)))
Copied!

7 数据整合

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

找共有基因

1
# 先看看各自多少基因
2
> lapply(all.sce, function(x) length(rownames(x)))
3
$pbmc3k
4
[1] 32738
5
6
$pbmc4k
7
[1] 33694
8
9
$pbmc8k
10
[1] 33694
11
12
# 找共同基因
13
universe <- Reduce(intersect, lapply(all.sce, rownames))
14
> length(universe)
15
[1] 31232
Copied!

对每个数据批量取子集

1
# 这个操作可以记下来,lapply批量取子集
2
all.sce2 <- lapply(all.sce, "[", i=universe,)
3
4
> lapply(all.sce2, dim)
5
$pbmc3k
6
[1] 31232 2609
7
8
$pbmc4k
9
[1] 31232 4182
10
11
$pbmc8k
12
[1] 31232 8157
13
14
# 同样的,对找高变异基因的结果取子集
15
all.dec2 <- lapply(all.dec, "[", i=universe,)
Copied!

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

1
library(batchelor)
2
normed.sce <- do.call(multiBatchNorm, all.sce2)
Copied!

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

这次是把3个批次放在一起再找的
1
combined.dec <- do.call(combineVar, all.dec2)
2
combined.hvg <- getTopHVGs(combined.dec, n=5000)
Copied!

对一个大数据进行降维

结合:单细胞交响乐10-数据集整合后的批次矫正 中的【第4部分 MNN矫正】
fastMNN先进行PCA降维,以加速下面的聚类环节
1
set.seed(1000101)
2
merged.pbmc <- do.call(fastMNN, c(normed.sce,
3
list(subset.row=combined.hvg, BSPARAM=RandomParam())))
4
5
merged.pbmc
6
# class: SingleCellExperiment
7
# dim: 5000 14948
8
# metadata(2): merge.info pca.info
9
# assays(1): reconstructed
10
# rownames(5000): ENSG00000090382 ENSG00000163220 ...
11
# ENSG00000122068 ENSG00000011132
12
# rowData names(1): rotation
13
# colnames: NULL
14
# colData names(2): batch label
15
# reducedDimNames(1): corrected
16
# altExpNames(0):
17
18
# 这时就出现了批次信息
19
> table(merged.pbmc$batch)
20
pbmc3k pbmc4k pbmc8k
21
2609 4182 8157
Copied!
检查一下结果,使用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.
1
metadata(merged.pbmc)$merge.info$lost.var
2
## pbmc3k pbmc4k pbmc8k
3
## [1,] 7.003e-03 3.126e-03 0.000000
4
## [2,] 7.137e-05 5.125e-05 0.003003
Copied!

对一个大数据进行聚类

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

可视化

1
set.seed(10101010)
2
merged.pbmc <- runTSNE(merged.pbmc, dimred="corrected")
3
# 查看3个数据混合后的聚类结果,以及有没有批次效应
4
gridExtra::grid.arrange(
5
plotTSNE(merged.pbmc, colour_by="label", text_by="label", text_colour="red"),
6
plotTSNE(merged.pbmc, colour_by="batch")
7
)
Copied!
最近更新 1yr ago