超详细 | 生物医学研究和临床应用中scRNA-seq的数据分析指南

随着高通量scRNA-seq（包括临床样本）能力的扩大，对这些海量数据的分析能力已成为进入该领域研究人员的必备技能。近日，《Military Medical Research》发表了一篇综述文章，回顾了典型scRNA-seq数据分析的工作流程，总结每个分析步骤的当前方法，并针对某些特定的分析任务和方法指出了建议和注意事项。 

scRNA-seq的典型数据分析步骤通常可分为三个阶段：  1）原始数据处理和质量控制；2）适用于几乎所有scRNA-seg数据集的基本数据分析：数据标准化和整合、特征选择、降维、细胞聚类、细胞类型标注和标记基因识别；3）应针对特定研究场景定制的高级数据分析：轨迹推断、细胞间通讯分析、调节子推断和TF活性预测以及代谢分析。

实验设计

scRNA-seq 实验需要精心设计，以优化解决科学问题的能力。在开始数据分析之前，需要收集以下与实验设计相关的信息。 (1) 物种。 对于生物医学研究和临床应用，通常收集来自患者的人体样本进行测序。在某些情况下，为了研究潜在的分子机制，还使用小鼠和其他模式生物。由于人类和其他物种的基因名称和相关数据资源不同，因此明确数据分析的物种非常重要。 (2)样本来源。 根据科学问题和样本可及性，不同研究的样本类型可以不同。例如，为了研究肝细胞癌等实体瘤，从患者身上收集肿瘤活检和肿瘤周围样本用于病例对照设计。尽管上述设计在某种程度上是可行的，但外周血单核细胞 (PBMC) 更容易获得并广泛用于 scRNA-seq。此外，来自患者的类器官细胞通常用于研究个人遗传变异对特定器官发育的影响，这也可能是特定疾病的起源。了解样本来源有助于特定分析，例如细胞聚类和细胞类型注释。 (3)实验设计。 为了研究疾病的发病机制和特定治疗的有效性，通常采用病例对照设计，如肿瘤与肿瘤周围设计。对于COVID-19等疾病，无法从同一患者身上获得正常样本，因此将年龄和性别匹配的健康人作为对照组。为了控制患者和对照组之间可能的协变量，需要仔细考虑每组中的个体数量。在（前瞻性）队列研究中，样本量通常相当大，因此 scRNA-seq 不能应用于来自个体捐赠者的每个样本；在这种情况下，通常会应用嵌套病例对照研究和样本多重分析。 一般来说，数据分析策略需要根据实验设计的类型进行调整。

原始数据处理和质量控制

原始数据处理

原始数据处理步骤包括：测序read QC、read mapping、细胞解复用和细胞唯一分子标识符（UMI）计数表生成。目前随着 scRNA-seq 平台的发布提供了标准化的数据处理管道，例如 10 × Genomics Chromium 的 Cell Ranger 和用于 Singleron 系统的 CeleScope，替代工具包括UMI-tools、scPipe、zUMIs、celseq2、kallisto bustools 和 scruff 也可用于此过程。无论如何，我们不建议在个人电脑上处理原始数据，因为这些管道需要大量的计算资源，并且针对高性能计算架构进行了优化。 第三方公司通常提供经过处理的数据，包括 UMI 计数矩阵和 QC 指标，这使研究人员能够专注于下游数据分析以解决科学问题。

QC和双胞去除

细胞质量控制的目的是确保所有被分析的“细胞”都是单个完整的细胞。受损细胞、死亡细胞、应激细胞和双重细胞需要丢弃。细胞质量控制最常用的三个指标是：总UMI计数（即计数深度）、检测到的基因数量和每个细胞条形码的线粒体衍生计数分数。Cell Ranger和CeleScope通常执行第一轮细胞质控，通过检查scRNA-seq库中计数深度的分布，将潜在的真实细胞与背景细胞条形码区分开来。

