OSCA单细胞数据分析笔记12—Intergrating Datasets

对应原版教程第13章http://bioconductor.org/books/release/OSCA/overview.html

无论是scRNA-seq，还是Bulk RNA-seq，批次效应都是一个很头疼的问题，如何有效地校正、并且正确地使用校正后的数据是很值得讨论的分析点。

笔记要点

1、单细胞数据的批次效应 batch effect
2、评判批次效应（示例数据）
3、两种消除批次效应的方法
4、进一步评价校正批次效应的结果
5、关于校正之后的细胞表达水平的思考

1、单细胞数据的批次效应 batch effect

当一个大型单细胞测序项目涉及多个相同来源的细胞样本时，不得不分批测序时，常常因各种各样的客观因素（操作者，实验试剂）造成不同批次间的测序结果存在一定的差异；
而这种批次间的差异是系统性的表达水平差异；如果直接合并这些多批次的数据，会引入批次效应造成的细胞表达异质性；
因此，可采用特定的计算方法对来自不同批次的细胞表达系统性差异进行校正。（传统的Bulk RNA-seq也会涉及同一意义上的批次效应，校正方法有所异同，具体见笔记后文）

2、评判批次效应（示例数据）

2.1 示例数据

来自TENxPBMCData包的两个批次的PBMC数据集

代码语言：javascript复制

#均已分别完成质控、hvg、降维、聚类。具体见原教程
pbmc3k
#class: SingleCellExperiment 
#dim: 31232 2609 
dec3k
# DataFrame with 31232 rows and 6 columns
# mean     total      tech          bio   p.value       FDR
# <numeric> <numeric> <numeric>    <numeric> <numeric> <numeric>
#   ENSG00000243485         0         0         0            0       NaN       NaN

pbmc4k
#class: SingleCellExperiment 
#dim: 31232 4182
dec4k
# DataFrame with 31232 rows and 6 columns
# mean       total        tech          bio   p.value       FDR
# <numeric>   <numeric>   <numeric>    <numeric> <numeric> <numeric>
#   ENSG00000243485 0.000000000 0.000000000 0.000000000  0.00000e 00       NaN       NaN

为了之后的批次效应评价与校正，需要进一步处理这两个批次的测序数据

（1）取测序基因交集

代码语言：javascript复制

universe <- intersect(rownames(pbmc3k), rownames(pbmc4k))
length(universe)

# Subsetting the SingleCellExperiment object.
pbmc3k <- pbmc3k[universe,]
pbmc4k <- pbmc4k[universe,]

# Also subsetting the variance modelling results, for convenience.
dec3k <- dec3k[universe,]
dec4k <- dec4k[universe,]

（2）使两个批次的测序深度相统一

multiBatchNorm()函数可校正因不同批次的测序深度不同而造成的表达水平差异，即消除了部分技术误差带来的差异。

代码语言：javascript复制

library(batchelor)
rescaled <- multiBatchNorm(pbmc3k, pbmc4k)
pbmc3k <- rescaled[[1]]
pbmc4k <- rescaled[[2]]

2.2 评判批次效应

想看看多个批次间的测序数据是否存在明显的批次效应，首先需要直接合并多个批次的数据集，然后降维聚类，观察聚类结果，是否存在仅仅由单个batch的细胞组成的cluster。若存在，则表明的确有批次效应。此外还可通过绘制t-SNE二维可视化，直接观察批次间的“分离度”

（1）合并数据集

代码语言：javascript复制

rowData(pbmc3k) <- rowData(pbmc4k)
pbmc3k$batch <- "3k"
pbmc4k$batch <- "4k"
uncorrected <- cbind(pbmc3k, pbmc4k)

（2）挑选高变基因combineVar()函数可以在考虑批次因素前提下，筛选出尽可能全面的hvg

代码语言：javascript复制

library(scran)
combined.dec <- combineVar(dec3k, dec4k)
chosen.hvgs <- combined.dec$bio > 0 #筛选标准放宽很多
sum(chosen.hvgs)
#13431

（3）PCA降维

代码语言：javascript复制

# Using RandomParam() as it is more efficient for file-backed matrices.
library(scater)
set.seed(0010101010)
uncorrected <- runPCA(uncorrected, subset_row=chosen.hvgs,
                      BSPARAM=BiocSingular::RandomParam())

（4）KNN聚类

代码语言：javascript复制

library(scran)
snn.gr <- buildSNNGraph(uncorrected, use.dimred="PCA")
clusters <- igraph::cluster_walktrap(snn.gr)$membership
tab <- table(Cluster=clusters, Batch=uncorrected$batch)
tab

基于uncorrected的聚类结果

如上图可以看到相当多的cluster的细胞仅来自batch-3k、4k
下面直接t-SNE进行可视化看一下

代码语言：javascript复制

set.seed(1111001)
uncorrected <- runTSNE(uncorrected, dimred="PCA")
plotTSNE(uncorrected, colour_by="batch")

基于uncorrected的t-SNE降维结果

如上结果，反应了来自pbmc3k与pbmc4的测序数据间存在十分明显的批次效应，接下来会介绍两种校正方法。

