3.8 批次效应处理

刘小泽写于2020.7.7

1 前言

大型的scRNA分析一般都需要整合多个批次的数据集，但每个批次的质量可能存在不可调和的差异，例如仪器的调整、试剂质量的差异。结果就导致了不同批次细胞中基因表达量的系统差异，称之为“批次效应”。之前在 3.2 归一化 这一部分中也介绍了一些批次效应的事情，它最大的隐患就是：如果它这个差异成为整个数据的主力军，那么真正的生物学差异就难以检测，对结果的解读变得更难。

因此现在出了一些多数据整合时批次效应处理工具，早期的方法 (Ritchie et al. 2015; Leek et al. 2012)是基于线性模型，它假设细胞群体的组成要么是已知的，要么是在批次间也是同类型的。但很明显这个假设有局限，我们不能保证这个假设成立，于是定制的算法产生(Haghverdi et al. 2018; Butler et al. 2018; Lin et al. 2019)，它不再需要对细胞组成有一个先验知识，保证了对未知的探索过程。

数据准备

将使用两个不同批次的10X PBMC数据进行整合，它们都是从 TENxPBMCData这个数据包获得的，分别进行了前期的数据处理。各自前期处理的一个好处，就是可以根据各自的特点先过滤掉一些低质量数据（比如各自根据QC结果对低质量细胞进行过滤）

常规的前期处理步骤还是：质控=》归一化=》找HVGs=》降维=》聚类 只是注意下面的lapply的使用，多个数据的批量处理用列表很高效

# 数据下载（数据我已做好，链接在：https://share.weiyun.com/iE0VjVD0）
library(TENxPBMCData)
all.sce <- list(
    pbmc3k=TENxPBMCData('pbmc3k'),
    pbmc4k=TENxPBMCData('pbmc4k'),
    pbmc8k=TENxPBMCData('pbmc8k')
)

# 批量质控
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]]]
}

# 批量归一化
all.sce <- lapply(all.sce, logNormCounts)

# 批量找HVGs
library(scran)
all.dec <- lapply(all.sce, modelGeneVar)
all.hvgs <- lapply(all.dec, getTopHVGs, prop=0.1)
> length(all.hvgs[[1]])
[1] 1001
> length(all.hvgs[[2]])
[1] 1352

# 批量降维(使用了一个更高级的mapply：Apply a Function to Multiple List or Vector Arguments)
library(BiocSingular)
set.seed(10000)
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")

# 批量聚类
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)
}

# 把数据提取出来
pbmc3k <- all.sce$pbmc3k
dec3k <- all.dec$pbmc3k
pbmc4k <- all.sce$pbmc4k
dec4k <- all.dec$pbmc4k

为了方便下面的批次矫正，先做几件事

1 选出两个数据集的基因交集，并各自取子集

这里两个数据都是Ensembl ID，因此这个过程很简单

universe <- intersect(rownames(pbmc3k), rownames(pbmc4k))
> length(rownames(pbmc3k));length(rownames(pbmc4k));length(universe)
[1] 32738
[1] 33694
[1] 31232

然后取子集

pbmc3k <- pbmc3k[universe,]
pbmc4k <- pbmc4k[universe,]

dec3k <- dec3k[universe,]
dec4k <- dec4k[universe,]

2 矫正批次间测序深度的影响

利用multiBatchNorm()函数，将两个批次放在一起，重新计算log-count，可以提高下游批次矫正的质量

还记得之前的size factor这个名词吗？ 它是为了去除一个批次数据中的各个细胞之间的文库差异 而这里，是为了去除两个批次之间的差异

library(batchelor)
rescaled <- multiBatchNorm(pbmc3k, pbmc4k)
> class(rescaled)
[1] "list"
pbmc3k <- rescaled[[1]]
pbmc4k <- rescaled[[2]]

3 将之前单独挑出来的模型构建结果再次整合

当然还也可以把每个批次的HVGs取交集或并集，选多少基因这个没有成文规定。

这里只是因为scran提供了这么一个函数，并可以提供更多一些的基因，毕竟现在基因多一点没坏处，省的后面捉襟见肘

library(scran)
combined.dec <- combineVar(dec3k, dec4k)
chosen.hvgs <- combined.dec$bio > 0
sum(chosen.hvgs)
## [1] 13431

2 检查批次效应

在进行批次矫正之前，肯定要先看看有没有批次效应以及有多严重吧

