NBIS系列单细胞转录组数据分析实战（八）：拟时序细胞轨迹推断

Posted 2023-04-19

参考技术A 在本节教程中，我们将学习如何通过拟时序分析推断细胞分化轨迹。 slingshot 包可以对单细胞RNA-seq数据进行细胞分化谱系构建和伪时间推断，它利用细胞聚类簇和空间降维信息，以无监督或半监督的方式学习细胞聚类群之间的关系，揭示细胞聚类簇之间的全局结构，并将该结构转换为由一维变量表示的平滑谱系，称之为“伪时间”。

运行slingshot的至少需要两个输入文件：即细胞在降维空间中的坐标矩阵和细胞聚类群的标签向量。通过这两个输入文件，我们可以：

这里，我们模拟了两个不同形式的数据集进行谱系推断分析。
接下来，我们构建第一个测试数据集，用于表示单个谱系，其中三分之一的基因与过渡相关。我们将该数据集包含在一个 SingleCellExperiment 对象中，并将其用于演示相应的分析流程。

接下来，我们构建第二个示例数据集，该数据集由降维空间坐标矩阵（如通过PCA, ICA和diffusion map等获得）和细胞聚类群标签（如K-means聚类等生成）组成。该数据集表示一个分叉的轨迹，可以用来演示slingshot包的其他功能。

要开始对单个谱系数据集进行分析，我们需要降低数据的维数，并过滤掉无信息的基因。这将大大提高下游分析的速度，同时将信息损失降至最低。

对于基因过滤步骤，我们保留了至少在10个细胞中表达，且基因表达count数至少大于3的基因。

由于我们使用的是模拟数据，因此无需担心批处理效应或其他潜在的混杂因素。在这里，我们将使用分位数归一化的方法对表达数据进行归一化处理，它可以使每个细胞具有相同的表达值分布。

slingshot基本假设是，转录相似的细胞在某些降维空间中会彼此靠近。由于我们在构造谱系和测量伪时间时使用欧几里得距离，因此对数据进行降维处理很重要。

在这里，我们将演示两种降维方法：主成分分析（PCA）和均匀流形逼近和投影（UMAP，通过uwot包调用）。

接下来，我们可以将降维的结果添加到 SingleCellExperiment 对象中。

slingshot的最终输入是细胞聚类簇的标签向量。如果未提供此参数，slingshot则将数据视为单个细胞簇并拟合标准主曲线。但是，我们建议即使在仅期望单一谱系的数据集中也应该对细胞进行聚类分群，因为它可以潜在地发现新的分支事件。

在此步骤中得到的聚类群集将用于确定潜在分化谱系的全局结构。这与聚类单细胞数据的典型目标不同，后者的主要目的是识别数据集中存在的所有生物学相关的细胞类型。例如，在确定全局谱系结构时，无需区分未成熟神经元和成熟神经元，因为两种细胞类型可能会沿着谱系的同一段落入。

在这里，我们采用了两种聚类方法：高斯混合建模和 k-means聚类，它们类似地假设低维空间中的欧几里得距离反映了细胞之间的生物学差。前者可以使用 mclust 包进行实现（Scrucca et al，2016年），它可以提供一种基于贝叶斯信息准则（BIC）的方法自动确定细胞聚类的簇数。

k-means聚类没有类似的功能，因此我们需要人为的指定聚类的细胞簇数。如果细胞簇数选择的太低，我们可能会错过一个真正的分支事件；如果它太高或有很多小簇，我们可能会开始看到虚假的分支事件。这里，我们选择聚类簇数为4进行k-means聚类分析。

现在，我们已经得到了slingshot分析所要用的两个文件：细胞降维的空间坐标矩阵和聚类分群的标签向量。接下来，我们就可以运行slinghsot进行伪时间推断分析。这是一个两步的过程，包括用基于细胞群集的最小生成树（MST）识别全局谱系结构，以及拟合主要曲线来描述每个谱系。

这两个步骤可以使用 getLineages 和 getCurves 函数分开运行，也可以使用 slingshot 封装函数一起运行（推荐）。这里，我们将使用封装函数对单轨迹数据集进行谱系分析，稍后将在分叉数据集上演示分步函数的用法。

为了使用PCA降维产生的空间坐标信息和高斯混合建模识别的细胞聚类簇标签来运行slingshot，我们将执行以下操作：

如上所述，如果未提供聚类分群信息，slingshot则假定所有细胞都是同一个聚类群，并且将构建一条拟合曲线。如果未提供降维空间坐标信息，slingshot将使用由返回的列表的第一个元素 reducedDims 中的信息。

