BatchBench比较scRNA批次矫正方法
当你的才华还撑不起你的野心时，请潜下心来，脚踏实地，跟着我们慢慢进步。不知不觉在单细胞转录组领域做知识分析也快两年了，通过文献速递这个栏目很幸运聚集了一些小伙伴携手共进，一起成长。
文献速递栏目通过简短介绍，扩充知识面，每天关注，希望你也能有所收获！
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的一致性表现是牺牲细胞类型数量得到的