这篇综述今年6月1日发表在Journal of the American Society of Nephrology[IF=9.274]上。
在过去的5年里,单细胞方法使得能在单个实验中监测数千个单细胞的基因和蛋白质表达、遗传和表观遗传变化。随着测量方法的改进和测序成本的降低,这些数据集的规模正在快速增加。主要的瓶颈仍然是对单细胞实验产生丰富信息的分析。在本综述中,研究者对目前在这个领域广泛使用的分析流程进行了简要概述(Figure 1)。
1. 获得数据矩阵和质控
单细胞分析的一个关键技术进展是条形码(barcoding)的发展,它在允许大规模测序的同时将成本最小化。条码在逆转录时被加入RNA分子,用于识别单个细胞和特异分子。单细胞分析的第一个步骤是得到数据矩阵,它代表原始测序文件中转录数据库的条形码(细胞)。对于10X Genomics的数据来说,CellRanger是最常用的pipeline。它包括了对基因组的测序读长的拆分和比对、注释和定量。其他可选的技术包括UMI tools、zUMIs、kallisto、STAR和STARsolo等。
每个barcode代表一个单细胞,一个双细胞或不包含细胞只包含环境RNA的空液滴。一个重要的问题是标准的流程是将测序数据和转录组(如处理过的成熟mRNA)进行比对。然而,单核RNA数据或表观基因组数据(ATAC-seq)应该与全基因组比对,因为细胞核主要含的是pre-mRNA,其序列包括内含子区域。原始读长计数通常还会过滤掉极少数细胞中检测到的基因,从而有效地降低数据矩阵的大小。
分析流程的下一步是质量控制(QC),比如识别每个barcode的count数,每个barcode的基因数和每个barcode中的线粒体基因百分比。
基因数偏低和线粒体读数百分比较高通常提示细胞质量不佳。然而,一些细胞比如肾近曲小管和远曲小管细胞也富含线粒体。异常的高reads数和基因数常常代表双细胞。对于这种情况可以使用一些双细胞检测工具比如DoubletDecon、Scrublet和DoubletFinder等。
然而,双细胞鉴定中一个重要的问题是处于从一种细胞向另一种细胞分化过程中的包含两种marekr基因的细胞可能会被标记为双细胞,这有时会导致检测的假阳性。此外,这些工具对于同型细胞所组成的双细胞检测效果往往不佳。
控制环境RNA污染也是很重要的。环境RNA是指存在于单细胞悬液并在封装的时候被整合进油包水液滴中的RNA。我们通常使用 SoupX估计空液滴中的环境 RNA 污染。另外一个替代性的包是CellBender,它从基于UMI的单细胞表达矩阵(raw)中移除由于环境RNA分子和随机条形码交换造成的计数。
2. 标准化
单细胞数据需要不同类型和水平的标准化。一个常用的方法是假设每个细胞都有相同的初始转录本数,只需将数据标准化为每百万计数。Scran、Seurat、SCtransform、SCnorm、BayNorm等都提供了数据标准化方法。在标准化之后,数据进行了 log(x11)转换。在分析时通常需要通过regress out掉细胞周期相关变异的影响,一些pipeline如Seurat和Scanpy就将这一步整合进了标准分析。当然这些平台也可以通过回归消除掉其他技术或生物学因素的影响。
3. 批次校正和数据整合
大多数情况下,会生成多个数据集,因此需要额外的批次校正和数据整合方法(图 3)。
包含多个不同实验和不同方法的数据集通常使用非线性方法来整合。在Seurat中有一个基于参考的整合选项,它使用的是典型相关分析或互反主成分分析。Scanorama是Scanpy中使用的另一种普遍且高效的方法。最近,Harmony变得越来越流行,快速成为了单细胞数据整合最常用的方法。首先,PCA衍生的嵌入矩阵和批次元数据用于缩放,以便为每个细胞单元赋予一个长度参数。然后,用常规的k-均值聚类对缩放后的数据进行聚类中心初始化。最后,通过迭代地将特定批次的中心拉到聚类中心,来消除批次效应。见Harmony原理介绍和官网教程
4. 可视化和聚类
可视化的第一步是特征选择(feature selection),也就是保留提供有用信息的基因(1000-5000)并过滤掉其他基因,这在Seurat (Figure 4) 和Scanpy(Figure 5)中均可实现。可视化是将数据整合进低维以方便观测。通常,降维是通过线性和非线性方法来实现的,PCA是聚类和轨迹分析的基础。它是一种线性变换,在整个PCA中保留了细胞之间的欧氏距离。在常见的Seurat流程中,PCA被用于预处理阶段,可以将主成分映射到技术和生物协变量中,以了解其效能(图4)。使用jackstraw图和elbow图(Figure4 A B)有助于主成分数的挑选。单细胞数据可视化主要使用其他非线性降维方法,如t-SNE。这种方法侧重于以舍弃全局结构为代价来捕捉局部相似性。UMAP方法也因其运行速度快而广受推崇,它似乎能更好地捕捉潜在的数据结构,并能在两个以上的维度上汇总数据,因此,它现在最常用于单细胞数据可视化。UMAP和t-SNE的关键缺陷是它们高度依赖于用户选择的参数,而且结果对这些参数高度敏感。最值得注意的是这些可视化降维方法都没有保留细胞之间的距离,因此不能直接用于后续的下游分析。
基于基因表达相似性而形成的细胞群(cluster)是分析的第一个直接结果。细胞聚类允许基于基因表达的相似性对细胞进行分组来推断细胞类型。聚类是一种基于距离矩阵的无监督机器学习过程,默认的聚类方法是在单细胞KNN图上实现的Louvain社区发现算法(Louvain算法)。(注:Louvain算法已经成为Scanpy和Seurat单细胞分析平台中默认聚类的方法。)细胞在图中以节点的表示。每个细胞连接着它的K最近邻细胞,细胞之间的距离是基于PC降维表达空间计算的欧几里得距离。一个重要的问题是Louvain聚类的分辨率由使用者决定,分辨率决定了细胞群的数目。我们推荐进行次分群(subclustering),也就是对原有数据集聚类后再挑选感兴趣的细胞群重新分群。也有助于更精细的细胞型的划分。被整合进Leidenbase包中的Leiden社区检测算法是Louvain算法的另一替代性选择,被用于Monocle轨迹分析的默认算法。新的聚类方法使用神经网络和人工智能正在被开发。
值得特别注意的是,cluster并不一定意味着细胞类型。用户定义的resolution参数可以决定cluster的多少。分群和注释是高度相关的,通常需要反复进行,较为费时。目前,关于最佳聚类参数还没有统一的标准,因此,可以接受对同一数据进行多种版本的聚类和解释。Wilcoxon秩和检验用于通过组间表达差异对基因进行排序。
经典的细胞注释使用了真实正确的外部的数据集,可用的外部数据集包括Susztak Lab,Humphreys Lab,Tabula Muris,Human Cell Atlas,ImmGen Consortium,Garnett,SingleR,CHETAH和MOANA等。对同一细胞类型,不同数据集提供的marker基因可能不同,需要综合判断来进行注释。
5. 细胞水平的分析:细胞组分变化、分解和轨迹分析
细胞组分的变化(数据集中各种细胞类型的比例)与疾病状态有很强的关联性,这是单细胞分析最简单的结果之一。这些数字可提供不同条件下的相对估计,但是由于单细胞文库制备时细胞捕获偏差等因素,从单细胞数据得到的细胞成分变化可能不准确。此外,同一器官不同部位的样本细胞组分可能也存在差别。为了推断bulk RNA-seq数据的细胞类型组成,最近开发的MuSiC可以以单细胞表达数据为参考,对组织细胞类型进行去卷积分析。MuSiC使用加权非负最小二乘回归来估计细胞类型组分,可选的替代方法包括 CIBERSORT、BSEQ-sc和BisqueRNA。
仅仅通过离散分类系统如聚类并不足以充分描述细胞多样性。轨迹分析捕捉了细胞在转化过程中的显著特征。驱动这种可观测到的异质性的生理过程是连续性的。因此,捕捉细胞特性之间的转变、分支分化过程或生物功能的渐进的非同步化的变化需要基因表达的动态模型。Monocle 是一种机器学习方法,可以重建每个细胞在从一种状态转换到另一种状态时必须执行的基因表达变化过程。它基于反向图嵌入(reverse-graph embedding),一种高度可扩展的非线性流形学习技术。在这个方法学习了分化轨迹之后,它可以将每个细胞放置在分化轨迹的正确位置,这被称为拟时序列(pseudotime)。另一种更新的观测细胞轨迹的方法是使用velocyto包的RNA velocity(RNA速率)分析。RNA速率是基因表达状态的时间导数,可以通过区分剪接和未剪接的mRNA来直接估计,它是一个高维向量,可在数小时的时间尺度内预测单个细胞的未来状态。TradeSeq是基于Slingshot的一种新方法,其性能优于其他简单轨迹分析方法。另一个有用的包是PHATE,这是一种使用数据点之间的信息几何距离来捕获局部和全局非线性结构的可视化方法。推断的轨迹不一定代表生物过程,其他来源的信息将有助于轨迹分析结果的解读(Figure 6)。
6. 基因水平分析:差异表达、基因调控网络、驱动通路和细胞互作
差异表达(DE)分析输入的是未经校正的数据集。Seurat可以使用不同的模型进行差异表达分析(Figure 6)。MAST使用了一种栅栏模型(hurdle model)来计数drop-out。为了将scRNA-seq数据集信息与其他表型变量相关联,基于回归的模型可以结合多个样本及其相关的表型特征,以将某些细胞类型(如近端肾小管细胞)中的基因表达变化与相应的定量测量表型(例如GFR、蛋白尿)相关联。尽管DE分析工具通常允许用户数用灵活的整合混杂因素,用户必须警觉哪些变量可以被加入模型。基因水平分析也可以与基因集富集分析方法相结合,如基因集富集分析或加权相关网络分析。
为了解读DE的结果,我们通常基于基因参与的常见生物学进程将基因分组。生物学进程标签储存在MSigDB、GO、KEGG和Reactome等数据库中。
单细胞分析领域的最新进展是使用成对基因标签进行配体-受体分析,细胞群之间的相互作用是通过受体及其同源配体的表达来推断的。配体-受体对标签可以从相应的数据库,如CellPhoneDB或Connectome获得,并使用统计模型来解释跨群的高表达基因。
7. 单细胞分辨率下的基因调控
snATAC-seq(但细胞核ATAC)允许通过分析染色质可及性来分析单细胞中的表观基因组景观(Figure 7)。在多种snATAC-seq分析的工具中,比较好的包括Ren Lab开发的SnapATAC,Satija Lab开发的Signac和Greenleaf Lab开发的ArchR。我们更推荐SnapATAC,这是一种非线性降维方法。在CellRanger中生成barcode-细胞矩阵后,使用SnapATAC对矩阵进行预处理。为了识别异质性组织中的细胞类型,SnapATAC利用扩展映射中获得的低维嵌入来消除批次效应。然后,使用从k最近邻算法中选择的k值作为输入,使用Louvain算法执行聚类。为了识别不同细胞类型中的富集基序,可以使用 HOMER或chromVAR对snATAC-seq数据进行转录因子分析。为了研究开放染色质变化如何与细胞分化和细胞命运决定相关联,使用Monocle3进行轨迹分析,通过使用潜在语义索引减少维度,并通过UMAP进行可视化。为了了解开放染色质和靶基因表达的变化,主要通过分析在Cicero中实现的两个峰的共同可及性来进行峰-峰相关性研究。值得一提的是,Satija实验室开发了数据分析自动化平台Azimuth,它允许将单个数据集映射到参考数据集中。研究人员可以上传本地生成的数据集,该软件包会自动执行上述步骤,并与人类血液参考数据进行聚类。