需要注意的是，当损坏的细胞或细胞碎片在文库中占据相当大的比例时，有效细胞的最小计数深度阈值就很难确定。可能的解决方案包括同时考虑多个 QC 指标，以及应用更复杂的方法来排除背景和低质量细胞。 通常，检测到的基因数量少和计数深度低表明细胞受损，而线粒体衍生计数的比例高则表明细胞正在死亡。相比之下，检测到的基因太多和计数深度太高可能表明存在双联体。虽然像 Seurat 和 Scater 这样的 R 包实现了促进细胞QC的功能，但QC指标的阈值在很大程度上取决于所研究的组织、细胞解离方案、文库制备方案等。参考具有类似实验设计的文献将有助于确定阈值，研究人员也可检查QC指标的联合分布。值得注意的是，编码核糖体蛋白的基因的累积表达不是典型的QC指标，因为核糖体蛋白表达的变化可能具有生物学意义。

此外，在质量控制步骤中需要考虑和控制各种污染源。最近开发了用于估计和消除此类污染的方法和工具，包括SoupX、DecontX、FastCAR和 CellBender。

在高通量 scRNA-seq 实验中，观察到双细胞的比率并不少见，最高可达细胞条形码的 40%。在最近的一项研究中，使用一组综合的合成和真实数据对可用的计算doublet检测方法进行基准测试，推荐使用工具 Doubletfinder，因为它在下游分析中实现了最高的检测精度和最佳性能。

基本数据分析

数据归一化

每个细胞的总 UMI 计数的可变性取决于一系列技术和生物学参数。为 scRNA-seq 量身定制的标准化方法，包括SCDE和MAST，可以在 scRNA-seq 数据的差异表达分析中专门模拟dropout events。另一种方法Scran通过汇集具有相似基因表达谱的细胞来克服比例因子估计的问题（受太多零计数的影响）。此外，Census 在没有加入对照的情况下估计每个细胞的RNA分子总数，并将这些估计值用作比例因子。虽然Vallejos等人进行的模拟研究表明Scran的池化策略在缩放因子估计方面优于比较工具，但TPM-/count 深度缩放方法在实践中得到了广泛应用。为了更好地稳定方差，Seurat团队最近开发了SCTransform，它将正则化负二项回归应用于scRNA-seq数据归一化和方差稳定。

一些已知的生物效应，例如细胞周期和细胞应激（以线粒体基因的过度表达为特征），可能会阻碍对特定感兴趣的生物信号的表征。因此，针对已知生物学规范化或校正表达谱可能有助于解释数据。 如纠正细胞周期的影响可以改善发育轨迹的重建。可以通过对相关生物学特征（例如细胞周期分数）进行评分，然后根据Seurat中实施的计算分数进行简单线性回归，来实现生物效应的程序。此外，专用工具，如scLVM/f-scLVM和cgCorrect也可用于此目的。值得注意的是，校正一种特定分析（如细胞分化）的生物效应可能会无意中阻碍另一种分析（如细胞增殖）的信号；在为特定分析任务选择数据标准化策略时，应谨慎。

数据整合

如“实验设计”部分所述，生物医学研究通常进行病例与对照比较，通常下游分析之前需要进行整合，这个时候要考虑去除批次效应。批次效应校正方法可分为以下几类：1）为批量表达分析开发的工具，包括ComBat和limma；2） 基于高维基因表达空间或其子空间中的相互最近邻（MNN）的方法，如mnnCorrect、fastMNN、Scanorama和BBKNN；3） 尝试在降维空间中使细胞与相关/共享特征对齐的方法，包括CCA、Harmony和LIGER；以及4）基于深度生成模型（如scGen）的方法。此外，根据整合锚点的选择，算法也可以分为不同的类型，例如基因组特征作为锚点和细胞作为锚点。