3、两种消除批次效应的方法

3.1 Linear regression

线性回归是去除Bulk RNA-seq的批次效应的常见思路。其原理简单理解就是以不同批次的同一基因的表达量为y，不同批次为自变量x，进行回归拟合。然后根据拟合结果，反过来计算去除批次因素后（校正到同一“批次”水平）的基因表达量。例如limma包的removeBatchEffect()函数、sva包的comBat()函数.....
对于单细胞数据的特殊性，应用线性回归需要假设两个前提:
- 首先不同批次的细胞类型组成比例(population composition)需要大致相同:
- 其次批次效应对于不同细胞类型表达水平的“回归”影响需要保持基本一致。
可以使用batchelor包的函数对多批次的单细胞数据进行线性回归校正。

This(rescaleBatches()) is roughly equivalent to applying a linear regression to the log-expression values per gene, with some adjustments to improve performance and efficiency. For each gene, the mean expression in each batch is scaled down until it is equal to the lowest mean across all batches.

代码语言：javascript复制

library(batchelor)
rescaled <- rescaleBatches(pbmc3k, pbmc4k)
rescaled
#class: SingleCellExperiment 
#dim: 31232 6791

对校正后的结果，再次进行批次效应评判

代码语言：javascript复制

# To ensure reproducibility of the randomized PCA.
set.seed(1010101010) 
rescaled <- runPCA(rescaled, subset_row=chosen.hvgs, 
                   exprs_values="corrected",
                   BSPARAM=BiocSingular::RandomParam())

snn.gr <- buildSNNGraph(rescaled, use.dimred="PCA")
clusters.resc <- igraph::cluster_walktrap(snn.gr)$membership
tab.resc <- table(Cluster=clusters.resc, Batch=rescaled$batch)
tab.resc

基于线性回归校正批次效应的聚类结果

代码语言：javascript复制

rescaled <- runTSNE(rescaled, dimred="PCA")
rescaled$batch <- factor(rescaled$batch)
plotTSNE(rescaled, colour_by="batch")

基于线性回归校正批次效应的t-SNE降维结果

如上：从聚类结果可以看出，大部分cluster的细胞组成均来自两个batch；从t-SNE可视化结果可以看出两个批次有较好的重合度。但是仍存在batch-specific cluster，提示批次效应可能未完全去除。

在评判批次效应时，要时刻记住batch-specific cluster可能是由于批次效应造成的异质性，同时也有可能是每个batch具有specific subpopulation。因此需要分析者对测序样本有一定的把握。

3.2 MNN(Mutual nearest neighbour)

关于MNN校正的核心是找到不同batch之间的cell-pairs。首先根据细胞的表达特征，计算batch1中的细胞A到batch2中所有细胞的距离，取Top k个最邻近细胞集S1；然后再计算batch2中的细胞A'到batch1中所有细胞的距离，同样取Top k个最邻近细胞集S1'。如果细胞A的S1包含细胞A'，细胞A'的S1'也包含细胞A，则将A，A'作为一对MNN pairs。将上述步骤扩展到batch的所有细胞，就可以得到一系列batch之间的MNN pairs

关于MNN pairs的理解就是，来自不同批次的同一对cell是具有相同的生物学意义。而它们之间的"距离"代表着批次效应的程度。
相比于上面提到的线性回归方法，MNN最大的优势就是不需要假设批次之间的细胞类型组成相同。MNN方法需要假设批次效应造成的误差与细胞的生物学状态无关(orthogonal to the biology)
batchelor包的fastMNN()函数：首先提供覆盖全部批次的hvg，然后进行PCA降维，用于计算不同批次细胞间距离，然后取Top k个近邻细胞，用于定位MNN pairs，最后据此去除批次效应。

代码语言：javascript复制

# Again, using randomized SVD here, as this is faster than IRLBA for
# file-backed matrices. We set deferred=TRUE for greater speed.
set.seed(1000101001)
mnn.out <- fastMNN(pbmc3k, pbmc4k, d=50, k=20, subset.row=chosen.hvgs,
                   BSPARAM=BiocSingular::RandomParam(deferred=TRUE))
mnn.out
# class: SingleCellExperiment 
# dim: 13431 6791 
# metadata(2): merge.info pca.info
# assays(1): reconstructed
# rownames(13431): ENSG00000239945 ENSG00000228463 ... ENSG00000198695 ENSG00000198727
# rowData names(1): rotation
# colnames: NULL
# colData names(1): batch
# reducedDimNames(2): corrected TSNE
# altExpNames(0):

对校正后的结果，再次进行批次效应评判

代码语言：javascript复制

library(scran)
snn.gr <- buildSNNGraph(mnn.out, use.dimred="corrected")
clusters.mnn <- igraph::cluster_walktrap(snn.gr)$membership
tab.mnn <- table(Cluster=clusters.mnn, Batch=mnn.out$batch)
tab.mnn

基于MNN校正批次效应的聚类结果

代码语言：javascript复制

