写在前面
自2009年第一篇单细胞转录组文章发表以来，经过十多年的发展，各研究机构和生物公司相继推出单细胞RNA测序(scRNA-seq)服务，这使得单细胞测序已经从“王谢堂前燕”飞入了”寻常百姓家“，逐渐成为了科研人员使用的常规技术。但是，测序公司提供的往往只是基础的数据处理，想要根据自己的课题进行特定的高级分析，还需要我们有一定的生信基础，至少要了解在拿到数据之后，能做的分析有哪些，这对很多生物医学研究人员和临床医生来说，无疑有一定的难度。
最近，小编发现了一篇发表在Military Medical Research期刊上的综述，关于MMR这本国产杂志，小编不得不多说几句，MMR创刊于2014年，2019年被SCI收录，2020年首个IF为2.3，2022年直接飙升到了34.9分，多少可能是有点猫腻的（小声哔哔）。然鹅目前来看，该期刊走的还是优质路线，专业认可度也是较高的。好了，回归正文，在这篇综述中，作者针对有科研需求但更偏向于临床应用的生物人员，总结了生物医学方面单细胞数据分析的流程以及每一步的算法进展，并推荐了更适合生物医学的软件工具。这篇综述更nice的地方在于，它所有的软件环境设置和数据分析的脚本都是直接奉上了，链接在：https://github.com/WXlab-NJMU/scrna-recom，对于有编程基础的小伙伴，完全可以进行复现。
单细胞分析流程
scRNA-seq数据的分析步骤可分为三个阶段：1、原始数据处理和QC；2、几乎适用于所有scRNA-seq数据集的基础分析；3、针对特定科学问题量身定做的高级数据分析。基本的数据分析步骤包括数据归一化和整合、特征选择、降维、细胞聚类、细胞类型注释和标记基因识别，高级数据分析通常是指轨迹推断、细胞通讯分析、调控子和TF活性预测以及代谢通量分析。
实验设计
一个科研项目始于科学问题的提出，而实验设计则关乎到项目成立与否，也是后续分析的前提和基础。很多临床医生，往往是能够在工作中发现问题，但至于能不能研究、怎么研究，就摸不着头脑了。其实，对于生物医学领域来说，前期最关键的就是样本的收集，通常是指人体组织或血液，可能还会涉及小鼠等模式生物。样本的临床信息同样需要采集，比如肿瘤分期分级、治疗记录、并发症以及各项检查结果等，这个就是多多益善了。为了研究疾病的发病机制和特定治疗的有效性，通常采用病例对照设计，比如肿瘤活检和癌旁组织，或者患者与健康人，在单细胞采样中，最重要的就是要保证细胞的活性，这一部分，就需要实验人员的参与了。
数据处理
1.1、原始数据处理
原始数据处理的步骤包括：测序reads QC、reads映射、文库拆分和count矩阵的生成。各大scRNA-seq平台一般都会提供标准化的数据处理管道，例如10x Genomics Chromium的Cell Ranger和Singleron的CeleScope，一些替代工具如UMI-tools、scPipe、zUMIs、celseq2、kallisto bustools和scruff等也可用于原始数据的初步处理，对结果的影响并不是很大。如果是送去公司测序，我们拿到手的都是他们处理过的数据，包括count矩阵和质控报告，可以直接进行下游的数据分析。
1.2、细胞QC和双细胞去除
对细胞QC的目的是确保所有被分析的细胞都是单个且完整的细胞，因此需要过滤掉那些受损细胞、垂死细胞、应激细胞和双联体。细胞QC最常用的三个指标是：UMI总数、检测到的基因数量和线粒体基因比例，阈值的设置没有统一的标准，在很大程度上取决于所研究的组织、细胞解离方案、文库制备方案等，我们可以参考其他文章，根据QC结果来定。不过，保险起见，最好是尽可能的多保留细胞。需要注意的是，当受损的细胞或细胞碎片在文库中占据相当大的比例时，count阈值就很难确定，需要考虑多个QC指标，或者用更复杂的方法来去除背景和低质量细胞。通常来说，线粒体基因占比高表明细胞正处于死亡状态，当检测到的基因数量少、测序深度低时意味着细胞受损，反之，则考虑双细胞。在高通量scRNA-seq实验中，观察到高双联体比例的情况并不少见，目前也有很多方法用以检测，如Scrublet、doubletCells、bcds、DoubletDetection、DoubletFinder、Solo、DoubletDecon等。此外，在质控时还需要考虑和控制各种污染源，例如，源自PBMC和实体组织的文库可能会被红细胞污染，因此高表达血红蛋白基因的细胞通常会被丢弃，再比如无细胞和环境RNA，可以利用SoupX、DecontX、FastCAR和CellBender等软件去除此类背景信号。在最近的一篇方法测评文章中，作者发现Doubletfinder在下游分析中表现出了最高的检测精度和性能。
基础分析
1.1、表达归一化
每个细胞测得的count数取决于一系列技术和生物学因素，技术因素包括RNA捕获、逆转录、cDNA扩增和测序深度等，生物学因素则是指细胞大小和细胞周期阶段。因此，在比较细胞间的基因表达谱时，通常用的都是RNA相对表达丰度。由于单细胞表达矩阵的稀疏性，bulk测序使用的归一化方法，如TPM、TMM、上四分位数法和DESeq等都不适用。scRNA-seq数据使用的归一化方法，包括单细胞差异表达(SCDE)、基于模型的单细胞转录组学分析(MAST),以及通过汇集具有相似基因表达谱的细胞进行归一化的Scran。在实操中，我们会先将单细胞数据Normalization，以消除技术差异，确保同一基因在不同样本中的表达具有可比性，然后通过Scale添加z-score计算，使后续分析不受极值影响。Seurat团队最近开发的SCTransform，将正则化负二项式回归用于scRNA-seq数据的标准化和归一化处理。
一些已知的生物学效应，例如细胞周期和细胞应激，可能会影响我们感兴趣的特定生物学信号的表征，这时，针对已知的生物学过程来标准化基因表达谱有助于我们更好的解释数据。例如，校正细胞周期的影响可以改善发育轨迹的重建，不过，在校正一种生物学效应时，可能又会无意中阻碍另一种生物学信号的表征，所以在选择策略时需要小心谨慎。
1.2、数据整合
当研究样本来自不同医疗中心时，需要在下游分析前进行整合，对于使用类器官的研究，在不同时间点采集的细胞也需要整合，即在保留生物变异的同时，去除批次效应，在scRNA-seq中，这种批次效应可能是非线性的。目前，已有的校正方法可以分为以下几类：1) 针对bulk数据开发的工具，包括ComBat和limma；2) 基于高维基因表达空间或其子空间中相互最近邻（MNN）的方法，如mnnCorrect、fastMNN、Scanorama和batch balanced k nearest neighbors（BBKNN）; 3) 将细胞与降维空间中的相关/共享特征对齐的方法，包括CCA、Harmony和LIGER; 4) 基于深度生成模型的方法，如scGen。Tran等人使用10多个数据集比较了14种整合方法，发现Harmony、LIGER和CCA的整体性能最佳，根据实际经验，推荐依次用Harmony、Seurat3/4-CCA 和LIGER。Harmony运行速度最快，适合初步探索，Seurat3/4-CCA表现适中，而LIGER将数据混合得最好，相应的可能会损失细胞类型纯度。为了评估批次效应的程度或校正效果，我们可以将校正后的数据聚类并可视化，与直接合并样本的结果进行比较，通过计算kBET等测试指标来评估。
1.3、特征选择和降维
细胞QC是为了去除背景细胞和有问题的细胞，特征选择则是为了筛选有代表性的基因。在人类基因组中，超过20,000个基因被注释，但并非所有基因都可以提供表征细胞间异质性或区分细胞类型/状态的信息，而且在单细胞数据中，一个细胞只表达几千甚至几百个基因。因此，考虑到scRNA-seq数据的高噪声比，特征选择通常会识别高度可变基因(HVG)。为下游分析选择的基因数量在理论上取决于所研究样本中细胞组成的复杂性，通常HVG数量在1000到5000之间,Seurat默认为2000。有研究表明，下游分析对HVG的数量并不敏感。特征选择之后，需要通过降维技术进一步降低表达矩阵的维度，并保留细胞间变异性的生物学信息。广泛使用的降维方法包括PCA、NMF、MDS、t-SNE和UMAP，PCA已广泛应用于scRNA-seq数据分析，通过将原始表达矩阵线性投影到其子空间，按重要性顺序给出主成分 (PC)，NMF与PCA 类似，基于线性投影来降维，在scRNA-seq的细胞聚类中表现出稳健的性能。对于单细胞数据的可视化，常用的方法是t-SNE和UMAP，然而，当数据量比较大时，t-SNE运行速度会很慢，此时，UMAP更胜一筹。
1.4、识别细胞亚群
单细胞转录组学的关键应用之一就是对细胞聚类或分类，以识别细胞亚群。
PCA常用于bulk RNA-seq，但少数的PC可能无法区分细胞亚群中基因表达的生物学差异，而NMF则能很好的分离单细胞转录组数据中的亚群，并且已被证明比PCA具有更高的准确性和稳健性。此外，基于k-means算法开发的SC3，以及基于最近邻的Seurat包中提供的聚类方法也广泛用于scRNA-seq数据集。以上都是无监督聚类，还有一些监督聚类方法，包括CellAs-sign、scmap、SingleR、CHETAH和SingleCellNet，该类方法依赖于具有已知细胞类型注释的参考数据库，受批次效应、细胞类型数量和细胞群组成不平衡的影响较小。但是，无监督方法可以识别未知的细胞类型，计算效率也更高，因此，通过Seurat实现单细胞聚类是多数情况下的首选。此外，还有一些方法专门用于识别丰度较低的稀有细胞类型，如RaceID、GiniClust、SINCERA和DendroSplit聚类算法。
1.5、细胞类型注释
众所周知，细胞类型注释分为手动注释和自动注释，手动注释非常耗时且具有主观性，对注释人员的专业要求也很高，在这里，我们重点了解一下自动注释算法。
这些自动注释方法可以大体分为三大类，第一种是基于marker基因的，依赖于公共数据库或文献中细胞类型的特异性marker，CellMarker和PanglaoDB是常用的单细胞数据库，用于存储人类和小鼠组织中多种细胞类型的标记。此外，TF-Marker数据库提供了人类细胞或组织特异性转录因子和相关marker。同时，已有许多使用marker基因进行细胞类型注释的方法，例如ScType、scSorter、SCINA、scCATCH和CellAssign。第二类方法是基于参考转录组的，它使用细胞类型标记的scRNA-seq数据集作为细胞类型注释的参考，通过搜索查询数据和参考数据之间的最佳相关性进行注释，包括CHETAH、scmap、scMatch和SingleR。其中，SingleR根据参考数据中细胞类型之间的HVG，将每个未注释的单细胞转录组与已知细胞类型的参考转录组相关联，以迭代方式的进行分配。第三组方法基于监督机器学习，然后应用由参考marker训练的分类器来预测未注释细胞的细胞类型,例如，SingleCellNet使用多类随机森林分类器，ACTINN使用人工神经网络，scPred使用支持向量机SVM，而scClassify使用集成学习进行细胞类型注释。最近的一项基准研究表明，在不同的情况下，不同的计算方法都具有特定的优势，因此整合多个工具的注释结果可能实现更准确的细胞类型注释。比如，最近开发的ImmCluster集成了七种基于参考和四种基于标记基因的计算方法，用于免疫细胞聚类和注释，其结果比单个方法更准确和稳定。
1.6、识别marker基因
Marker基因不仅可以用于细胞类型的注释，也可用于表征特定细胞簇或细胞类型的功能。识别marker基因的典型方法是基于统计检验识别细胞簇中的差异表达基因(DEGs)，例如scRNA-seq分析软件Seurat和SINCERA，使用非参数Wilcoxon秩和检验来识别特定细胞类型的高表达基因。有研究表明，Wilcoxon秩和检验的假阳性率低于DESeq2和edgeR等基于测序的分析方法。此外，SC3方法采用非参数Kruskal-Wallis检验来比较两组以上的细胞，还有MAST、SCDE和DESingle等许多其他方法可以来识别细胞特异的marker基因。有一类方法在细胞聚类过程中同时识别细胞特异性基因，包括BackSPIN，它在对细胞进行聚类时将高表达的基因聚集在一起，还有ICGS通过识别引导基因进行迭代聚类，以及DendroSplit在识别子簇的同时计算了标记基因的显著性水平。对于这么多种识别scRNA-seq中DEG的方法，也有研究做了比较，总的来说，非参数Wilcoxon秩和检验在大多数情况下都是不错的选择。
高级分析
在对测序数据做了基本的分析之后，想要挖掘有意义的生物信息，还得靠高级分析。目前，应用最广泛的单细胞分析策略有功能富集、拟时序、细胞通讯和转录因子分析。富集分析是指对与生物学问题显著相关的基因功能类别进行统计分析，相信小伙伴们都已经做过N次了，拟时序分析可以根据细胞之间表达模式的相似性对单个细胞进行轨迹排序来模拟细胞动力学，细胞通讯分析通过确定不同细胞类型中受体和配体的表达和配对来推断不同细胞之间的相互作用，转录因子分析则是识别转录因子和潜在靶基因之间的共表达模块。
1.1、功能富集分析
为了解释每种细胞类型中识别到的标记基因，我们通常都要进行功能富集分析。在bulk转录组分析中，常用的就是DAVID或者基因集富集分析(GSEA)，估计大家对此也都不陌生了。在单细胞富集分析中，我们广泛使用的是单样本GSEA(ssGSEA)和基因集变异分析(GSVA)，它们类似于GSEA。此外，考虑到scRNA-seq数据的特征，Vision、Pagoda2、AUCell、SCSE和JASMINE等工具被开发出来专门用于单细胞，且通常更适合对scRNA-seq进行特征打分，还可用于推断通路活性。
1.2、轨迹推断和RNA速率分析
在整个发育过程中，细胞会从一种功能“形态”分化到另外一种功能“形态”，不同形态的细胞会表达不同的基因，以实现它们特定阶段的功能。基于单细胞数据分析，我们可以表征这些处于中间形态的细胞，拟时序分析，即轨迹推断根据每个细胞的时序基因表达，将每个细胞按照拟时间排列在对应轨迹上，反映了细胞的发育轨迹或细胞状态转变。在过去几年里，轨迹推断的方法大约开发了一百多种，根据轨迹的类型，可以分为线性方法（如SCORPIUS、TSCAN、Wanderlust），单分叉方法（如DPT、Wishbone），多分叉方法（如FateID、STEMNET、MFA），树形方法（包括Slingshot、scTite、Monocle）和图形方法（PAGA、RaceID、SLICER)。目前，常用的方法有PAGA、Monocle、RaceID和Slingshot。单细胞轨迹分析的另一种方法是RNA速率分析，它基于同一细胞中未成熟（未剪接的）和成熟（剪接的）mRNA丰度之间的比值来获得基因特异性速度，得出可能的细胞状态变化，从而追溯细胞的起源和潜在的命运。第一个RNA速率分析方法是Velocyto，scVelo在其基础上进行了改善。此外，还有一些将RNA速率与轨迹推断相结合的方法，例如CellRank和CellPath。
1.3、细胞通讯(CCC)分析
细胞通讯在生物体发育和体内平衡以及疾病的发生和发展中起着重要作用，CCC分析可以帮助我们了解细胞与细胞之间的互作关系，解析细胞间通信网络，探索肿瘤免疫微环境，挖掘疾病潜在的治疗靶点。
细胞之间的通讯通常取决于配体-受体相互作用(LRIs)，迄今为止，已经搭建了非常多的LRI数据库，包括CellPhoneDB、ICELLNET、CellTalkDB、SingleCellSignalR和Omnipath。对scRNA-seq数据完成注释之后，将其与已知的LRI整合，以计算特定样本的LR分数，从而量化相互作用的可能性。基于LR共表达，LR评分函数可以分为表达阈值、表达相关性、表达产物和差异表达的组合。例如，SingleCellSignalR就是基于LR基因表达的产物。根据特征，CCC分析工具可分为三大类，即基于网络、基于机器学习和基于空间信息。基于网络的方法利用基因之间的连接网络来预测细胞通讯，包括NicheNet、CCCExplorer、scConnect和NATMI。基于机器学习的方法采用了各种类型的机器学习算法，例如SingleCellSignalR、SoptSC和PyMINEr，此外，RCA-CCA、线性回归和决策树分类器也用于细胞通讯的预测。随着空间转录组学的快速发展，许多CCC方法将scRNA-seq数据与空间转录组学和/或图像数据相结合，用于识别CCC，如CellTalker、Squidpy和histoCAT，通过细胞邻近度或邻域分析来研究细胞间通讯，此外，CellChat在预测互作方面的表现也很好。分析结果通常借由热图、环形热图、桑基图和气泡图来进行可视化。
1.4、调控子推断和TF活性预测
转录因子在基因表达调控中起着重要作用，并参与人类的各种生理和病理过程。被同一TF调控的基因集合，称为调控子，基于scRNA-seq数据，我们可以绘制细胞类型特异性调控子的图表，重建每个细胞的调控网络。
调控子的识别离不开TF靶标数据库，AnimalTFDB、JASPAR、TRRUST、KnockTF和Cistrome DB是应用广泛的TF注释数据库。基于这些数据库，构建细胞类型特异性转录调控网络的一种简单方法，是识别上调的转录因子和/或差异表达的转录因子靶基因。结合单细胞基因表达和TF靶标信息，已经开发了许多用于推断调控子和TF活性的方法。共表达分析，例如WGCNA已广泛用于bulk数据中TF调控模块的挖掘，该方法也同样适用于scRNA-seq。SCENIC方法是最早基于scRNA-seq数据进行调控子推断的方法，目前已被用于研究癌症和COVID-19等多种疾病的调控网络。其他方法，包括SCODE和SINCERITIES，利用拟时序信息在scRNA-seq中基于常微分方程或随机微分方程模型，来重建TF-靶基因调控网络。此外，机器学习技术也已应用于转录调控分析，例如SIGNET和DeepDRIM。虽然有许多基于scRNA-seq的基因调控分析方法可供选择，但由于转录调控的复杂性和scRNA-seq数据的信息不足，需要我们对推断的结果进行严格判断，必要的话可以通过实验进行验证。
1.5、单细胞代谢分析
代谢是所有生物过程的核心，代谢失调是包括癌症、糖尿病和心血管疾病等的标志。单细胞代谢组学技术尚不成熟，因此，基于单细胞转录组学的代谢分析更常见。例如，研究人员可以使用scRNA-seq监测关键代谢基因在不同处理下或重要生理/病理过程期间的表达变化。
基于scRNA-seq的代谢分析工具可分为两大类：基于通路分析和基于通量平衡分析(FBA)的方法。对于第一类，通常使用功能富集分析方法，尤其是R包scMe-tabolism提供了一个用于定量分析scRNA-seq中代谢通路活性的集成框架，可以与ssGSEA、Vision和AUCell等单细胞功能富集分析工具兼容。第二类方法，利用基于约束的数学模型重建代谢网络，代谢网络的重建通常基于KEGG和Reactome等代谢数据库。scFBA是第一个结合scRNA-seq数据和FBA来估计单细胞通量的计算工具，然后Compass和scFEA被相继提出。Compass基于Recon2对人体新陈代谢的重建，并通过线性规划解决基于约束的优化问题，对单个细胞中每个代谢反应的潜在活性进行评分。scFEA则通过概率模型来考虑通量平衡约束，通过神经网络来模拟通量变化和酶基因表达变化的非线性，并利用图形神经网络来解决优化问题，scFEA的分析结果可以进行多种具有生物学意义的下游分析，如细胞间代谢通讯等。
说在最后
scRNA-seq在生物医学研究中的应用促进了我们对疾病发病机制的理解，并为新的诊断和治疗策略提供了有价值的见解。随着包括临床样本在内的高通量scRNA-seq能力的扩大，对这些海量数据的深入分析，对进入该领域的非生信研究人员均是一个挑战。在篇综述中，作者回顾了典型的scRNA-seq数据分析的工作流程，包括原始数据处理和质量控制，几乎适用于所有scRNA-seq数据集的基本数据分析，以及针对特定科学问题的高级数据分析。在总结每个分析步骤的当前方法的同时，作者还提供了软件的在线存储库和打包的脚本来实现支持。对一些具体的分析任务和方法提出了建议和注意事项。作者希望这一资源将有助于非生信研究人员参与scRNA-seq数据分析，特别是对新兴的临床应用上有所帮助。
对于单细胞转录组学数据，我们能做的分析基本就是这些，根据自己课题的实际情况，我们可以选择其中某些分析策略，而并不一定全都要做。当然，单细胞水平的其他组学也在迅速发展，如ATAC-seq、单细胞DNA甲基化测序、Hi-C等技术，也有不少研究通过整合这些多组学数据，来更好的解释生物学问题。等到时机成熟，小编会再为大家总结其他的单细胞分析方法，今天的分享就到这了，祝大家都科研顺利~
[参考文献]
1、Su M, Pan T, Chen QZ, Zhou WW, Gong Y, Xu G, Yan HY, Li S, Shi QZ, Zhang Y, He X, Jiang CJ, Fan SC, Li X, Cairns MJ, Wang X, Li YS. Data analysis guidelines for single-cell RNA-seq in biomedical studies and clinical applications. Mil Med Res. 2022 Dec 2;9(1):68. doi: 10.1186/s40779-022-00434-8.
2、Zhao Ruohan, Bai Yicheng, Zhao Jingying, Hu Mei, Zhang Xinyan, Yang Min, Dou Tengfei & Jia Junjing (2023) Advanced analysis and applications of single-cell transcriptome sequencing, All Life, 16:1, 2199140, DOI: 10.1080/26895293.2023.2199140.