最近的一项研究在5种不同的整合场景下，在10个数据集上比较了当前可用的14种批次效应校正方法：Seurat 3中实现的Harmony、LIGER、CCA根据综合性能推荐。 Harmony 运行速度比其他工具快，适合初步探索；Seurat3/4-CCA 在混合不同批次的细胞方面表现适中，而 LIGER 在批次混合方面做得最好，有时会以牺牲细胞类型纯度为代价。 值得注意的是，如果想要评估批次效应校正的有效性或评估数据中批次效应的程度，可以通过比较基于批次效应校正分析的聚类或可视化结果，以及直接合并来自多个样本的细胞的结果如Seurat 中的合并函数），并通过计算测试指标，例如kBET。

功能选择

“特征选择”一词是从统计学和机器学习领域借来的，用来描述为下游分析选择生物信息基因的过程。这个过程通常是无监督的，这意味着不需要与细胞类型或其他感兴趣的生物过程相关的信息。考虑到scRNAseq数据中相对较高的噪声水平，特征选择通常识别具有比技术噪声更强的生物变异性的基因。一些不利的协变量（如批次效应）可能会扭曲高度可变基因（HVG）识别。因此，HVG选择应在协变量校正后进行。在存在批量效应的情况下，也可以在数据集成之前在单个样本中进行特征选择。

降维和可视化

选择1000–5000个HVG后，表达数据的维度仍然很高，为此，可以通过降维技术进一步降低表达矩阵的维数。广泛使用的降维方法包括PCA、NMF、MDS、t-SNE和UMAP。PCA是一种用于降维和去噪的通用技术，然而，PCA无法在分析中考dropout events，这导致开发了几种新方法。ZIFA是一种基于因子分析的方法，它明确地模拟了dropout特征，并优于比较方法。与PCA相似，NMF是一种用于降维的线性投影方法，并在基于scRNA-seq的细胞聚类中表现出稳健的性能。

对于可视化，非线性降维方法更适合，它允许在二维/三维空间中进行全局非线性嵌入。MDS是一种非线性降维方法，它保持了原始空间中细胞之间的距离。然而，MDS不能扩展到大规模scRNA-seq数据。新的证据表明，t-SNE和UMAP更适用于scRNA-seq数据，这些数据已广泛用于单细胞分析，用于数据可视化和细胞群体识别。然而，t-SNE通常受到限制，例如大规模scRNA-seq数据集的计算时间慢，并且全局数据结构没有得到保存。UMAP具有上述两个方面的优势，目前成为最流行的降维选择。

细胞亚群的鉴定

scRNA-seq的关键应用之一是基于细胞聚类或分类来确定细胞亚群。NMF被用于在单细胞转录组数据中分离亚群，并已被证明以更高的准确性和鲁棒性优于PCA。同样的，开发 SinNLRR 是为了通过非负和低秩表示提供强大的基因表达子空间聚类。

最先进的聚类方法，如k-means算法，也已应用于scRNA-seq数据集，并基于此应用，开发了SC3方法。scRNA-seq中另一类常用的细胞聚类方法是基于细胞最近邻网络的聚类检测方法，并在Seurat R包中采用和实施。此外已开发了多种细胞聚类方法，如BackSPIN等。然而，两项独立的基准测试研究表明，SC3和Seurat中的聚类方法表现相似，优于所有其他比较方法。

相似性或距离度量对于scRNA-seq中的细胞聚类至关重要，这可能是实验平台或特定样本特有的。研究表明，与无监督聚类方法相比，有监督的细胞类型识别方法较少受到批量效应、细胞类型数量和细胞群体组成失衡的影响。 CellAssign、scmap、SingleR、CHETAH和SingleCellNet都是有监督的方法。尽管监督方法具有明显的优势，但无监督方法通常更擅长识别未知细胞类型并且具有更高的计算效率。因此，在 Seurat 中实现的聚类方法具有最佳的整体性能，并被建议作为细胞类型识别的首选。

单细胞聚类分析的另一个重要问题是稀有细胞类型的检测，这些细胞在复杂疾病中起着重要作用，但丰度较低。RaceID、GiniClust、SINCERA 和 DendroSplit 是专门设计用于在 scRNA-seq 数据分析中识别稀有细胞类型的聚类算法。

细胞类型注释

