scRNA-seq 的原始数据格式和目前大多数 scRNA-seq 分析过程都基于 FASTQ 文件(或压缩格式 fq.gz)。Illumina 平台测序数据默认生成 BCL 格式文件,可以通过 CellRanger mkfastq 进行转换。scRNAseq 的分析流程包括数据预处理、处理和扩展下游分析(图 4),其中数据预处理包括质控、read 比对和表达量化;数据处理包括标准化、批次效应校正、归一化、特征选择(HVG 选择)、降维与聚类、细胞分型注释、差异表达分析(DEGs)、可视化;扩展下游分析包括拟时序、细胞间相互作用(CCI)、通路富集分析、基因调控网络(GRN)等下游分析。整体来看,scRNAseq 分析方法层出不穷,没有绝对完美适用于所有场景的方法,分析工具重要的是获取生物学信息,难点在于选择最合适的方法。本文中,QY球友会将提出总结常见的单细胞转录组分析方法,并对其优缺点和适用范围提出建议。
图 4:单细胞分析概览。A. 在预处理阶段,基于测序数据,细胞 - 基因矩阵读数通过单细胞读数校正和定量产生;B. 分析使用的高质量细胞矩阵通过原始的基因表达矩阵获得,通过去批次效应矫正批次,通过标准化降低生物学差异,补充未检测到的基因;C. 依照或不依照先前的参考信息对细胞类型进行注释;D. 转录组特征相似的细胞被归为一类,称为“细胞簇(cluster)”,细胞的可视化通过降维方法实现,差异基因分析对组间差异进行检验;E. 拟时序分析重建细胞转录水平变化的动力学过程;F. 细胞间转录组调控关系可以通过胞间互作分析进行推断。
数据预处理
将原始测序数据通过滤除低质量 reads 和环境干扰与参考基因组进行比对和量化。从而得到每个细胞的特征计数矩阵和记录其他信息的辅助文件,用于下游的数据分析(图 4 A)。
1、质控
由于测序仪器问题、人为操作、细胞自发情况,或存在空液滴、双细胞、死细胞等,不可避免地会产生低质量的测序数据(Chen 等,2019a ; Hao 等,2021b)。空液滴通常出现在液滴捕获细胞外背景转录本而不是细胞时(Ilicic 等,2016 ; Kolodziejczyk 等,2015)。一种高度主观的方法是根据曲线的膝点确定一个 UMI 阈值,并过滤掉 UMI 计数低的细胞。随后使用 DropEst (Petukhov 等,2018)、EmptyDrops (Lun 等,2019) 和 DIEM (Alvarez 等,2020) 增强过滤效果。DropletQC (Muskovic and Powell, 2021) 量化未剪接前 mRNA 含量的核分数得分。MT 基因阈值虽然是衡量死细胞的标准,但它的选择需要综合考虑细胞生理因素 (Subramanian 等,2022)。近年来,基于深度学习的方法也应运而生,例如基于神经网络的 EmptyNN (Yan 等,2021) 和基于深度生成模型的 CellBender (Fleming 等,2019),能够有效识别空液滴中的背景转录本。
双细胞是指两个细胞包含在一个液滴中的情况,根据转录分布可分为同源双峰和异源双峰,均服从泊松统计量(Bloom, 2018)。绝大多数方法基于基因表达计算,利用先验知识或深度学习获取单峰与双峰细胞的差异,然后训练分类器进行筛选,例如基于最近邻的 DoubletFinder (McGinnis 等,2019a)、Scrublet (Wolock 等,2019);基于反卷积的 DoubletDecon (DePasquale 等,2019)、基于变分自编码器的 Solo (Bernstein 等,2020) 和基于集成算法的 Chord (Xiong 等,2021a)。此外,Scds 是另一种筛选方法,它依赖于基于共表达的双联体打分和基于二分类的双联体打分策略,实现 scRNA-seq 表达数据的双联体分离 (Bais and Kostka,2020)。一些方法使用其他特征,例如 demuxlet 它使用自然遗传变异信息指导实验并通过计算进行过滤 (Kang 等,2018)。
合理的质控需要综合考虑技术性和生物性因素,这也是当前研究的主要方向。最近一种由生物数据驱动的自学习无监督质控方法 ddqc 被提出来,用于确定各种 GC 指标的具体阈值 (Macnair and Robinson,2023)。
2、reads 比对和定量
质控后剩余的高质量细胞需要将这些短 reads 映射到特定的参考基因组上进行比对,以此对基因表达水平进行定量。RNA 比对通常分为两步:比对 reads 以建立索引和映射 RNA 剪接序列,前一步与 DNA reads 比对共用,解决错配问题并设置索引参考;后一步是 RNA reads 比对所特有的,提供连通性信息。
早期二代测序结果是几十对长度的碱基 reads。Seed-to-extend (Buhler,2001)(包括 MAQ (Li 等,2008a)、SOAP (Li 等,2008b)、CloudBurst (Schatz,2009)、ZOOM (Lin 等,2008))、BurrowsWheeler 变换方法 (Burrows and Wheeler,1994)(包括 SOAP2 (Li 等,2009)、Bowtie (Langmead 等,2009)、BWA (Li and Durbin,2009))、Needleman-Wunsch 方法(包括 Novocraft (Hercus,2009))和 suffix-tree 算法方法(包括 MUMmer 2 (Delcher 等,2002))都是百万级短链 DNA 测序 reads 比对的有效工具。Bowtie 采用了一种依赖于 Burrows-Wheeler Transforming 的 FM-index 方法,如果 reads 有多个准确匹配则结果只报告一个,与 MAQ(Ferragina and Manzini,2001)相比,大大优化了运行内存和比对速度。BWA 是另一种基于 BWT 的比对方法,使用新的 SAM(Sequence Alignment/Map)格式输出比对结果。基于 MAQ 和 Bowtie 两种短链 DNA 比对算法,Cole Trapnell 于 2009 年提出了第一个针对 NGS 数据的 RNA-seq 比对方法 TopHat,它使用 2 -bit-per-base 编码实现 reads 与哺乳动物基因组中剪接位点的有效比对,而无需事先知道剪接位点的具体信息(Trapnell 等,2009)。
上述方法在碱基对长度超过 50 bp 时比对精度急剧下降(Gupta 等,2018 ; Lebrigand 等,2020)。NGS 单细胞测序分析主要采用两类方法:基于 Bowtie2 的方法和基于 seed 策略的方法(Langmead and Salzberg, 2012)。Bowtie2 是 Bowtie 的升级版,保留了 FM-index 依赖的 BWT 算法核心,允许有间隙比对,并使用单指令多数据(SIMD)扩展到长测序比对,同时提高运行速度。Daehwan Kim 在 Bowtie2 基础上,先后提出了 TopHat2(Kim 等,2013)和 HISAT(Kim 等,2015)。种子策略主要有 STAR(Dobin 等,2013)和 Subread(Liao 等,2013)。STAR 基于最大可映射前缀(MMP)的思想,采用顺序检索的策略,将与参考匹配的最长部分 reads 设为种子 1,其余 read 继续匹配,依次从种子 2 调用至种子 n。值得注意的是,Rsubread 完全基于 R 语言平台实现了第一次 read 比对和基因量化的过程(Liao 等,2019)。
基因表达量化又可分为伪比对量化和基于 read 比对的量化。伪比对是指不采用上述严格的两步法将所有 reads 比对到参考基因组上,包括选定的 k-mers 比对方法(Sailfish(Patro 等,2014)、Kallisto(Bray 等,2016)、Salmon(Patro 等,2017)、RapMap(Srivastava 等,2016)和 Barcode-UMI-Set (BUS) 比对方法 BUStools(Melsted 等,2019)。Kallisto-BUStools 是最新的工作流程,它使用 BUS 文件格式进行初始数据预处理,与 BUStools 一样,伪比对结果和量化计数都保存在 BUS 文件中(Melsted 等,2021)。另一方面,基于 reads 比对的方法依赖于 RNA reads 比对方法的结果来量化基因。CellRanger 是 10x Genomic 公司指定替代 Longranger 的官方开源数据预处理软件(Zheng 等,2017)。STARsolo 是替代 Cellranger 的 mapping/quantification 功能的工具,可实现多平台测序数据的分析和基因表达之外的转录组特征的量化(Kaminow 等,2021)。其他基于 reads 比对的基因表达定量方法如 UMItools (Smith 等,2017)、zUMIs (Parekh 等,2018)、Alevin-fry (He 等,2022)、DropEst (Petukhov 等,2018)、RainDrop (Niebler 等,2020)、baredSC (Lopez-Delisle and Delisle, 2022)、BCseq (Chen and Zheng, 2018) 使用各种质量过滤器和 barcode/UMI 处理策略在一定程度上提高了 CellRanger 的性能。
CellRanger 和 STARsolo 在处理包括 10x Chromium 在内的各种单细胞转录组数据集时都具有良好的运行速度,并且准确度极高。但在获得几乎相同结果的前提下,后者相比前者提升了至少 5 倍的运行速度,这也验证了 Alexander Dobin 等人使用 STARsolo 取代 CellRanger 的目的(Brüning 等,2022 ; Chen 等,2021a ; You 等,2021)。
数据处理
在对表达矩阵进行必要的调整(Normalization、Batch Effect Correction、Imputation)后,即可从单细胞转录组数据中充分挖掘出生物信息进行分析。Seurat 和 Scanpy 分别基于 R 和 Python 对上述过程进行模块化、可扩展的处理,是目前单细胞转录组数据的主流分析流程(Satija 等,2015 ; Wolf 等,2018)。常规分析流程和预期处理结果可参见总分析框架(图 4 B-D)。
1、标准化
在测序过程中,技术原因或者细胞本身的生物学差异可能造成同一样本内(细胞之间)或者不同样本之间的文库大小差异(Marinov 等,2014)。无限数方法按照文库大小进行处理,按照具体原理大致可以分为基于全局缩放的标准化、spike-in 标准化和其他数据变换模型标准化。
全局缩放方法最初是为 bulk RNA 分析而发展起来的,通过特定的缩放因子对全局数据进行缩放(Finak 等,2015)。每万计数(CPT)变换和每百万计数(CPM)变换是常见的线性缩放方法,在不考虑 spike-in count 的情况下,它们都对每个 UMI/ 总 UMI count 等距缩放。其他标准化方法包括每百万 reads 数(RPM)(Mortazavi 等,2008)、修剪均值 M 值(TMM)、DESeq(Robinson and Oshlack, 2010)、上四分位缩放(Bullard 等,2010)、FPKM(Trapnell 等,2010)、RPKM(Tu 等,2012)等,它们对于极值的稳定性比线性缩放更好,因此与 RPKM/FPKM 一样具有更广泛的应用范围。但单独使用该类方法进行单细胞转录组的标准化时,由于数据的稀疏性和假阳性率虚高,效果并不可接受(Evans 等,2018),与特定方法结合时往往需要改进。SCnorm 使用分位数回归方法来评估不同测序深度依赖细胞组之间的尺度因子(Bacher 等,2017)。bayNorm 基于基因原始计数与真实计数服从负二项(NB)分布的假设,使用集成贝叶斯模型对 scRNA-seq 数据进行标准化(Tang 等,2020)。
spike-in 标准化方法可以看作是全局尺度方法的另一种扩展,因为尺度因子是根据 spike-in 基因计算出来的。需要注意的是,将 RNA spike-ins 的信息添加到其他方法中也可以提高 SCnorm 等标准化的效果。GRM 是一种基于 spike-in ERCC 分子浓度伽马分布的方法,其中 ERCC 是测序中常用的校准材料(Ding 等,2015)。BASiCS 是一种自动贝叶斯标准化方法,将泊松分层模型应用于 spike-in(技术)基因,以推断细胞特定的标准化常数(Vallejos 等,2015)。
以上方法都是在细胞内 RNA 数量恒定的假设下对基因进行缩放,而这可能具有欺骗性,因此其他转化模型采用了不同的策略。由于单细胞转录组数据存在零膨胀问题,一些模型就是为此而设计的,例如相对对数表达(RLE)方法 ascend(Senabouth 等,2019)和基于 NB 的模型,如 Dino(Brown 等,2021)、scTransform(Hafemeister and Satija, 2019)。其他转化模型归一化方法如 MUREN 使用最小二乘(LTS)回归算法(Feng and Li, 2021);Sanity 使用从 UMI 计数推断出的对数转录商(LTQ)作为贝叶斯框架的输入,以避免泊松波动,因为 LTQ 向量的变化估计了基因表达值(Breda 等,2021);PsiNorm 是一种基于无监督帕累托分布尺度参数的方法,用于提升标准化效率和准确率(Borella 等,2021)。Charles Wang 比较了 sctransform、TMM、DESeq 等共 8 种标准化方法,其中 sctransform 和 logCPM(Seurat 的内置处理方法)受数据影响最小,在可变数据集上最稳定(Chen 等,2021a)。
2、批次效应校正
由于实验设计、测序平台、测序时间、人员操作流程等原因,不同的单细胞转录组测序数据在 mRNA 捕获效率、测序深度等会存在明显差异,从而产生样本间的批次效应(Chen 等,2019a;Hwang 等,2018;Tung 等,2017)。理论上可以通过实验策略消除技术变异,但由于实验过程的客观限制以及测序仪器误差,不可避免地会或多或少地引入批次效应。利用计算方法进行校正是解决不完善实验设计的必要手段,通常使用的方法可以分为相互最近邻(MNN)方法、基于潜在空间的方法、基于图的方法、DL 方法和其他方法。
MNN 首先识别出不同批次之间同一细胞类型的最相似细胞,然后利用这些细胞进行批次效应校正,包括 batchelor(Haghverdi 等,2018)、Scanorama(Hie 等,2019)、Canek(Loza 等,2022)。另一类使用 MNN 的方法是基于降维后的潜在空间,如 Seurat (Satija 等,2015)、BEER (Zhang 等,2019b)、SMNN (Yang 等,2021a)、iSMNN (Yang 等,2021b)。例如,Seurat 使用典型相关分析 (CCA) 潜在空间中的 MNN 对 (称为“锚点”)来匹配相似细胞,而 BEER 使用主成分分析 (PCA) 子空间来筛选相似性较差的子群。SMNN 和 iSMNN 分别采用监督机器学习和迭代监督机器学习来细化从预校正细胞聚类或迭代细胞聚类信息中训练出的 MN 对。
基于潜在空间的方法是指在隐藏空间或降维后的嵌入中进行批次效应校正的方法,除了基于 MNN 聚类的策略外,还有与 PCA 相关的空间方法 Harmony(Korsunsky 等,2019)、FIRM(Ming 等,2022)、Monet(Wagner, 2020);t 分布随机邻域嵌入 (tSNE) 空间方法 sc_tSNE(Aliverti 等,2020)和 ZINBWaVE(Gao 等,2019)。Harmony 被广泛用于去除样本间的批次效应,使用 PCA 方法将排序的细胞输入到单个公共嵌入中,然后在最大多样性聚类和线性批次校正之间迭代循环,直到为每个细胞分配一个特定的校正因子,可用于后续的批次效应去除。Sc_tSNE 方法引入梯度下降算法对传统 t -SNE 算法进行优化,随后采用线性校正(Aliverti 等,2021)。ZINB-WaVE 最初设计用于在单细胞数据中进行基因提取,Risso et al.(2018)将该方法扩展至小批量优化。
基于图的方法利用细胞基因表达矩阵将数字信息转化为空间构造的图,其中节点代表不同类型的批次,边的权重基于不同的计算方法。BBKNN 利用 k 近邻细胞构建图(KNN 图),通过使用均匀流形近似与投影(UMAP)方法合并不同数据集间单个细胞的图实现批次效应校正,这也是 Scanpy 工作流程中的默认方法(Pola ński 等,2020 ; Wolf 等,2018)。王波在 OCAT 中提出“幽灵细胞”(默认为 k-means 算法聚类中心)来制作细胞连接的二分图(Wang 等,2022a)。
近年来,深度学习方法的快速发展也为批次效应校正提供了新思路,实现高效、大通量的数据处理,如 INSCT(Simon 等,2021)(三重态神经网络)、CLEAR(Han 等,2022)(自监督对比学习)、BERMUDA(Wang 等,2019e)(迁移学习)、iMAP(Wang 等,2021a)(VAE-GAN)、ResPAN(Wang 等,2022e)(Wasserstein GAN),一些新方法被证明在批次效应校正方面有更好的效果;例如,基于从 SciBet 学习到的带注释数据集中的生物学先验知识,SSBER 可以在大型 RNA 测序数据集中去除批次效应(Zhang and Wang,2021)。建议在整合单细胞转录组数据之前,应根据数据的实际情况先测试多种方法,然后选择最合适的批次效应去除方法。例如,Jinmiao Chen 团队和 Charles Wang 团队分别于 2020 年和 2021 年对本综述 2.2 中提到的前三种方法的大部分进行了基准测试,证明了 Harmony 和 Seurat V3 在大多数情况下都能达到良好的批次效应校正效果,这符合这两种方法如今仍然被广泛使用,但对于深度学习方法来说仍然缺乏好的指标这一事实(Chen 等,2021a;Tran 等,2020)。
3、填补
测序过程中会引入大量 0 值(高通量大规模 10x 基因组测序数据中零值可能超过 90%)(Stegle 等,2015 ; Talwar 等,2018),这会干扰下游生物学差异分析,因此必须对原始基因表达矩阵中的缺失数据值进行填补,同时有效区分技术噪音零值与生物学零值。
基因 / 细胞分离方法主要应用于早期的填补,其分别考虑细胞相似性(MAGIC (van Dijk 等,2018)、Sclmpute (Li and Li, 2018)、VIPER (Chen and Zhou, 2018)、RESCUE (Tracy 等,2019)、scRMD (Chen 等,2020a)、scRoc (Ran 等,2020))或基因间关系(SAVER (Huang 等,2018a)、SAVER-X (Wang 等,2019a)、G253 (Wu 等,2021e)、DCA (Eraslan 等,2019)、DeepImpute (Arisdakessian 等,2019))。总体而言,这些方法缺乏对数据集整体的考虑,容易导致过度插补或者引入误差(Zhang 等,2019d)。综合方法综合考虑细胞与基因之间的联系:CMF-Impute 和 netNMF-sc 是最早有效利用细胞与基因之间的关联进行插补的方法(Elyanow 等,2020;Xu 等,2020a)。scIGANs 通过特定的 GAN 模型处理基因表达矩阵,利用生成的细胞训练 GANs 模型来插补 dropout(Xu 等,2020b)。近年来,新的方法还在不断被提出,以更好地解决 dropout 之外的技术噪声对数据的影响,并实现对生物零值的更好的区分。AutoClass(Li 等,2022c)实现了无监督处理,而 ALRA 方法主要针对生物零值(Linderman 等,2022)。scMOO 进行了根本性的改变,利用数据的潜在结构来学习细胞相似性垂直结构和总低秩结构中的深度关联,从而取得了比单一基因表达矩阵作为输入更好的插值效果,但对内存的要求也更高(Jin 等,2022a)。sc-PHENIX 利用 PCA-UMAP 初始化方法,实现了基因表达的非线性插值(Padron-Manrique 等,2022),目前哪种插值能取得最佳效果尚无明确定论。由于数据集本身的原因,下游分析的目的会有不同的选择,但毫无疑问最好的填补方法将能够以较低的计算要求有效区分技术噪声零值和生物零值(Jiang 等,2022a;Wen 等,2022)。
4、特征选择
为了降低数据维数以提升计算分析效率、减少技术噪声干扰和模型过拟合的风险,QY球友会常常采取特征选择策略,选取不同细胞中差异较大的基因,而非整个数据集基因进行聚类等后续分析(Brennecke 等,2013;Jackson and Vogel,2022;Svensson 等,2017)。
在 bulk RNA-seq 分析中,寻找差异基因的方法通常包括基于倍数变化(FC)的方法、基于统计检验的方法和 FC- 统计检验方法,显然后者的筛选结果和可信度最好(Chung and Storey,2015)。
早期的单细胞特征选择方法缺乏平均表达量与方差之间的校正,导致结果中高表达基因的比例过高(Brennecke 等,2013)。EDGE 采用大量弱学习器的集成学习方法来学习细胞间相似性概率,提取基于信息熵的显著贡献作为高可变基因(Sun 等,2020c)。同样,SAIC 基于迭代聚类最终输出实现了最优细胞簇分离(Yang 等,2017)。近期,一些新的特征提取策略被提出并证明了其稳定性和有效性,但它们之间的性能权威验证尚缺乏:包括基于基因表达分布矩阵的方法 SCMER(Liang 等,2021b)、RgCop(Lall 等,2021)、scPNMF(Song 等,2021a)、SIEVE(Zhang 等,2021e);基于熵的方法 IEntropy(Li 等,2022g)、infohet(Casey 等,2023);综合考虑聚类的方法有 Triku(Ascensión 等,2022)、FEAST(Su 等,2021)等。由于上述方法绝大多数忽略了整体的依赖于基因表达的特征,因此提出了综合的方法,如 Triku 使用 k 最近邻图的方法对基因表达模式进行综合探索和分类,实现无偏差地筛选出更有生物学意义的特征基因;FEAST 在共识聚类上通过 f 检验对特征进行排序,并基于特征评估算法提取 HVG(Wang 等,2022c)。
其他一些方法使用高可变基因以外的特征来表示数据集,例如 scVEGs 和 scSensitiveGeneDefine 方法,使用高变异系数(CV)作为特征提取;BASiCS 方法利用 spike-in 基因的信息(Chen 等,2016b;Chen 等,2021b)。总体来看,基于准确性、生物学可解释性等角度,当前特征选择的主要目标是有效提取 HVG,以便对高维转录组数据进行有效的下游分析。
5、降维
由于单细胞转录组通常包含数万个甚至更多的基因,不利于直接提取有效信息,在实际分析过程中,通常需要对原始测序数据进行降维。除了利用前文提到的特征选择方法处理高维单细胞转录组测序数据外,降维也是一种有效的方法,根据降维策略可分为线性降维(基于潜在狄利克雷分配(LDA)的方法、基于 PCA 的方法)和非线性降维(基于 t -SNE 的方法、基于 UMAP 的方法)(Andrews and Hemberg,2018;Becht 等,2019;Laurens and Hinton,2008;Peres-Neto 等,2005)。
在线性降维中,LDA 和 PCA 是两种广泛使用的算法,LDA 从分类最大的角度区分特征,而 PCA 则从方差最大的角度正交提取主成分。尽管有 JPCDA、LDA-PLS 等改进算法,但是 LDA 模型在单细胞转录组数据中的降维效果仍然不是最优的(Tang 等,2014 ; Zhao 等,2020)。PCA 是另一种线性变换,Seurat 通常根据标准差 -PC 图的拐点或者 PC 的比例检验结果 P 值(ScoreJackStraw 函数)来确定 PC 数量的多少。其他基于 PCA 的降维方法的变体包括 pcaReduce(Žurauskien ė和 Yau,2016),GLM-PCA(Townes 等,2019),RPCA(Gogolewski 等,2019),tRPCA(Candès 等,2011),scPCA(Boileau 等,2020),PCAone(Li 等,2022l)。GLM-PCA 将传统 PCA 分析扩展到非正态分布,通过引入指数家族似然策略直接处理原始矩阵,使 PCA 摆脱正态化限制,然后使用偏差对基因实现进行排序和提取(Collins 等,2002)。ScPCA 使用对比 PCA 和稀疏 PCA 分别去除技术噪音和数据,进一步增加了 PCA 的稳定性(Abid 等,2018 ; Zou 等,2006)。由于大多数 scRNA-seq 数据集难以通过简单的线性降维进行有效表示,解决这一问题的第一个策略是基于快速 PCA 分析方法。PCAone 提出了一种新的快速随机奇异值分解(RSVD)策略,在 35 分钟内完成 130 万小鼠脑细胞单细胞数据的分析(Li 等,2022l)。
非线性降维是另一种解决方案,如非参数降维方法 t -SNE 和 UMAP,都需要预先设置聚类的超参数;而在分类效果上,前者倾向于离散数据中细胞的形成。在合理使用参数设定的情况下,UMAP 与 t -SNE 并无明显差异,即在使用相同的信息初始化方法后,二者可以在保留数据集全局结构的同时,产生近似的分析效率(Do and Canzar,2021;Kobak and Linderman,2021)。针对 t -SNE 的改进方法包括 net-SNE、qSNE、FItSNE、联合 t -SNE(Cho 等,2018a;Linderman 等,2019;Wang 等,2022b),而 UMAP 的改进主要来自于 Leland McInnes 课题组对该方法的自我改进(McInnes 等,2018)。为了更好地可视化 t -SNE 或 UMAP 的降维结果,Hyunghoon Cho 提出了基于局部半径依赖优化的转录组变异信息 den-SNE/densMAP 方法,以迭代优化传统 t -SNE/UMAP 的功能;Stefan Canzar 提出了 j -SNE/jUMAP 来改善多模态组学数据联合可视化结果,减少可视化的误导性(Do and Canzar,2021;Narayan 等,2021)。
6、聚类
在单细胞转录组数据分析中,通过聚类将细胞划分为亚群,从而能够表征多细胞生物中不同细胞类型,这有助于QY球友会从细胞异质性的角度准确地分析不同的组织或发育过程。聚类的实际效果会受到数据预处理步骤的影响,例如浴效应归一化、归纳、降维等。
在特征基因选择和降维之后,绝大多数单细胞是基于距离进行聚类的。K 均值聚类算法的概念被用于 SCUBA、SC3 和 RaceID 等应用(Grün 等,2015;Kiselev 等,2017;Macqueen, 1967;Marco 等,2014)。在参数选择改进方面,SAIC 通过 Davies-Bouldin 指数迭代优化多个初始中心 K 和 P 值,以获得最优解;LAK 将参数选择算法应用于数据集,实现参数的自动选择(Davies and Bouldin, 1979;Hua 等,2020;Yang 等,2017)。在超高维数据的操作中,LAK 添加 Lasso 惩罚项进行标准化,mbkmeans 使用小批量 k 均值实现百万细胞级别的快速聚类(Hicks 等,2021)。SMSC 应用谱聚类方法来提高聚类性能,但对于超高维数据会损失一定的准确性(Qi 等,2021)。另一大类广泛使用的距离聚类方法依赖于共享最近邻图结构和图聚类,其中使用最广泛的是 Louvain 或 Leiden(Blondel 等,2008;Xu and Su, 2015)。稀有细胞的识别需要结合特定方法进行改进,例如 dropClust 使用局部敏感哈希工作流筛选最近邻,然后是 Louvain 聚类,它使用指数衰减函数来保留更多稀有细胞的转录组特征(Sinha 等,2018)。其他基于距离的聚类方法使用不同的算法核心:SIMLR 使用高斯核学习模型为数据集中潜在的 C 细胞群体构建核矩阵,而 Conos 提出联合相互最近邻(mNN)图聚类来实现对多个不同单细胞转录组样本的整合分析(Barkas 等,2019 ; Wang 等,2017a)。基于密度的聚类利用样本分布的接近程度进行聚类,DBSCAN 是最经典的算法(Ester 等,1996 ; Fukunaga and Hostetler, 1975)。对于单细胞测序,densityCut 和 FlowGrid 就是基于此原理设计的(Ding 等,2016 ; Fang and Ho, 2021)。层次聚类是一种自下而上的聚类方法,通过无监督学习,不断重复计算细胞与细胞的相似性进行分类,直至完成预设的聚类数(Guo 等,2015)。随后,RCA 聚类采用常规的层次聚类方法,对映射到全局参考面板的细胞进行聚类;HGC 在 SNN 图上构建层次树(Li 等,2017;Zou 等,2021)。为了解决常规层次聚类方法难以对某一组细胞进行聚类、只允许同一组特征基因进行聚类的缺陷,K2Taxonomer 采用约束 K 均值算法扩展到样本组,基于多个基因集递归进行积分计算,以捕获各种分辨率下的亚组(“类似分类学的细胞”)(Reed and Monti, 2021)。Mrtree 将层次聚类的策略应用于平面簇的多个划分,并构造多分辨率协调树用于细胞聚类(Peng 等,2021a)。最近,Zelig 和 Kaplan(2020)提出了一种 KMD 聚类方法,通过平均链接层次聚类模型在聚类时消除了超参数 K,大大减少了主观性带来的判断误差。
深度学习聚类方法是将机器学习方法与上述单细胞转录组聚类策略相结合,可以以无监督、监督或半监督的形式实现更高效的聚类结果。这些方法倾向于学习一种非线性变换,通过将原始高维数据映射到较小的潜在空间中来获得最佳的低维表示。总体而言,这种方法避免了传统聚类方法对聚类前数据处理方法选择的影响。无监督聚类方法包括 ADClust、DESC、SAUCIE、VAE-SNE 等,通常不需要预设聚类个数等参数,以自主学习的方式完成对数据集的分析处理(Amodio 等,2019;Graving and Couzin,2020;Li 等,2020c;Zeng 等,2022c)。虽然无监督聚类方法避免了手动输入聚类个数等参数,可以延伸到超高维细胞聚类,但有时利用高质量标注数据集或其他先验知识辅助约束进行监督或半监督聚类,可以实现更为准确的细胞类型分类,提高聚类性能(Bai 等,2021)。基于迁移学习的 ItClust、基于互监督 ZINB 自编码器和图神经网络(GNN)的 scDSC、基于软 K 均值卷积自编码器的 ScCAE、基于 Cramer-World 距离最大均值惩罚高斯混合自编码器的 SeGMA、基于时间序列聚类网络 STCN 都是广泛使用的监督聚类(Gan 等,2022 ; Hu 等,2022a ; Hu 等,2020a ; Ma 等,2021b ; Smieja 等,2021)。此外,Zhang 团队(Yang 等,2023b)利用分层 GAN 设计了另一种广泛使用的深度学习方法 IMDGC,用于单细胞转录组数据分析,以生成的方式构建细胞嵌入簇。
针对聚类中的特殊情况,设计了有针对性的聚类方法:GiniClust(Jiang 等,2016)(更新为 GiniClust 3(Dong and Yuan,2020)、MicroCellClust(Gerniers 等,2021)用于稀有细胞亚群聚类;EDClust(Wei 等,2022)、ENCORE(Song 等,2021b)和 MLG(Lu 等,2021)用于降噪和消除批次效应;ClonoCluster(克隆起源信息)(Richman 等,2023)、IsoCell(可变剪接信息)(Liu 等,2023)使用附加信息进行聚类。Wu 和 Yang 从特征选择对聚类的影响的角度对聚类方法进行了评估,他们证明更具代表性的特征选择会提高细胞聚类的水平,基于“聚类相似性”的方法(QY球友会综述中提到的大多数基于距离的聚类方法)通常具有广泛的高聚类类型性能;然而,高精度和高运行速度需要根据实际数据集进行有针对性的选择(Su 等,2021;Yu 等,2022)。双重浸入 (double dipping) 是一个显著的问题,即在细胞聚类和差异表达基因中使用相同的表达数据,导致在细胞聚类不正确时 DE 基因的错误发现率 (FDR) 过高。例如,如果只存在一个特定的细胞聚类,则不应将任何基因视为差异基因。为了系统地解决这个问题,ClusterDE 采用了聚类对比学习策略进行聚类后 DE 测试。该方法与截断正态分布 (TN) 检验和 Countsplit 方法相比,在不同阈值范围内具有更好的 FDR 控制 (Song 等,2023a)。
7、细胞类型注释
细胞分型注释是指利用特定的信息对单细胞测序数据集中的细胞或细胞亚群进行注释,作为后续生物学分析的基础。最常用的策略是对细胞进行无监督聚类,然后根据标记基因进行注释,例如 scCATCH、SCSA (Cao 等,2020b;Shao 等,2020a),但它难以处理复杂的高维数据集 (Franzén 等,2019;Luecken and Theis, 2019;Zhang 等,2019c)。目前已经开发了多种自动细胞分型方法,大致可分为两类,即依赖参考和无参考的注释方法。
依赖参考信息的注释方法要求用户提供预先注释的高质量单细胞转录组数据集或来自 PanglaoDB 数据库、ScType 数据库等的先验知识进行比对(Ianevski 等,2022)。根据方法原理的不同,可分为基于层次树的方法(CHETAH(de Kanter 等,2019)、Garnett(Pliner 等,2019)、HieRFIT(Kaymaz 等,2021)、scHPL(Michielsen 等,2021)、scMRMA(Li 等,2022e)、基于相似性的方法(SingleR(Aran 等,2019)、scmap(Kiselev 等,2018)、deCS(Pei 等,2023)、scID(Boufea 等,2020)、scMatch(Hou 等,2019)、Symphony(Kang 等,2021)、基于签名基因的方法(Cellassign(Zhang 等,2021)、基于特征基因的方法(Cellassign(Zhang 等,2022)。al., 2019a )、Cell-ID(Cortal 等,2021)、scMAGIC(Zhang 等,2022g)、SciBet(Li 等,2020b)和其他 DL 方法。作为早期方法,ACTINN 是一种使用 3 个隐藏层神经网络进行注释分类的深度学习方法(Ma and Pellegrini, 2020)。SCPred 随后提出了一种基于嵌入的无偏特征选择的机器学习概率预测方法(Alquicira-Hernandez 等,2019)。其他方法如 Seurat 在 PCA 空间中投影查询细胞并通过加权投票分类器训练细胞分型注释;scSorter 采用高斯混合模型,GraphCS 使用虚拟对抗训练 (VAT) 损失修改的 GNN 来扩展到多物种、大规模细胞注释数据集(Guo and Li,2021;Zeng 等,2022a)。
不依赖参考信息的注释方法使用预先训练的深度学习模型,可以直接使用查询数据集作为输入进行细胞分类。scDeepSort 使用来自人类细胞图谱 (HCL) 和小鼠细胞图谱 (MCA) 数据库的单细胞图谱作为预训练加权 GNN 模型的输入,该模型适用于人类和小鼠细胞注释并取得良好的效果(Han 等,2018b;Han 等,2020;Shao 等,2021b)。类似地,Pollock 是一个预训练的人类癌症参考 VAE 模型,用于对癌症环境中的多模态细胞进行分类(Storrs 等,2022)。虽然使用起来更方便,但对于差异显著的查询数据集难以达到更好的细胞注释效果,而且由于准确性和预训练参考数据集的数量也难以扩展应用。还有一些其他用于有针对性领域研究的细胞注释工具,例如,用于人类肾细胞注释的 DevKidCC(Wilson 等,2022),用于识别癌症和正常细胞的 ikarus(Dohmen 等,2021)。总体而言,无参考注释方法的性能受到预训练参考数据集的覆盖率和准确性的制约。
目前,改进细胞注释工具以在大平台和多细胞模式下统一分配细胞类型是细胞注释研究的主流方向,最新的 Cellar 和 ELeFHAnt 方法在这方面做了一些尝试并取得了初步成果(Hasanaj 等,2022 ; Thorner 等,2021)。总体而言,基于相似性的注释方法计算量大,在应用于非常大的查询和参考数据集时,往往会在准确率和速度之间做出权衡,因此一般只适合在较小的数据集中进行细胞分类;对于较大规模的数据集,建议使用 F 检验特征选择或 MLP 分类器(Hu 等,2020a ; Huang and Zhang, 2021 ; Ma 等,2021c)。此外,半监督迁移学习的方法,如 Itclust,在发现新的细胞亚型方面也有不错的效果。近年来,基于上述参考注释方法分类的新方法不断完善,VAE 等深度学习模型也在该领域得到应用。
8、差异表达分析(DEG)
统计检验是 Bulk RNA-seq 的差异基因分析中常用到的方法,类似章节 2.4HVG Selection 算法:通常以 P 值和对数倍变化量作为重要参数。统计检验包括 t 检验(两个样本为基础),Wilcoxon 检验,Kolmogorov-Smirnov 检验(KS 检验),Kruskal-Wallis 检验(KW 检验),其中一些在单细胞转录组 DEGs 的检验中也被广泛使用。基于此,发展了相应的检测工具:limma(Ritchie 等,2015),edgeR(Robinson 等,2010),DESeq2(Love 等,2014)。limma 和 edgeR 算法均由 Smyth GK 提出,前者基于正态或近似正态分布模型,后者基于过度离散的泊松分布模型。DESeq2 基于 NB 分布模型进行假设检验,对 DEG 采用经验贝叶斯程序。目前 limma 由于特定的分布模型假设,在 RNA 计数分析中误差较大,虽然 edgeR 和 DESeq2 都利用贝叶斯模型对过度离散进行归一化,但后者通过数据集 reads 的平均值和异常值检测促进了 CPM 阈值的筛选,分析效果更好。
单细胞转录组 DEG 按照时间和分析方法大致可以分为早期零值参数检验、非参数检验和其他方法。由于 scRNA-seq 数据中存在大量零数,早期的方法大多基于此观察做参数检验,例如 Monocle (Trapnell 等,2014)、SCDE (Kharchenko 等,2014)、MAST (Finak 等,2015)、scDD (Korthauer 等,2016)、D3E (Delmans and Hemberg, 2016)、TASC (Jia 等,2017)、DEsingle (Miao 等,2018) 和 HIPPO (Kim 等,2020b)。对以上一些方法的评测表明,虽然它们在单细胞数据集的分析中普遍取得了不错的效果,但对于批量数据(Soneson and Robinson, 2018)相比 DEA 方法并没有明显的性能提升。对于不同的数据集,有可能没有最好的分布模型,因此一种替代解决方案是考虑非参数 DEA 方法。
非参数检验或无分布检验不需要对数据分布形式做事先假设,因此适用于多数据集的分析,常用方法有 Swish(Zhu 等,2019a)、IDEAS(Zhang 等,2022d)、ccdf(Gauthier 等,2021)、distinct(Tiberi 等,2022)。Swish 通过 Salmon Gibbs 评估转录本水平,然后用 Mann-Whitney Wilcoxon 检验计算 FC 值。IDEAS 是一种使用 Jensen-Shannon 散度(JSD)或 Wasserstein 距离(Was)进行基因差异表达测量的伪 F 统计量检验,P 值由基于 PERMANOVA 的距离测试器基于核的回归生成。Ccdf 是一种依赖条件累积分布函数的条件独立性检验,通过多元回归模型预测 DEG。Distinct 提出了一种分层非参数置换方法,使用经验累积分布函数 (ECDF) 的总距离进行 DEG 识别。替代方法包括深度学习策略 MRFscRNAseq (Li 等,2021a)、基于拟时序推断的 PseudotimeDE (Song and Li, 2021)、基于非预聚类的 singleCellHaystack (Vandenbon and Diez, 2020)、基于多重评分的 MarcoPolo (Kim 等,2022)。建议不同的单细胞转录组数据集应采用数据特定的 DEGs 检测策略,以优化 DEGs 分析,基于 scCODE 工作流程,可以使用涉及 CDO(DE 基因顺序)和 AUCC(一致性曲线下面积)的指标找到最优化的 DEGs 方法(Zou 等,2022)。此外,研究方法在不同的研究背景下会有特定的研究取向,例如在给药后的剂量反应研究中,DEGs 分析、LRT 线性检验和贝叶斯多组检验均比其他方法有更好的结果(Nault 等,2022)。
9、可视化
单细胞转录组数据分析可视化是指将上述分析结果以图形的形式直观地呈现,ggplot2 是 R 中最广泛的可视化工具,在 R 中被广泛使用,可以大大增强绘图能力(Wickham,2009)。ARL 是另一个专门显示标记基因关联图并可显示其在每个簇中的特征的 R 包(Gralinska 等,2022)。此外,还有其他专门用于标记基因可视化的包,如 Complex Heatmap,本文不再详细介绍。HVG 可视化通常以火山图的形式呈现,默认情况下,图的左侧和右侧部分分别是代表性不足的基因和代表性过高的基因,而中间是恒定基因。Enhanced Volcano 是一个专门用于绘制火山图的 R 包,默认情况下也可以使用 ggplot2 来获得更好的结果。簇可视化通常以 PCA 图、t-SNE 图和 UMAP 图呈现,但值得注意的是,可视化的结果非常具有欺骗性,因为一些小的细胞亚群可能代表 UMAP 图中显示的大量细胞。为了解决这些问题,提出了 den-SNE/densMAP、j-SNE/j-UMAP 等改进方法(Macqueen,1967;Marco 等,2014)。此外,FastProject 可以输出注释簇的 2D 显示,PieParty 可以在簇 2D 图中为每个基因绘制颜色图(DeTomaso and Yosef,2016;Kurtenbach 等,2021)。
同时,单细胞转录组数据的交互式可视化是目前的热门领域,诸如 Single Cell Explorer 等软件可以一定程度上实现交互式可视化,但仍需增加交互自由度,以提供更全面的单细胞转录组数据 3D 呈现(Cakir 等,2020;Feng 等,2019)。为此,CellexalVR 利用 VR 理论进行交互可视化;CellView 是一个基于 Web 的工具,包括用于不同用途的探索选项卡、共表达选项卡、子簇分析选项卡模块;Cellxgene VIP 是一个基于 cellxgene 框架的插件,并扩展到基于多个模块组合的 ST 数据的交互式可视化(Bolisetty 等,2017 ; Legetth 等,2021 ; Li 等,2022f)。
10、单细胞模拟
随着单细胞转录组方法的不断扩展,基准测试成为了重要挑战,关键问题是需要稳定可靠的数据,因为单细胞转录组的直接测序可能缺乏基本事实。真实的单细胞模拟数据为基准测试提供了已知的事实,允许使用真实数据进行训练,同时匹配实际数据的特征。此外,模拟数据比真实数据提供了更大的灵活性,使分析师能够根据特定的测试方法调整诸如 dropout rate 等参数。
Splatter 是一个两步模拟框架,首先模拟来自真实数据的估计参数,然后合并来自用户的额外参数(Zappia 等,2017)。其六个预先设计的管道模块接口确保了数据生成的可重复性。最近的更新侧重于专业化和泛化。在专业化领域,splaPop 生成具有遗传效应(数量性状基因座)的人口规模数据,而 dyngen 模拟动态细胞过程,如发育轨迹(Azodi 等,2021;Cannoodt 等,2021)。在泛化领域,Li 的团队介绍了理想模拟的六个概念,包括真实性、基因的保存、基因相关性的捕获、稳健性、参数可调性和效率(Song 等,2023b;Sun 等,2021)。随后,scDesign2 提出来满足所有 6 个属性(Sun 等,2021),接着是 scDesign3,解决单细胞组学统计模拟的空白(Song 等,2023b)。模拟准确性和透明度的提高增强了不同单细胞数据处理方法之间的基准测试,指导选择最合适的方法以满足特定数据和许可需求。