010-BatchBench比较scRNA批次矫正方法

刘小泽写于2020.5.29 这篇文章做了一件事，就是帮助我们区分不同的批次矫正方法，然后比较了一下优劣

文章题目：Flexible comparison of batch correction methods for single-cell RNA-seq using BatchBench

文章在：https://www.biorxiv.org/content/10.1101/2020.05.22.111211v2

上传时间是：2020.05.22

BatchBench ：https://github.com/cellgeni/batchbench

BatchBench is a Nextflow workflow for running the following scRNA-Seq data batch effect correction methods:

• mnnCorrect
• limma
• ComBat
• Seurat 3
• Scanorama
• Harmony
• FastMNN
• BBKNN

前言

首先为什么要进行批次矫正？

单细胞分析经常会整合一些公共数据，不同的实验时间、文库制备、测序方案，都会产生一些技术误差，如果太多，可能会干扰真实的生物信号。因此来自这些非生物因素的干扰就称作批次效应

作者将8种常用的批次矫正方法分为3类：

• mnnCorrect、limma、ComBat、Seurat 3、Scanorama：产生一个整合、矫正后的表达矩阵
• Harmony、FastMNN：不是直接操作原始表达矩阵，而是对降维后的结果操作（they operate on a low-dimensional embedding of the original expression matrices），因此如果下游分析如果要用到原始表达矩阵的话，这类方法就会受限
• BBKNN：基于表达矩阵构建k-nearest neighbor graph（KNN），只能进行后续基于细胞的分析（如聚类、分群可视化），不能进行基于基因的分析（如marker基因鉴定、基因网络）

关于这8种方法：

以及这三类针对什么进行分析以及后续可以做什么，作者也作图说明：

但真正使用哪种方法，还是要依赖一个评测结果。但传统的评测只能针对已发表的方法，并且评测缺少一些高质量的数据集（比如尽可能多的包含批次效应的因素）

作者使用BatchBench，针对3种研究深入的数据集，对8种方法进行评测，这个方法的流程是：

三个数据集

Pancreas dataset

Baron (GSE84133)、Muraro (GSE85241)、Segerstolpe (E-MTAB-5061) 分别由inDrop, CEL-Seq2 和 Smart-Seq2产生。过滤细胞：细胞中基因表达量少于200；过滤基因：在少于3个细胞中表达。另外只保留有注释的细胞类型（去掉了unclassified这类的细胞）

Mouse Cell Atlas datasets

数据来自：https://figshare.com/s/865e694ad06d5857db4b

按照组织进行整合，得到了包含37个器官的数据集，其中选取了18个数据集（它们中包含大于1个批次并且有合理的细胞类型分布）。过滤细胞：基因表达量少于250；过滤基因：在少于50个细胞中表达；过滤细胞类型：细胞数量少于整体1%的类型；过滤批次：细胞数量少于总体5%的批次

Tabula Muris datasets

数据来自：https://www.google.com/url?q=https://figshare.com/projects/Tabula_Muris_Transcriptomic_characterization_of_20_organs_and_tissues_from_Mus_musculus_at_single_cell_resolution/27733&sa=D&ust=1589187433512000&usg=AFQjCNFC_0CGNwum-u2nka-OvFAmxoECtA

来自两个平台的同一组织的不同数据混合，得到11个器官的数据集。过滤细胞：基因表达量少于1000；过滤基因：在少于50个细胞中表达；过滤细胞类型：细胞数量少于整体1%的类型；过滤批次：细胞数量少于总体5%的批次。结果得到4168个基因，60828个细胞（40,058 from 10X and 20,770 from Smart-Seq2）

直接上结论

Seurat的整体效果最好，它既正确地整合了批次，又没有丢失不同细胞类型；

Harmony在pancreas和MCA的数据中表现也不错，但在矫正Tabula Muris数据时失败；Scanorama 和 fastMNN表现也算良好；

这里使用的熵评估方法，可能不太适用BBKNN，因此它需要额外的评测方法；

另外对于处理大量的细胞数量和批次，Harmony表现优秀，并且计算资源分配合理。除了Harmony和BBKNN，其他方法当遇到上百个批次的处理时（即使一个批次中的细胞数量不多）也会捉襟见肘，因此未来的批次效应处理方法应该向数据可扩展性（scalability）上发展。