将细胞身份分配给细胞亚群，这一过程称为细胞类型注释，是 scRNA-seq 数据分析中的关键步骤。细胞类型的手动注释非常耗时且可能具有主观性。因此，已经开发了用于自动细胞类型注释的新兴计算工具。这些计算方法通常可以分为三大类。

第一类是基于标记基因的，这依赖于公共数据库或文献中细胞类型特异性标记的可用性。CellMarker和PanglaoDB是常用的在线资源，此外，开发了TF标记物数据库，用于为人类提供细胞或组织特异性TF和相关标记物。同时，已经开发了许多工具来将标记基因用于细胞类型注释，例如ScType、scSorter、SCINA、scCATCH和CellAssign。

第二组方法是基于参考转录组的，通过搜索查询数据和参考数据之间的最佳相关性，使用细胞类型标记的scRNA-seq数据集作为细胞类型注释的输入。该组的流行工具包括CHETAH、scmap、scMatch和SingleR。值得注意的是，参考转录组学数据的全面性对于这组方法至关重要。常用的参考数据来自 Blueprint、Encode 和 Human Primary Cell Atlas。

第三组利用基于监督机器学习的方法，其中由标记参考训练的分类器然后被应用于预测未标记细胞的细胞类型。例如SingleCellNet、scPred和scClassify等。此外，基于OnClass还可以通过识别由细胞本体和未标注的细胞产生的低维嵌入中的最近细胞类型来准确标注训练数据集中缺失的细胞类型。

细胞类型标注的自动化方法已应用于广泛的生物医学研究，包括癌症研究。然而，最近的一项基准研究表明，在不同的场景下，每种计算方法都比其他方法具有特定的优势，因此临床用户很难选择合适的工具。集成来自多个工具的注释结果可能是上述问题的解决方案，并可能实现更精确的细胞类型注释。因此，ImmCluster最近被开发用于免疫细胞聚类和注释，它集成了七种基于参考和四种基于标记基因的计算方法，并由人工选择的标记基因集支持。比较研究表明，ImmCluster提供了比单个方法更准确和稳定的细胞类型注释。

标记基因鉴定

特定细胞簇或细胞类型的标记基因是表征其功能的重要资源。反过来，标记基因也可以用于细胞类型注释。常用方法包括Seurat、SINCERA、SC3、MAST、SCDE和DEsingle等。

还有一类方法，它在细胞聚类过程中同时识别细胞特异性基因，而不是此后的一个步骤。例如BackSPIN、ICGS、DendroSplit、scVDMC和BCseq等。

关于 scRNA-seq 中 DEG 工具的最佳选择，最近的一项研究比较了 36 种方法，发现所比较的方法之间存在根本差异。有人指出，低表达基因的预过滤可能有助于 DEG 分析，用于 bulk RNA-seq 分析的方法通常与专门为 scRNA-seq 开发的方法具有相当的性能。总体而言，非参数 Wilcoxon 秩和检验在大多数应用场景中排名靠前，但复杂的实验设计除外。

高级数据分析

功能富集分析

为了便于解释和组织每种细胞类型中识别的标记基因，通常进行功能富集分析。为批量转录组学开发的计算方法可以轻松应用于此分析，例如DAVID。这类分析需要在统计意义上进行严格的切割，以定义标记基因；相比之下，广泛使用的基因集富集分析（GSEA）是一种无边界的方法，为了促进 GSEA 分析，MSigDB数据库提供了一系列带注释的基因集，包括通路和标志基因特征。

除了基于标记基因或两组细胞之间的差异表达执行功能注释的上述场景之外，该分析也可以在单细胞水平上进行。常用方法包括Vision、Pagoda2、AUCell、SCSE和JASMINE，其更适合于scRNA-seq中的特征评分。此外，这些特征评分方法也可用于通路活动推断。

轨迹推断和RNA速度分析