最简单的方法是把两个数据合并，然后PCA+t-SNE

合并数据

如果简单的cbind(pbmc3k, pbmc4k)，会报错

> cbind(pbmc3k, pbmc4k)
Error in FUN(X[[i]], ...) : 
  column(s) 'Symbol_TENx' in ‘mcols’ are duplicated and the data do not match
# 看似问题出现在mcols中的Symbol_TENx，看一眼

看样子的确存在不匹配的问题，第一行就暴露了，同一个Ensembl ID的Symbol就不同。这个问题先不深究。先看看两个数据的Ensembl 是否一致吧

> identical(mcols(pbmc3k)[,1],mcols(pbmc4k)[,1])
[1] TRUE
# 一致，那么就把它们的变成一样就好了，symbol的问题我们不考虑
rowData(pbmc3k) <- rowData(pbmc4k)
pbmc3k$batch <- "3k"
pbmc4k$batch <- "4k"
uncorrected <- cbind(pbmc3k, pbmc4k)

进行PCA

使用的HVGs也是合并两个数据集后得到的chosen.hvgs

library(scater)
set.seed(0010101010)
uncorrected <- runPCA(uncorrected, subset_row=chosen.hvgs,
    BSPARAM=BiocSingular::RandomParam())

再看聚类分群结果

如果说，这两个批次是真的重复，差异不大的话，那么分群的cluster应该在两个批次中包含差不多数量的细胞。但是下面发现基本clusters主要是仅仅包含了一个批次

library(scran)
snn.gr <- buildSNNGraph(uncorrected, use.dimred="PCA")
clusters <- igraph::cluster_walktrap(snn.gr)$membership
tab <- table(Cluster=clusters, Batch=uncorrected$batch)
tab
##        Batch
## Cluster   3k   4k
##      1     1  781
##      2     0 1309
##      3     0  535
##      4    14   51
##      5     0  605
##      6   489    0
##      7     0  184
##      8  1272    0
##      9     0  414
##      10  151    0
##      11    0   50
##      12  155    0
##      13    0   65
##      14    0   61
##      15    0   88
##      16   30    0
##      17  339    0
##      18  145    0
##      19   11    3
##      20    2   36

看t-SNE结果

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

虽然说t-SNE图上的点大小、远近说明不了问题，但是两个批次分的这么开，也提供了两个批次的证据

当然，还可以继续找证据。利用上一次的细胞类型注释，看看两个批次中的细胞类型，是不是差异很大。

3 矫正批次效应之线性回归

bulk mRNA转录组中常用的矫正批次效应方法就是线性回归，对每个基因表达量拟合一个线性模型。例如limma的removeBatchEffect() (Ritchie et al. 2015) 、sva的comBat() ( Leek et al. 2012)。如果要使用这类方法，就需要假设：批次间的细胞组成相同。另外的一个假设是：批次效应的累积的，对于任何给定的基因，在不同亚群中经过任何因素诱导的表达变化倍数是相同的。（其实，从这两个假设就看出来，这个方法不适合我们的单细胞数据，但还是要继续了解下去）

下面将使用batchelor 包中的rescaleBatches()函数进行处理

它也是对每个基因的log表达量进行了线性回归，并提高了一些运行性能。另外与removeBatchEffect()不同的是，rescaleBatches()保持了数据的稀疏性，而removeBatchEffect()会破坏稀疏性

library(batchelor)
rescaled <- rescaleBatches(pbmc3k, pbmc4k)
# 依然进行PCA、聚类
set.seed(1010101010) 
rescaled <- runPCA(rescaled, subset_row=chosen.hvgs, exprs_values="corrected")

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
##        Batch
## Cluster    1    2
##      1   278  525
##      2    16   23
##      3   337  606
##      4    43  748
##      5   604  529
##      6    22   71
##      7   188   48
##      8    25   49
##      9   263    0
##      10  123  135
##      11   16   85
##      12   11   57
##      13  116    6
##      14  455 1035
##      15    6   31
##      16   89  187
##      17    3   36
##      18    3    8
##      19   11    3

效果还是有的，现在大部分clusters都包含了两个批次的细胞，但是依然有存在批次差异的cluster（看cluster9，在batch2中细胞数为0），表明处理有效果但不彻底 。影响效果的原因是：我们的数据违背了这个方法的假设

再做个图看看

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

注意