如果想使用处理批次效应后的表达矩阵进行下游分析（如鉴定marker基因），这些方法都会遇到问题。因为marker基因并不是保守存在的，任何基于基因的分析（例如 找差异基因或者鉴定marker基因），都是基于基因表达量，而批次矫正方法需要保证不会干扰表达量的变化，这一点也是未来需要改进的。

结果

1 测试批次整合与细胞分群

使用了人类胰腺癌的3个scRNA数据集，原始数据的UMAP结果是：

左边是三个数据集，右边是各种细胞类型 但不得不说，两个图例使用的颜色太相近，容易引起混淆

可以看到，所有的方法都能将不同数据集的细胞混合起来，而依然可以分离不同的细胞类型

然后为了评估处理（先整合不同批次的细胞，然后分离不同的细胞类型）的效果，提出了计算一个”熵“：normalized Shannon entropy。如果批次方面的熵比较高，说明混合的批次之间更接近，也就是混合效果更好；如果细胞类型方面的熵比较低，说明细胞类型依然可区分

可以看到，不同的方法都保持较低的细胞类型方面的熵，因此它们都能够保证分离不同类型的细胞；但批次方面的熵差别较大。其中Seurat和Harmony整体表现较好，汽其次是Scanorama和fastMNN；而mnnCorrect, Limma 和 ComBat的表现较差

并且大部分方法对MCA（Mouse Cell Atlas）数据集的整合效果更好

图例：pancreas data (red), Mouse Cell Atlas (green), and Tabula Muris (blue)

2 当细胞数量或批次数量增加时，批次矫正变得困难

利用 Tabula Muris数据集（总共60,828 cells），取了它的1%、3%、5%、10%、20%、50%作比较

当细胞数量从608（1%）增至60828时，除了Scanorama、Harmony、Seurat，其他方法的批次熵都下降了50%左右。但是Scanorama在混合批次的同时，也混合了细胞类型（可以看到蓝色的虚线基本不变，说明细胞类型熵不变，也就是没有分离细胞类型）

Harmony是唯一一个在增加细胞数量后，批次熵增加的（图a）。除了Scanorama，其余方法的细胞类型熵都降低，说明细胞数量增多，细胞分群更容易

批次数量增加时，BBKNN, Seurat 和 Harmony表现最稳定（图d）

在时间方面，mnnCorrect和fastMNN随细胞数量增长，运行时间也呈现指数增长，mnnCorrect运行最慢。不过时间消耗在大部分软件中差别不大

在内存方面，所有的方法随细胞数量增长，内存消耗都呈现指数增长，其中Seurat消耗内存最多。综上，Seurat, mnnCorrect, ComBat 和 fastMNN是比较消耗资源的，而Harmony, Scanorama 和 BBKNN资源需求最小

a、d：熵的变化；b、e运行时间的变化；c、f：消耗内存的变化

3 批次处理对非监督聚类和marker基因鉴定的影响

使用三种非监督聚类方法：Leiden、Louvain、SC3，然后比较矫正前后的数据聚类结果。这个结果相似性的量化是利用Adjusted Rand Index (ARI)，图a可以看到：MCA数据集利用不同的方法聚类后结果相似，但组织之间的差异比较大；Tabula Muris数据集也是如此，不过与MCA不同的是，利用SC3方法得到的聚类结果会比Louvain或Leiden结果普遍（11个组织中有7个）有更高的ARI。奇怪的是，对于心脏和乳腺组织，最佳的聚类结果发生在：SC3+非批次矫正的数据。对于pancreas数据，SC3倾向于得到更大的ARI，而且不想MCA数据，Seurat和Harmony的聚类结果与之前的熵分析结果也一致。

对整合后的数据进行marker基因鉴定，只有ComBat和Limma的结果可以找到大部分细胞类型的marker基因，Seurat只能对少部分细胞类型进行鉴定（图b），但如果检测单个数据的marker基因与混合矫正后的marker基因之间的一致性，Seurat的一致性更强（图c）。Seurat的一致性表现是牺牲细胞类型数量得到的