轨迹推断、伪时间估计和RNA速度建模都有助于揭示细胞分化或激活过程中的分子特征和调控机制。轨迹推断是过去几年的热门研究领域，开发了大约一百种计算工具。根据轨迹的类型，轨迹推断方法也可分为不同的类别，包括线性方法[例如SCORPIUS、TSCAN、Wanderlust]、分叉方法[例如DPT、Wishbone]、多分叉方法[例如FateID、STEMNET、MFA]、树方法（例如Slingshot、scTite、Monocle）和图方法[例如PAGA、RaceID、SLICER]。基于最近的一项基准研究，给出了实际应用指南，以便生物医学研究人员可以根据数据中预期拓扑的先验知识选择适当的方法；否则，建议使用PAGA、Monocle、RaceID和Slingshot进行初步调查。

捕获转录组动力学的另一种方法是使用RNA 速度，相关方法包括Velocyto和scVelo等。最近开发的方法将RNA速度与轨迹推断相结合，例如CellRank和CellPath等。

细胞间通讯（CCC）

CCC 事件在生物体发育和体内平衡以及疾病的发生和发展中起着重要作用。为此，已经开发了各种计算工具来使用 scRNA-seq 数据推断 CCC。细胞之间的通讯通常取决于配体-受体 (LR) 相互作用，通常通过 LR 共表达来量化。

为了促进上述研究，已知的配体-受体相互作用 (LRIs) 已被手动整理并存储在数据库中。迄今为止，LRIs数据库已经相当多了，包括 CellPhoneDB、ICELLNET、CellTalkDB、SingleCellSignalR 和 Omnipath。

最近，基于scRNA-seq数据预测CCC的计算方法不断发展。CCC推理工具可根据其特殊性分为三大类：即基于网络的（NicheNet、CCCExplorer、scConnect和NATMI）、基于机器学习的（SingleCellSignalR、SoptSC和PyMINEr）和基于空间的（CellTalker、CellChat）方法。此外，参考成分分析 (RCA)-CCA、线性回归和决策树分类器也用于 CCC 预测。

识别CCC的新兴计算方法提高了我们对疾病发展微环境的理解。然而，所有的方法都依赖于LRIs的先验知识和统计或机器学习模型来预测潜在的CCC事件。另外，选择LRIs资源和预测方法可能会产生不同的结果，但选择对结果的影响在很大程度上是未知的。为了解决这个问题，最近的一项研究系统地比较了CCC推断的16种资源和7种方法，以及比较方法的共识。比较表明，不同的LRIs资源涵盖了不同比例的集体先验知识，预测的CCC在很大程度上彼此不一致，这表明需要继续努力改进CCC推理资源和工具。

调节子推断和TF活动预测

识别调节子的一个重要资源是TF目标数据库。AnimalTFDB、JASPAR、TRRUST、KnockTF和Cistrome DB揭示的转录调控关系是广泛应用的TF注释数据库，涵盖了大多数人类和小鼠TF。

结合单细胞基因表达和综合的TF靶信息，已经开发了许多推断调节子和TF活性的方法，包括共表达分析（如WGCNA）、SCENIC、SCODE、SINCERITIES和机器学习的方法（如DeepDRIM）等。

代谢分析

基于scRNA-seq的代谢分析的计算工具可分为两大类：基于途径的分析和基于通量平衡分析（FBA）的方法。对于第一类，通常使用标准功能富集分析方法。特别是，R 包 scMetabolism 提供了一个集成框架，用于定量分析 scRNA-seq 中的代谢途径活动，能够解释丢失，并与为单细胞功能富集分析设计的多种工具兼容，包括 ssGSEA、Vision、 和 AUCell。scFBA是第一个结合scRNA-seq数据和FBA来估计单细胞通量的计算工具后来提出了Compass和scFEA。

本文中所涉及的软件及代码均可在如下链接获取：https://github.com/WXlab-NJMU/scrna-recom

参考文献Su M, Pan T, Chen Q Z, et al. Data analysis guidelines for single-cell RNA-seq in biomedical studies and clinical applications[J]. Military Medical Research, 2022, 9(1): 1-24.