输出的结果是一个包含了slingshot结果的SingleCellExperiment对象。大多数的输出以列表的形式添加到元数据中，并可以通过metadata(sce)$slingshot来进行访问。此外，所有推断的伪时间变量（每个谱系一个）都添加到colData中。如果要提取slingshot单个对象中的所有结果，我们可以使用 SlingshotDataSet 函数。这对于可视化很有用，因为SlingshotDataSet函数中包含了几种对象绘制方法 。下面，我们对单轨迹数据谱系推断的结果进行可视化处理，并用伪时间对细胞进行着色。

通过使用type参数，还可以看到基于细胞聚类群的最小生成树是如何推断出谱系结构的。

运行完 slingshot 推断出细胞分化谱系后，我们可以进一步鉴定出随伪时间或发育过程中表达变化的基因。我们将使用 tradeSeq 软件包演示这种类型的分析（Van den Berge等，2020）。

对于每个基因，我们将使用负二项式噪声分布来拟合通用加性模型（GAM），以对基因表达与伪时间之间的（可能是非线性的）关系进行建模。然后，我们使用 associationTest 函数来鉴定与伪时间存在显著关联的基因。

然后，我们可以根据p值挑选出变化最显著的基因，并通过热图可视化其在整个伪时间或发育过程中的表达情况。这里，我们使用top250基因进行展示。

接下来，我们将重点介绍 slingshot 包的一些其他功能。我们使用该包内置的 slingshotExample 数据集进行说明。此数据集中包含了细胞降维空间的低维表示，以及k-means聚类后的细胞簇标签。

首先，我们 getLineages 函数构建细胞分化谱系的全局结构。该函数使用 n × p 的细胞降维空间矩阵和聚类分群的标签向量 n 作为输入，通过最小生成树（MST）算法映射相邻细胞群集之间的连接，并进一步通过这些连接关系鉴定出细胞分化谱系的路径。该函数的输出结果是一个 SlingshotDataSet 对象，其中包含了输入文件，以及推断出的MST（由邻接矩阵表示）和分化谱系（细胞聚类簇的有序向量）。

通过指定已知的起始簇和终点簇，我们还可以以完全无监督的方式或以半监督的方式执行此分析。如果我们不指定起点，则 slingshot 根据简约性选择一个起点，以使分割前谱系之间共享的细胞簇数最大化。如果没有分割或多个簇产生相同的简约分数，则可以任意选择起始簇。

对于我们的模拟数据，这里选择以“Cluster1”作为起始簇。但是，我们通常建议根据先验知识（样本采集的时间或已知的基因标记）确定初始簇的类别。

在此步骤中，我们还可以指定已知的终点簇。当构造MST时，被指定为终端状态的细胞簇将被约束为仅具有一个连接（即，它们必须是叶节点）。该约束可能会影响树的其他部分的绘制方式，如以下示例中所示，在该示例中，我们将Cluster3指定为终点簇。

这种类型的监督对于确保结果与先前的生物学知识一致可能是有用的。具体地，它可以防止将已知的终末细胞命运归为暂时状态。

我们还可以传递一些其他参数，以实现对 getLineages 函数更好的控制：

为了构建平滑的谱系曲线，我们使用 getCurves 函数对 getLineages 的结果进行处理。该函数遵循类似于（Hastie and Stuetzle 1989）中提出的主要曲线的迭代过程。如果只有一个分化谱系，则生成的曲线只是经过数据中心的主曲线，并进行了一次调整：初始曲线是通过聚类中心之间的线性连接而不是数据的第一个主成分而构建的。这种调整增加了稳定性，并通常加快了算法的收敛速度。

当有两个或多个分化谱系时，我们向该算法添加一个额外的步骤：对共享细胞附近的曲线求平均。两种谱系应该在尚未分化的细胞上都很好地吻合，因此在每次迭代中，我们均对这些细胞附近的曲线进行平均。这增加了算法的稳定性，并产生了平滑的分支谱系。

getCurves 函数的结果也是一个 slingshotDataSet 对象，其中包含了分化谱系的主曲线和有关其拟合方式的其他信息。我们可以使用 pseudotime 函数提取不同分化谱系中细胞的伪时间值，其中的NA表示未分配给特定谱系的细胞的值。使用 curveWeights 函数可以提取相似的权重矩阵，将每个细胞分配给一个或多个谱系。

可以使用 curves 函数访问曲线对象，该函数将返回 principal_curve 对象的列表。这些对象中包含以下信息：

