单细胞交响乐21-实战四批量处理并整合多个10X PBMC数据

刘小泽写于2020.7.19
为何取名叫“交响乐”？因为单细胞分析就像一个大乐团，需要各个流程的协同配合
单细胞交响乐1-常用的数据结构SingleCellExperiment
单细胞交响乐2-scRNAseq从实验到下游简介
 单细胞交响乐3-细胞质控
 单细胞交响乐4-归一化
 单细胞交响乐5-挑选高变化基因
 单细胞交响乐6-降维
 单细胞交响乐7-聚类分群
 单细胞交响乐8-marker基因检测
 单细胞交响乐9-细胞类型注释
 单细胞交响乐9-细胞类型注释
 单细胞交响乐10-数据集整合后的批次矫正
 单细胞交响乐11-多样本间差异分析
 单细胞交响乐12-检测Doublet
单细胞交响乐13-细胞周期推断
 单细胞交响乐14-细胞轨迹推断
 单细胞交响乐15-scRNA与蛋白丰度信息结合
 单细胞交响乐16-处理大型数据
 单细胞交响乐17-不同单细胞R包的数据格式相互转换
 单细胞交响乐18-实战一 Smart-seq2
单细胞交响乐19-实战二 STRT-Seq
单细胞交响乐20-实战三 10X 未过滤的PBMC数据

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与PCA：https://www.cnblogs.com/bjwu/p/9280492.html
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")
)

欢迎关注我们的公众号～_～　　
我们是两个农转生信的小硕，打造生信星球，想让它成为一个不拽术语、通俗易懂的生信知识平台。需要帮助或提出意见请后台留言或发送邮件到jieandze1314@gmail.com

Welcome to our bioinfoplanet!

单细胞交响乐21-实战四 批量处理并整合多个10X PBMC数据