除了这个函数，还可以尝试更常规的regressBatches() ，只不过它的表现可能会更差，因为会丢失数据稀疏性

4 MNN矫正

4.1 先了解一下这个方法

试想batchA中有一个细胞a，然后想根据挑选的特征基因（如HVGs）去在batchB中鉴定与a细胞空间相近的细胞们。同样的，对batchB中的细胞b 重复这个过程，也鉴定它在batchA中的近邻。于是这个MNN就是：Mutual nearest neighbors ，是对一个批次的细胞找到它们在另一个批次中对应的最相邻细胞集合。

Mutual nearest neighbors are pairs of cells from different batches that belong in each other s set of nearest neighbors. Haghverdi et al. (2018)

在进行批次矫正前，MNN找到的一对对细胞都是生物学状态相似的，因此如果再有不同，那么就姑且认为是外部的批次效应导致的，也就方便了去除。

与线性回归方法相比，MNN方法不会假设细胞群组成相同或者事先已知。MNN会自己学习细胞群的结构并进行估计。

4.2 还是应用到PBMC数据

batchelor 这个包除了线性模型，还提供了MNN 的函数fastMNN()，但这个函数与单纯的MNN算法还有不同。fastMNN()先进行PCA降维，以加速下面的检测细胞近邻。

set.seed(1000101001)
# 这里的d指的是前多少个主成分，k指近邻的数量（k增大会导致数据集合并更激进；一般不要设置太大，如果看到相同的细胞类型没有在批次之间充分合并，可以适当增加k）
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(1): corrected
## altExpNames(0):

其中mnn.out的每列表示某个批次中的一个细胞（具体存储在batch），行表示chosen.hvgs中的基因

head(mnn.out$batch) 
## [1] 1 1 1 1 1 1

看到降维相关模块reducedDimNames中多了一个corrected，它是根据50个主成分得到的所有6791在低维空间的坐标

dim(reducedDim(mnn.out, "corrected"))
## [1] 6791   50

> reducedDim(mnn.out, "corrected")[1:3,1:3]
            [,1]       [,2]         [,3]
[1,] -0.12783777  0.0409469 -0.000116194
[2,] -0.03364614 -0.1466529  0.161690348
[3,] -0.09895631  0.1219669  0.010321992

看到表达量模块assays() 中多了一个reconstructed矩阵，是每个细胞中每个基因矫正后的表达量

> dim(assay(mnn.out, "reconstructed"))
[1] 13431  6791

> assay(mnn.out, "reconstructed")[1:4,1:4]
<4 x 4> matrix of class LowRankMatrix and type "double":
                         [,1]          [,2]
ENSG00000239945 -2.522191e-06 -1.851424e-06
ENSG00000228463 -6.626821e-04 -6.724341e-04
ENSG00000237094 -8.077231e-05 -8.038006e-05
ENSG00000229905  3.838135e-06  6.179994e-06
                         [,3]          [,4]
ENSG00000239945 -1.198984e-05 -3.192269e-06
ENSG00000228463 -4.820230e-04 -6.330560e-05
ENSG00000237094 -9.630608e-05 -6.855333e-05
ENSG00000229905  5.432122e-06 -2.118592e-05

4.3 矫正后的检查

这次函数已经做过PCA了，那就继续分群

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
##        Batch
## Cluster    1    2
##      1   337  606
##      2   289  542
##      3   152  181
##      4    12    4
##      5   517  467
##      6    17   19
##      7   313  661
##      8   162  118
##      9    11   56
##      10  547 1083
##      11   17   59
##      12   16   58
##      13  144   93
##      14   67  191
##      15    4   36
##      16    4    8

可以看到这里的效果不错，没有哪个cluster在某个batch中细胞数为0

再作图看看

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")

注意

fastMNN()函数其实也自带了检测数据

metadata(mnn.out)$merge.info$lost.var
##             [,1]        [,2]
## [1,] 0.006617087 0.003315395

它的含义是：the proportion of variance within each batch that is lost during MNN correction。如果这个值很大（比如>10%），说明有些”矫枉过正“，把一些内在的生物学差异也给矫正了。

5 更为细致的检查

5.1 同一批次内clusters的比较

两个批次数据混合后进行矫正，还是应该保留每个批次的特性。一个很不想看到的情况是：看上去细胞混合效果很好，以为矫正的效果不错，但没想到矫正过了，细胞中生物学异质性也被抹除了。