对于大型数据集，我们建议在运行 slingshot 函数时使用 approx_points 参数。通过设置 approx_points ，用户可以指定曲线的分辨率（唯一点的数量），这可以大大缩短计算时间，并且精度损失最小。我们建议该值设置为100-200。

10X单细胞转录组整合、转录组 && ATAC整合分析之VIPCCA

参考技术A

单细胞测序在基因调控、细胞分化和细胞多样性研究中具有革命性意义 。 随着近年来技术的显着改进，每个实验检测的单细胞数量呈指数级增长，同时大规模研究产生的数据集也在快速增长和积累。 因此， 当前单细胞研究中的一个主要计算挑战是对来自多个不同样本或跨不同平台和数据类型的测量进行标准化，以进行有效的综合和比较分析 。 这种综合分析需要开发单细胞数据对齐方法，该方法可以消除批次效应并考虑跨数据集的技术噪声。

最近开发了许多单细胞数据对齐方法 。它们中的大多数，除了一些值得注意的例外，例如最近的 iNMF ，都针对小型和中型数据集。这些现有的方法可以概括为四类：（i） 基于参考的方法 ，例如 Scmap-cluster 和 scAlign，它们基于注释良好的参考数据集对齐新的查询数据集； (ii) 基于聚类的方法 ，例如 Harmony 、DESC，它们通过迭代优化聚类目标函数来消除批效应并在嵌入空间中对齐样本； (iii) 基于匹配的方法 ，例如 MNN 和 Scanorama ，它们应用相互最近的邻居策略来识别跨数据集的重叠单元格和 (iv) 基于投影的方法 ，使用统计模型将来自不同数据集的单个细胞投影到较低的维空间，包括对投影应用典型相关分析的 Seurat ，使用来自非负矩阵分解的潜在因子进行投影的 LIGER ， and scVI and others that use variational techniques for projection.

然而， 大多数现有的对齐方法都存在固有缺陷，无法成功应用于大型数据集 。具体而言， 基于参考的方法的对齐将受到参考数据大小和参考中可用的预选细胞类型注释的限制，因此当数据大小增加时，可能会导致错过新发现的机会增加 。像 MNN 这样的基于匹配的方法使用往返游走策略，该策略需要为具有两个以上样本的数据集生成所有成对对齐，这对于大样本量来说将是耗时的。具有复杂参数模型的方法（例如 LIGER 和 scAlign）或具有复杂事后数据处理的方法（例如 Seurat ）也难以扩展到大型数据集。 基于 ZINB 的方法（例如 scVI）在捕获多个数据集的复杂表达特征方面可能效率较低 。尽管一些现有的最新方法可以扩展到大型数据集，但由于复杂的参数模型，它们仍然有可能不准确地对齐细胞。因此， 迫切需要开发在计算上也有效的有效对齐方法 。

除了迫切需要开发可扩展的比对方法外，当前比对方法的另一个阻碍问题是它们的性能通常仅使用单细胞 RNA 测序 (scRNA-seq) 数据进行基准测试和优化。 因此，大多数现有的比对方法不适合整合其他单细胞测序数据类型，例如使用测序 (scATAC-seq) 进行转座酶可及染色质的单细胞测定。 此外， 现有的比对方法（如 Seurat）返回的结果只能保留真实的细胞间关系（或相似性），而不能代表基因表达水平，不适合进行差异表达分析或富集分析等下游分析 。

为了应对这些挑战， 作者提出了一个统一的计算框架 VIPCCA，它基于非线性概率典型相关分析，用于有效且可扩展的单细胞数据对齐 。 VIPCCA 利用来自深度神经网络的尖端技术对单细胞数据进行非线性建模，从而允许用户通过跨技术、数据类型、条件和模式的多个单细胞数据集的集成来捕获复杂的生物结构。此外，VIPCCA 依靠 变分推理 来进行可扩展计算，从而能够将大规模单细胞数据集与数百万个细胞有效集成。重要的是，VIPCCA 可以将多模态转换为低维空间，而无需任何事后数据处理，这是与现有对齐方法形成直接对比的独特且理想的功能。

加载

Loading data

该函数仅适用于 fit_integrate() 函数训练生成的 AnnData。 在基因表达矩阵中随机选择 2000 个位置。 x轴代表这些位置原始数据的表达值，y轴代表同一位置的vipcca整合后数据的表达值。

After filter, we converting Seurat Object to AnnData via h5Seurat using R packages. In this case, the atac.h5ad file will be generated in the corresponding path .

生活很好，有你更好