2.3.1 算法|2020-BatchBench比较scRNA批次矫正方法
刘小泽写于2020.5.29
这篇文章做了一件事,就是帮助我们区分不同的批次矫正方法,然后比较了一下优劣
文章题目:Flexible comparison of batch correction methods for single-cell RNA-seq using BatchBench
文章在:Flexible comparison of batch correction methods for single-cell RNA-seq using BatchBench
上传时间是: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
来自两个平台的同一组织的不同数据混合,得到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的一致性表现是牺牲细胞类型数量得到的
最后更新于