我们想在整合+矫正批次后，看到矫正后的一个cluster中既包含了原来的batch1信息，也包含batch2。当然表现的差异越大越好，体现在下图中就是：after1（即矫正后的cluster1）对应的PBMC3k 的clusters与PBMC 4k的clusters最好别重复，并且差别越大越好

library(pheatmap)
# 画出二者混合后的cluster与各自之前的cluster对比
# batch 1
tab <- table(paste("after", clusters.mnn[rescaled$batch==1]),
    paste("before", colLabels(pbmc3k)))
heat3k <- pheatmap(log10(tab+10), cluster_row=FALSE, cluster_col=FALSE,
    main="PBMC 3K comparison", silent=TRUE)

# batch 2
tab <- table(paste("after", clusters.mnn[rescaled$batch==2]),
    paste("before", colLabels(pbmc4k)))
heat4k <- pheatmap(log10(tab+10), cluster_row=FALSE, cluster_col=FALSE,
    main="PBMC 4K comparison", silent=TRUE)

gridExtra::grid.arrange(heat3k[[4]], heat4k[[4]])

就以after5来说，分别对应3k混合之前的cluster3、4、7；以及4k混合前的cluster1、10、3，差别蛮大，说明混合后依然保持了各自的个性

除了看图，还能通过数据体现——计算Rand index

这个指标是矫正后各个batch原来差异的保留度，这个数值越大越好（但如果某个批次矫正方法不给力，这个值也是大的，所以还是要慎重对待）

library(fossil)
# batch1
ri3k <- rand.index(as.integer(clusters.mnn[rescaled$batch==1]),
    as.integer(colLabels(pbmc3k)))
ri3k
## [1] 0.9335226

# batch2
ri4k <- rand.index(as.integer(clusters.mnn[rescaled$batch==2]),
    as.integer(colLabels(pbmc4k)))
ri4k
## [1] 0.9575746

5.2 利用marker基因检查

目的就是看数据整合+矫正后，原始batch的内部结构是否还完整

首先分别从原来的两个批次数据中寻找marker基因，因为marker基因的存在就保证了这个批次数据集中”生物学的有趣性“存在

stats3 <- pairwiseWilcox(pbmc3k, direction="up")
markers3 <- getTopMarkers(stats3[[1]], stats3[[2]], n=10)

stats4 <- pairwiseWilcox(pbmc4k, direction="up")
markers4 <- getTopMarkers(stats4[[1]], stats4[[2]], n=10)

然后用双方的marker基因合集，作为矫正的HVGs，看看算法会不会把这部分”有趣性“去除

marker.set <- unique(unlist(c(unlist(markers3), unlist(markers4))))
length(marker.set) 
## [1] 314

set.seed(1000110)
mnn.out2 <- fastMNN(pbmc3k, pbmc4k, subset.row=marker.set,
    BSPARAM=BiocSingular::RandomParam(deferred=TRUE))

最后作图看看

mnn.out2 <- runTSNE(mnn.out2, dimred="corrected")
gridExtra::grid.arrange(
    plotTSNE(mnn.out2[,mnn.out2$batch==1], colour_by=I(colLabels(pbmc3k))),
    plotTSNE(mnn.out2[,mnn.out2$batch==2], colour_by=I(colLabels(pbmc4k))),
    ncol=2
)

看样子根据这些marker基因进行的批次矫正，并没有去掉marker基因背后的生物学差异，因为两个批次数据还是能分出来群的，并且分群结果看似还不错

6 矫正后数据的应用

数据整合并矫正后，对合并后数据的cluster分析即可，不用再对每个batch的各个cluster单独分析。另一个好处是：batch整合后作为一整个数据集，细胞数量增加，相当于增加了群体结构的分辨率，方便下游的细胞层面分析（比如后面的轨迹推断）。

可能你还想利用合并后的数据进行基因层面的分析，例如差异分析marker基因鉴定。但一般不推荐，因为算法对多个批次进行合并时，不会保留每个基因表达的差异。如果要探索基因表达相关，最好还是要在未矫正数据基础上进行，并且还要把批次信息封锁掉（还记得之前findMarkers中的block参数吗？）

# 例如
m.out <- findMarkers(uncorrected, clusters.mnn, block=uncorrected$batch,
    direction="up", lfc=1, row.data=rowData(uncorrected)[,3,drop=FALSE])

TODO：对13.7 Using the corrected values继续补充