library(scater)
set.seed(0010101010)
mnn.out <- runTSNE(mnn.out, dimred="corrected")
mnn.out$batch <- factor(mnn.out$batch)
plotTSNE(mnn.out, colour_by="batch")

基于MNN校正批次效应的t-SNE降维结果

如上结果，可以看出：相对线性回归，MNN具有更理想的校正结果。

最后关于K值的意义与选择：从最理想的结果来看相同组织、不同批次的相同细胞类型的细胞间最有可能组成cp(MNN pairs)；K值越大，一个细胞就会参与到更多的MNN pairs当中，即适合放宽相同细胞类型的标准的分析情况。此外如果先前已对每个批次进行单独的全套流程的单细胞数据分析，然后再尝试多个批次的合并效果。fastMNN()的subset.row参数可以是每个批次的全部marker gene的组合，详见原教程13.7，这里就不多展开了。

4、进一步评价校正批次效应的结果

以下的示例，均以3.2MNN的校正结果为例进行分析

4.1 同一cluster里不同batch来源的细胞组成比例

最理想的情况下，校正批次效应的聚类结果里，对于每一个culster的不同batch来源的细胞组成比例应当近似等于batch细胞总数的比例。举例来说，batch1有10个细胞，batch2有20个细胞，那么校正批次效应的聚类结果cluster1--4:8,cluster2--6:12是符合预期的。
但实际情况可能并不如此，可以进行卡方检验，看看哪些cluster 最不符合我们的预期(p值越小)

代码语言：javascript复制

chi.prop <- colSums(tab.mnn)/sum(tab.mnn)
chi.results <- apply(tab.mnn, 1, FUN=chisq.test, p=chi.prop)
p.values <- vapply(chi.results, "[[", i="p.value", 0)
p.values
# 1            2            3            4            5            6            7            8 
# 9.047103e-02 3.093226e-02 6.699936e-03 2.626832e-03 8.423646e-20 2.774871e-01 5.546015e-05 2.273510e-11 
# 9           10           11           12           13           14           15           16 
# 2.136201e-04 5.479656e-05 4.019086e-03 2.971868e-03 1.538104e-12 3.936414e-05 2.197300e-04 7.172327e-01

4.2 校正后的结果相比原始batch的异质性是否丢失

如果单纯为了校正批次效应而盲目合并数据，可能会丢失掉原始batch水平的异质性。可以通过比较校正前后的聚类分群结果来进行判断。
值得注意的是：当校正合并多个批次后，由于细胞数目的增加，聚类分群数也会相应的增加。例如对于batch1的细胞单独分析，可能分为10个cluster，但当合并batch2之后，可能batch1的细胞被分为了15个cluster。而我们主要就是去评价，对于batch1来说，分为10 cluster与15 cluster的相似性程度。方法即采用之前单细胞分群笔记中介绍的ARI指标。
（1）总体ARI指标值越接近1，越理想

代码语言：javascript复制

library(bluster)
ri3k <- pairwiseRand(clusters.mnn[rescaled$batch==1], colLabels(pbmc3k), mode="index")
ri3k
#0.7360574
ri4k <- pairwiseRand(clusters.mnn[rescaled$batch==2], colLabels(pbmc4k), mode="index")
ri4k
#0.8300956

（2）具体两两cluster之间的ARI 如下图结果（左下三角未绘制）：方格的颜色越深(值越大)越理想。其中对于对角线上的方格越深表示同一cluster的cell pair再聚类之后仍属于同一cluster，对于非对角线上的方格颜色越深则表示原本不属于同一cluster的cell pair重新聚类之后仍不属于同一custer。

代码语言：javascript复制

# For the first batch.
tab <- pairwiseRand(colLabels(pbmc3k), clusters.mnn[rescaled$batch==1])
heat3k <- pheatmap(tab, cluster_row=FALSE, cluster_col=FALSE,
                   col=rev(viridis::magma(100)), main="PBMC 3K probabilities", silent=TRUE)
heat3k

5、关于校正之后的细胞表达水平的思考

校正批次效应后的单细胞数据集具有批次间表达水平的统一性。由于修改了原始数据集的表达水平（log-fold change value），作者不建议使用校正后的表达值进行DEA差异分析等基于基因表达的分析方法，而是使用原始的表达水平，并设置批次参数。

This strategy is based on the expectation that any genuine DE between clusters should still be present in a within-batch comparison where batch effects are absent.

代码语言：javascript复制

# uncorrected 为直接合并的原始表达矩阵
# clusters.mnn 为基于校正批次效应合并后的表达矩阵的聚类分群结果
m.out <- findMarkers(uncorrected, clusters.mnn, block=uncorrected$batch,
                     direction="up", lfc=1, row.data=rowData(uncorrected)[,3,drop=FALSE])

# A (probably activated?) T cell subtype of some sort:
demo <- m.out[["10"]]
as.data.frame(demo[1:20,c("Symbol", "Top", "p.value", "FDR")]) 

plotExpression(uncorrected, x=I(factor(clusters.mnn)), 
               features="ENSG00000177954", colour_by="batch")   facet_wrap(~colour_by)