按树维度推断轨迹存在和子树覆盖的子集特异性-厦门畜牧期刊杂志论文发表
洛夫莫尔·特尼亚，宋明州
出版日期： 2022年02月08日
抽象
细胞分化等生物过程的复杂性反映在细胞状态之间的动态转变中。轨迹推理使用单细胞生物学推动的方法将状态排列成一个过程。然而，目前的方法（所有返回最佳轨迹）都无法充分评估噪声模式的统计显著性，从而导致推断轨迹的不确定性。我们通过欧几里得最小生成树的维度度量、测试统计量和空分布，引入了多变量数据中轨迹存在的树维度检验。树维度量值可按树大小线性时间计算，比起对次级树枝无动于衷的全局不敏感的叶子数或树径，可以更有效地总结分枝的程度。检验统计量量化轨迹存在，其零分布是在数据中没有轨迹的原假设下估计的。在模拟和真实的单细胞数据集上，该测试优于直观的叶子数量和树直径统计数据。接下来，我们根据子集在最小生成树中子集的最小子树覆盖率，为子集动力学的组织特异性开发了一个度量。我们发现，在人类和小鼠发育中，途径基因表达动力学的组织特异性是保守的：包括钙和Wnt信号传导在内的几种信号转导途径最具组织特异性，而遗传信息处理途径如核糖体和错配修复则最不重要。树维度测试和子集特异性度量都没有任何要调整的用户参数。我们的工作打开了一扇窗，可以优先考虑发育和其他多变量动力系统中的细胞动力学和途径。
作者简介
现代生物学现在在发育过程中经常研究转录组图谱。这种做法需要计算方法来量化细胞状态的动态变化及其异质性。许多方法处理单细胞转录组数据以重建细胞轨迹，这是细胞从早期发展到晚期发育阶段的顺序。由于转录组数据中的噪声，因此非常需要量化观察到的数据由于偶然性而呈现轨迹状模式的可能性。为了满足这一需求，我们开发了一个树维测试，以量化基于图论概念的多变量数据中轨迹存在的证据。通过此测试，人们可能会由于数据质量低而拒绝轨迹存在，或者接受具有高统计显著性的轨迹。现在，人们可以根据生物途径的轨迹质量对它们进行排名。我们还引入了一个子集特异性测量，以量化细胞或途径动力学如何具有组织特异性。我们发现，通路组织特异性在人和小鼠之间高度保守。轨迹存在测试和子集特异性为研究发育生物学提供了一套独特的信息学工具集。
数字
Fig 10Fig 11Fig 1Fig 2Fig 3Fig 4Fig 5Fig 6Fig 7Fig 8Fig 9Fig 10Fig 11Fig 1Fig 2Fig 3
引文： Tenha L，Song M（2022）通过树维度推断轨迹存在和子树覆盖的子集特异性。PLoS Comput Biol 18（2）：e1009829。https://doi.org/10.1371/journal.pcbi.1009829
编辑 器： Lilia M. Iakoucheva，加州大学圣地亚哥分校，美国
收到： 七月 29， 2021;接受： 一月 12， 2022;发表： 二月 8， 2022
版权所有： ? 2022天岚，歌曲。这是一篇根据知识共享署名许可条款分发的开放获取文章，该许可允许在任何媒体上不受限制地使用，分发和复制，前提是注明原始作者和来源。
数据可用性： 我们已经在名为"TreeDimensionTest"的R软件包中实现了用于轨迹存在的TDT以及最小子树覆盖算法，该软件包可在 https://www.cs.nmsu.edu/~joemsong/TrajTest/TenhaSong2022TDT/trajectoryTest.zip https://CRAN.R-project.org/package=TreeDimensionTest 提供有关数据的其他源代码和信息。
资金： MS从美国国家科学基金会获得赠款1661331，https://www.nsf.gov MS从美国农业部获得2016-51181-25408 https://nifa.usda.gov 资助者在研究设计，数据收集和分析，出版决定或准备手稿方面没有任何作用。
相互竞争的利益： 作者宣布不存在相互竞争的利益。
介绍
识别细胞状态之间的动态转变可以更深入地了解生物系统内的发育，疾病过程或环境反应。单细胞生物学促进了对这种动力学的探索。轨迹推理利用细胞分辨率下的组学数据来识别细胞状态进展，为此已经开发了许多计算方法[1，2，3]。广泛的评估表明，许多方法可以在不使用先前的生物学信息的情况下，以高度的准确性推断某些类型的轨迹[4]。轨迹推理方法通常在嵌入在高维空间中的低维流形上运行，采用各种策略来捕获轨迹。TSCAN等方法利用构建在聚类质心上的最小生成树（MST）来捕获数据底层的轨迹结构[5]。SLICER使用局部线性嵌入和k-最近邻图来查找轨迹[6]。拓扑数据分析发现复杂高维数据的紧凑表示[7，8]。Vandaele等人开发了一种推断图数据中拓扑结构的方法，适用于轨迹推理[9]。他们强调了一些挑战[10]：大多数方法倾向于低估轨迹图表示中的叶子数量。他们还表明，拓扑信息与连续像元轨迹推理算法的性能相关，并且许多具有轨迹的数据集缺乏足够的拓扑信息来进行有效推理。
然而，据我们所知，现有方法大多假设数据中存在轨迹，因此无论统计显著性如何，总是推断轨迹。由于技术或生物学原因，即使预期到轨迹模式，实验也可能无法捕获轨迹模式。最好的轨迹可能只捕捉到偶然产生的随机变化。此外，没有可用的方法来优先考虑亚空间中由生物途径上的基因跨越的轨迹存在，其动力学可能与涉及所有基因的细胞轨迹截然不同。最后，没有已知的统计数据用于测量细胞或途径动力学的组织特异性。这种局限性严重阻碍了我们解释隐藏在多变量生物数据中的突出信号的能力。
为了解决轨迹存在问题，我们建立了一个树维测试（TDT），该测试利用图理论统计来表征指示观察到的数据中动态过程的存在和特异性的模式。我们的方法前提是，多变量数据中轨迹模式的存在可以通过原始数据的相应欧几里得MST（EMST）中的线性度和分支度来测量。我们引入了一个基于树维度量 T 的统计框架d的 EMST 来完成任务。它计算从树维度量值 T 派生的检验统计量 Sd以量化轨迹存在的统计证据，使用球面多元正态随机向量总体上对数正态的零分布。
为了在诸如组织类型之类的上下文中研究轨迹动力学的特异性，我们引入了一个子集特异性度量，该度量基于给定子集的最小子树覆盖率。组织特异性测量不同发育组织类型中途径的差异表达，从而实现途径行为表征。
我们在模拟和真实的单细胞数据上评估了我们的方法，其中TDT在识别数据中是否存在轨迹模式方面的表现大大优于树直径和叶子数量。我们应用我们的方法测试了发育中的哺乳动物组织类型的转录组数据中通路轨迹的存在，然后量化了通路动力学的组织特异性。事实上，在发育过程中，整体细胞轨迹与特定的路径轨迹在质量上是不同的。几种信号转导途径显示出高组织特异性，而在哺乳动物发育过程中观察到遗传信息处理途径的组织特异性较低。
尽管在多变量数据的图和拓扑表征方面做了大量工作，但我们的方法为动态模式的分析增加了独特而严格的统计视角。由于我们的方法解决了关于轨迹存在的更一般问题，因此不需要识别特定的轨迹，因此可以在轨迹查找方法的上游使用。我们的生物学发现支持树维框架的有用性，该框架适用于由于现代生物技术创新稳步提高时间和细胞分辨率而具有越来越复杂动力学的分子生物学数据。
方法
我们为多变量数据中的轨迹存在设计了树维测试（TDT）。该测试有四个主要步骤，如图 1 所示。第一步以紧凑的形式查找输入数据的EMST，同时保留基础信号（如果有）。数据中的样本被视为完整图形的顶点，每对顶点之间都有边。边的长度是其两个顶点之间的欧氏距离。从完整的图形中，EMST被计算为数据的紧凑表示形式。第二步利用 EMST 表示来测试数据中是否存在轨迹模式。目标是识别EMST的特征，以指示底层模式是动态形成轨迹还是只是随机的。因此，我们开发了树维度量度Td将 EMST 结构映射到轨迹模式的存在与否。Td总结了 EMST 的分支程度，对于高度分支的 EMST 而言，它很高，而 T 则较低d提示一个强大的轨迹模式。通过第三步，基于 T 计算检验统计量 Sd以集成统计支持。在最后一步中，我们从球面多元正态随机数据中推导出 EMST 零总体上 S 的对数正态零分布。利用空分布的上尾概率，我们最终计算出S的统计显著性，以确定轨迹的存在。
thumbnail 下载：-厦门畜牧期刊杂志论文发表
个人电脑幻灯片
巴新放大图片
断续器原始图像
图 1. 轨迹存在的树尺寸检验概述。
输入是多变量数据。在数据点上计算欧几里得最小生成树 （EMST）。轨迹推断检验统计量 S 是从树维度量值 T 计算得出的d的EMST。S 的对数正态零分布是从零总体派生而来的，该分布遵循没有轨迹模式的球面多元正态分布。计算观测到的统计量 S 的 p 值以量化轨迹模式存在的统计显著性。
https://doi.org/10.1371/journal.pcbi.1009829.g001
另外两个图论统计对于表征树的线性程度是直观的。首先是树中叶子L的数量。叶子是树中只有一个入射边的顶点。L 测量树中分支事件的频率。如果树是完全线性的，则 L 最小化，也称为路径图。在光谱的另一端，L在星树上最大化。然而，L严重受噪声的影响，倾向于降低全局强轨迹模式，该模式可能在局部具有嘈杂的小分支。第二个是树直径Dm，树中任意一对顶点之间的最大路径长度。Dm如果一棵树是线性的，则最大化，并且在星形树上最小化，星形树是最大分支的。Dm，可有效识别树的主干，对可以表示数据收集中的多分支事件的辅助分支不敏感。
树的维度及其度量
现在，我们引入一个树维度测量值，以结合树直径 D 的强度m和叶子的数量L，以稳健地表征树木的分枝程度。其基本原理是捕获全局轨迹模式，同时对树中的局部噪声分支不敏感。
定义 1（树维度）。我们将树的第一维定义为其最长的路径之一。树的 （k + 1）-第 （k > 0） 个维度是两片叶子之间的 （k + 1）-第 1 条最长路径，而不是位于维度 1 到 k 上。如果只有一个这样的叶子可用，则维度 k + 1 由将唯一叶子连接到前一维度中的顶点的最短路径跨越。如果没有这样的叶子可用，则树正好具有 k 个维度。但是，任何维度的中间顶点都可能位于前一个维度上。
很明显，树的每个维度的长度都可以不同，由一加上该维度上的顶点数来定义，但在较低维度上则不然。
定义 2（树维度度量）。树维度度量 Td，量化树中分支的程度，是每个维度的长度之和，由长度归一化（Dm+ 1）的第一个维度。该定义很直观，但计算起来不方便。不难看出，Td等效地由树中顶点 N 的数量及其直径 D 确定m和叶子的数量 L 由
(1)
对于 N = 1 的单例树的退化情况，我们将叶子数量 L 定义为 1，树直径 Dm为零。这组 Td= 1 表示单例树。
图2A显示了第一维的路径图。图2B和2C都有两个维度。图2D中的树有三个维度。它们对应的树维度量 Td与尺寸数相同，除了在图2B Td为 1.5，因为两个维度具有不同的长度。
thumbnail 下载：
个人电脑幻灯片
巴新放大图片
断续器原始图像
图 2. 树、维度和维度度量。
每个树维度都以不同的颜色突出显示。（A） 一维树或路径图，具有树维度量值 Td= 1。（B） 具有 T 的二维树d= 1.5。（C） 具有 T 的二维树d= 2。（D） 具有 T 的三维树d= 3。
https://doi.org/10.1371/journal.pcbi.1009829.g002
算法 1 三维测量-Measure（X） 计算树维度量 Td在多变量输入数据 X 上。在 X 上首次获得 EMST T 后，它计算树直径 Dm.我们从树中的任何顶点 v 开始对 T 应用广度优先搜索 （BFS），同时跟踪每个顶点与 v 的距离。在 BFS 完成后，我们选择与 v 距离最长的顶点 w。我们再次运行BFS算法，从顶点w开始，同时跟踪每个顶点与w的距离。BFS 完成后，顶点 z 与 w 的距离最长，z 到 w 的距离为 Dm.我们通过度数为 1 的顶点数得到叶子的数量 L。Dm然后使用 L 来计算 Td.
算法 1Tree-Dimension-Measure（X)
输入：多变量数据 X
1 T = 在 X 上查找欧几里得 MST
2 // 从树中的某个顶点 v 查找最远的顶点 w
按广度优先搜索：
w = Breadth-First-Search（T， v)
3 // 在树中查找从顶点 w 最远的顶点 z
按广度优先搜索：
z = Breadth-First-Search（T， w)
4 树直径 Dm= 从 w 到 z 的路径长度
5 L = 长颈鹿的铟（T)
6 计算树维度度量 Td通过均衡器 （1）-厦门畜牧期刊杂志论文发表
7 返回 Td
树维度量的数学属性
命题 1.树维度量 Td 当且仅当树是路径/线性图时，最小化为 1。
证明。我们将一个顶点（N = 1）的树视为路径图，它始终具有Td= 1（根据定义）。N ≥ 2 个顶点的路径图正好有两个叶子 （L = 2），所有其他顶点的度数为 2。树直径 Dm因此为 N ? 1。因此，Td= 1 乘以均衡器 （1）。
非路径图的 L > 2。它必须遵循?L/2 ? 1? > 0。作为直径 Dm小于任何树中的顶点 N 个数，N ? Dm? 1 ≥ 0.因此，我们到达 N ? Dm? 1 + ?L/2 ? 1? > 0，表示 Td> 1 by Eq （1）。
因此，Td当且仅当树是路径图时，才会最小化为 1。
命题 2.树维度量 Td 当且仅当我们有一个星形树，在所有N≥3个顶点的树中，最大化。
证明。N ≥ 3 个顶点的星树具有 N ? 1 个叶子和一个 N ? 1 度的顶点。因此，树直径 Dm为 2。因此。
根据定义，N ≥ 3 的非星形树的直径必须为 x > 2。设 y 为叶子的数量。由于任何 N ≥ 3 个顶点的树最多有 N ? 1 个叶子，因此我们有 .树直径 x > 2 表示 x + 1 > 3。因此，通过 Eq （1） 可以看出 .
因此，Td当且仅当我们有一个星形树，在所有 N ≥ 3 个顶点的树中，才最大化。
树木尺寸测量、直径和叶子数量的示例
树维度量 Td利用树直径 D 的优势m同时通过合并叶子L的数量来减轻其缺点。图3说明了T的根本区别d从而获得比 D 更好的性能m和 L.图3A显示了一个具有T的多furcating树d= 2， Dm= 6 和 L = 4。图3B显示了一棵有T的树d= 2.28，Dm= 6 且 L =8。两棵树具有相同的 Dm.图3B中的树有更短的树枝，但只是膨胀了它的Td值稍轻。Dm对于两棵树都是一样的，对结构差异不敏感。这凸显了D的不敏感。m到二次分支，而 Td捕获全局和局部树属性。图3B中树上的L值急剧加倍至8，突出显示了L可以是局部极值的情况。T的优越性能d在 D 上mL在结果中进一步证明。
thumbnail 下载：
个人电脑幻灯片
巴新放大图片
断续器原始图像
图 3. 两个具有相同顶点数的树及其树维度度量 Td，直径 Dm和叶子数量 L.
顶点的颜色表示唯一的树尺寸。（A）树是二维的，Td=2， Dm=6 且 L=4。（B）树是四维的，Td=2.28， Dm=6 和 L=8。
https://doi.org/10.1371/journal.pcbi.1009829.g003
轨迹存在的树尺寸检验
现在，我们描述了树维数检验 （TDT），这是一种针对轨迹存在的统计检验。树维度量 Td是轨迹模式的有用数学描述符，尽管在统计上并不完全授权。因此，我们引入轨迹检验统计量 S：
(2)
其中 N 是树中的顶点数，以及 T 的最小值和最大值d由命题 1 和 2 给出。S 与 T 呈负相关d.强轨迹模式将具有较大的 S 值。最重要的是，S 随样本数量 N 的增加，以促进未反映在 T 中的统计支持d.
我们还通过 S/N 定义检验效应大小，S/N 是树维度度量 T 的对数d负和线性缩放到[0，1]，0表示没有轨迹模式，1表示完美的轨迹模式。
为了测试轨迹存在，我们定义了原假设，即数据中不存在轨迹模式。我们选择零总体来遵循球面多元正态分布（MVN），因为它呈现了一个与均值衰减的各向同性点云，代表了人们不会期望的轨迹模式。
接下来，我们获取检验统计量 S 的空分布。给定一个具有 N 个样本和 d 个维度的输入数据集，我们从具有单位协方差矩阵的球面 MVN 分布中获得 N 个 d 维随机向量样本。我们在这些样本上计算一个S。多次重复该过程，我们得到 S 的经验零分布（Alg. 2 Null-Distri distributionion）。
实验结果表明，S的零分布可以通过基于贝叶斯信息准则（BIC）和赤池信息准则（AIC）的对数正态分布来最好地近似模型选择。我们还执行了 Kolmogorov-Smirnov （KS） 检验，以确定 S 统计量零分布与三个候选分布的相等性：对数正态分布、伽马分布和正态分布。每个都适合从八维球面MVN分布的1000个数据样本中获得的S的经验零分布。通过最大似然估计获得每个分布的两个参数。
对数正态分布，达到最低（最佳）BIC和AIC分数（图4），以及KS检验p值与数据最小显着偏差，被认为是S零分布的最佳选择。使用空分布的上尾概率，我们可以计算观测到的检验统计量 S（Alg. 3 Test-Trajectory-Presence）的统计显著性（p 值）。因此，可以在给定的I型错误率（例如，0.05）下决定轨迹存在。
thumbnail 下载：
个人电脑幻灯片
巴新放大图片
断续器原始图像
图 4. 近似轨迹检验统计量 S 的经验零分布。
BIC、AIC 和 KS 检验 p 值 （A） 对数正态 （lnorm）、（B） gamma 和 （C） 正态（范数）分布（红色曲线）与模拟的零检验统计量 S 值（直方图）拟合后。
https://doi.org/10.1371/journal.pcbi.1009829.g004
算法 2Null-distribution（X， B)
输入：X是d维中N个样本的多变量数据;
B 是空总体的随机样本数
1 for i = 1 到 B
2 // 获取 n 个 d 维球面 MVN 随机样本
蒙特卡罗采样向量：
3 倍π= 球面-MVN（N， d)
4 吨d= Tree-Dimension-Measure（Xπ)
5 秒0[i] = Tree-Dimension-Statistic（Td)-厦门畜牧期刊杂志论文发表
6 通过最大似然估计拟合对数正态零分布
?
7 返回对数正态零概率密度函数
轨迹存在的树维数检验的时间复杂度
现在，我们分析树维度算法的运行时。给定具有 d 维和样本数量 N 的多变量数据 X，Alg. 1 Tree 维使用在 "mlpack" 软件包 [12] 中实现的双树 Boruvka 算法 [11] 计算大约 O（N log N） 中的 EMST。对于带有 B 模拟的 Alg. 2 Null-Distribution，运行时为 O（BN log N）。因此，Alg. 3 Test-Trajectory-Presence 的整体运行时复杂度约为 O（BN log N）。
通过最小子树覆盖测量组织对基因表达动力学的特异性
给定来自多种组织类型的观察到的基因表达数据，我们检查来自一种组织类型的数据点在基因表达中是否彼此接近或与其他组织类型混合。这允许人们评估组织特异性。如果基因组中的所有基因都包含在表达数据中，则相应的组织特异性是关于涉及组织类型所有基因的细胞状态;如果只检查通路中的基因，则组织特异性是关于组织类型通路的动力学。在这里，我们使用给定基因表达数据的EMST T研究组织特异性。
算法 3Test-Trajectory-Presence（X)
输入：多变量数据 X
1 空分布 = 无分布（X， B)
2 树维度量 Td= Tree-Dimension-Measure（X)
3 S = Tree-Dimension-Statistic（Td)
4
5 返回 p 值
现在，我们使用顶点集 V 在 T 上定义一个组织特异性度量，其中每个顶点都用组织类型标记。让Wt是 V 中由组织类型 t 标记的顶点子集。关键步骤是找到一个子树 Tt的 T 使得 Tt封面 Wt具有最小数量的顶点。让 Vt是 T 的顶点集t，其中可能包含未由 t 标记的顶点。
我们通过以下方式定义树 T 对组织 t 的特异性：
(3)
其中集合的绝对值是其基数。如果最小子树覆盖的所有顶点都由相同的组织类型 t 标记，则组织 t 的组织特异性为 1。当最小子树覆盖包含许多不属于组织类型 t 的顶点时，组织特异性降低。我们用M种组织类型的平均组织特异性来表示树的组织特异性：
(4)
标识覆盖顶点集 V 的步骤t对于顶点子集 Wt，我们引入算法 4 Minimum-Subtree-Cover，它获得最小子树 Tc顶点集 Tc和边缘集 Ec覆盖树 T 中的顶点子集 W。算法 5 Cover-DFS 使用 DFS 遍历查找要包含在最小子树覆盖中的所有顶点。Cover-DFS 要探索的第一个顶点必须位于顶点子集 W 中。在 DFS 遍历期间，如果顶点 u 位于最小子树覆盖中，则其发现顶点 v 必须包含在覆盖中。
算法 4Minimum-Subtree-Cover（tree T， 顶点子集 W)
输入：T，一棵树（V，E），顶点集V和边集E
W，顶点集 V 的子集，用于查找最小子树覆盖
1 设 w 为子集 W 中的顶点
2 Cover-DFS（T， W， w） // 识别子树覆盖中的顶点
3 最小子树覆盖顶点集 Vc = {v |inCover[v] = true， v ∈ V}
4 最小子树覆盖边缘集 Ec = {(u， v） |u， v ∈ Vc， （u， v） ∈ E}
5 返回最小子树 Tc = (Vc、Ec）， 覆盖 W 中的所有顶点
我们现在分析Alg. 4 Minimum-Subtree-Cover的运行时间。对于具有 N 个顶点和 N ? 1 条边的输入树，树上 Cover-DFS 的运行时间复杂度为 O（N + （N ? 1）） = O（N）。追踪最小子树也采用 O（N + （N ? 1））。因此，查找给定顶点子集的最小子树覆盖的总运行时间复杂度为 O（N），与输入树中的顶点数成线性关系。
图5说明了7个发育中的小鼠组织中两个途径的组织特异性。Wnt信号通路具有0.933的高组织特异性，其中不同的组织样品在基因表达中分离良好（图5A）。错配修复途径的组织特异性低，为0.542（图5B）。
thumbnail 下载：
个人电脑幻灯片
巴新放大图片
断续器原始图像
图 5. 在七种小鼠组织发育过程中形成对比组织特异性的两种途径。
连接组织样品的EMST来自每个途径上基因的表达水平。顶点是按组织类型着色的样本。（A）Wnt信号通路在基因表达动力学中具有高组织特异性，同一组织类型的发育样品形成独特的轨迹片段。（B）基因表达动力学中组织特异性低的错配修复途径显示了不同组织类型的混合样品。
https://doi.org/10.1371/journal.pcbi.1009829.g005
算法 5覆盖-DFS（树 T，顶点子集 W，顶点 v)
输入：T，一棵树（V，E），顶点集V和边集E-厦门畜牧期刊杂志论文发表
W，顶点集 V 的子集，用于查找子树覆盖
v，V 中要探索的顶点
1 访问过的马克五
如果顶点 v 是子集 W 的成员，则为 2
3 内覆盖[v] = 真
4 其他
5 ?inCover[v] = false
对于树 T 中与 v 相邻的每个顶点 u，6
7 如果 u 未被访问
8 覆盖-DFS（T， W， u)
9 如果 inCover[u] 为真
10 内覆盖[v] = 真
结果
模拟数据树维试验的性能评价
我们首先评估了Alg 3的性能。Test-Trajectory-Presence（X）识别轨迹模式与更直观的替代方案L和树直径Dm在合成数据上，其中包括来自"Dynverse"轨迹推理项目的229个模拟单单元数据集[13]。这些数据集具有各种类型的轨迹：分岔，多分叉，树和循环，构成地面真实轨迹;我们对每个数据集进行了洗牌，以生成另外229个数据集，这些数据集代表了没有轨迹的地面真实示例。
为了预处理数据，我们使用主成分分析（PCA）执行降维。大多数轨迹推理方法在找到轨迹之前采用降维。PCA可能无法保留轨迹，但它捕获构成数据中主要动态的变化。我们采用CNG碎石测试来选择最佳数量的主成分[14，15]。TDT方法与PCA无关，其他数据预处理程序可以按照我们随附软件的设计进行说明。
应用每种方法生成 458 个分数，每个数据集一个分数。我们使用 1 ? L/N 和 D 的归一化分数m/N 分别表示叶子数和树直径。TDT 效应大小 S/N 用作树维度检验的分数。归一化用于使不同样本大小的数据集具有可比性以进行性能评估。然后，使用评分绘制接收器工作特性（ROC）和精度召回（PR）曲线，并计算每种方法的ROC曲线下面积（AUROC）和PR曲线下面积（AUPR）作为其性能（图6）。在图6A中，TDT的最佳AUROC得分为0.996，而DmAUROC评分为0.951，L的AUROC评分最低为0.826。类似地，图6B显示了PR曲线的AUPR方法的分数。TDT的最高AUPR分数为0.997，而Dm得分为0.954，L的得分为0.756。TDT的高AUC分数证明了树维度测量在识别各向同性模式的轨迹方面比其他更直观的方法更有效。
thumbnail 下载：
个人电脑幻灯片
巴新放大图片
断续器原始图像
图 6. 在 229 个模拟的单单元数据集上进行轨迹存在测试。
（A） TDT效应大小S/N、叶片数量L和树径D的ROOC曲线和AUROC评分m.（B）三种方法的PR曲线和AUPR分数。（C）PCA图为多分岔轨迹模式，点代表细胞和轴基因表达，EMST表示轨迹模式，每个点代表一个细胞。文本显示应用于数据集时每个方法的 p 值。（D） 循环轨迹模式的 PCA 和 EMST 图。（E） 多岔路轨迹的常设仲裁法院和EMST图。（F）–（H） 没有显著轨迹模式的数据集的 PCA 和 EMST 图。
https://doi.org/10.1371/journal.pcbi.1009829.g006
图6C–6E显示了实验中使用的不同类型的轨迹模式，以及由每种方法的p值表示的轨迹存在的统计显著性。每个数据集都由主成分分析 （PCA） 图和相应的 EMST 表示进行汇总。EMST 表示形式保留数据中的结构，即数据点之间的空间关系。我们将 EMST 中的边设置为长度为 1 以捕获拓扑。此 EMST 表示形式足以用于识别轨迹存在的后续步骤。图6C显示了具有轨迹的模拟数据，TDT正确识别了数据中存在轨迹，由1.2e-06的有效p值表示。同样，TDT还识别图6D中p值为1.01e-07的轨迹，以及图6E中p值为9.34e-08的轨迹。叶子的数量L无法识别这些模式，如相对较高的p值所示。树直径 Dm但是，可以识别有效 p 值所指示的轨迹模式。图6F-6H给出了没有轨迹的图案示例。所有方法都正确地认识到这里没有轨迹模式，如相对较高和不太显着的p值所表示的那样。
单单元数据树维试验性能评价
我们在110个已知具有轨迹的真实单细胞数据集上评估了这三种方法，这些数据集也用于"Dynverse"项目[13]。它们构成了具有轨迹的地面真实模式。为了获得对照，我们对基因进行了洗牌，使它们在统计上独立，以生成另外110个在地面真实中没有轨迹模式的数据集。图 7 显示了结果。图7A显示了ROC曲线和每种方法的相应AUROC评分。TDT效应大小具有最佳性能，AUROC得分为0.997，而L的AUROC得分为0.878和Dm在AUROC得分最低为0.846。同样，图7B显示了PR曲线，TDT效应尺寸的最佳AUPR为0.998。L 的最低 AUPR 分数为 0.818 和 DmAUPR 得分为 0.858。这些结果表明了树维检验在从真实单细胞数据中的随机模式识别轨迹方面的有效性和鲁棒性，这些模式具有噪声并受到辍学效应的影响。
thumbnail 下载：
个人电脑幻灯片
巴新放大图片
断续器原始图像
图 7. 在 110 个真实的单细胞数据集上进行轨迹存在测试。
（A） TDT效应大小S/N、叶片数L和树径D的ROC曲线和AUC评分m.（B）三种方法的PR曲线和AUC分数。（C）单细胞人肺腺癌细胞系数据中轨迹模式的PCA图和EMST表示。点表示单元格。文本显示应用于数据集时每个方法的 p 值。（D）人类雌性生殖系单细胞数据中轨迹模式的PCA和EMST图（E） PCA和平面单细胞数据中轨迹模式的EMST图。（F）–（H） 没有显著轨迹模式的数据集的 PCA 和 EMST 图。
https://doi.org/10.1371/journal.pcbi.1009829.g007
图7C–7E显示了由PCA图和EMST表示的真实轨迹模式。图7C是人肺腺癌细胞系数据与轨迹。图7D是人类雌性种系单细胞数据[16]的轨迹。图7E是带有轨迹的平面数据。所有方法都能够识别这些数据集中由显着p值表示的轨迹模式的存在，除了肺癌数据（图7C），其中L的p值不显着。图7E强调了仅使用二维PCA可视化数据中模式的局限性，而EMST表示显示了数据中有趣的动态模式。图7F-7H显示了没有轨迹的图案。没有任何方法在统计上验证这些数据集中是否存在由相对较高的 p 值表示的轨迹。
模拟大型数据集上树维检验的经验运行时
一个大型单细胞测序数据集可以包含一百万个细胞（N = 106).一个精确的EMST解，采用二次时间O（N2） 是不切实际的。因此，我们的研究使用已建立的双树Boruvka算法[11]来查找时间复杂度O（N log N）中的近似EMST。图8显示了TDT在具有不同数量单元的模拟数据上使用近似Boruvka算法的经验运行时间。此运行时包括 100 个模拟，用于建立可以预先计算的空分布参数。结果表明，实际运行时间大致线性地扩展到N。数十万个细胞的效率是合理的，正如人们所期望的单细胞分析一样。
thumbnail 下载：
个人电脑幻灯片
巴新放大图片
断续器原始图像
图 8. TDT对具有不同细胞数的模拟单细胞数据的经验运行时间。
运行时包括 100 个模拟，用于估计空分布的参数。水平轴表示像元数。垂直轴是以分钟为单位的运行时间。时间是在2015年Apple Macbook Pro笔记本电脑上使用单个线程的2.2 GHz四核Intel Core i7处理器上记录的。
https://doi.org/10.1371/journal.pcbi.1009829.g008
此外，在没有模拟的情况下，在输入数据上计算轨迹存在的TDT效应大小，因此在计算上比p值快得多（在显示时间的约1%内）。在实际应用中，可能只需要计算具有相当强效应的数据的 p 值。
通路在发育过程中表现出与细胞轨迹不同的轨迹模式
测试轨迹模式是否存在的能力开辟了检查细胞动力学之外的通路动力学的可能性。通路轨迹由该通路中的基因跨越。我们研究了小鼠胚胎干细胞单细胞数据集中的通路轨迹是否可能与细胞轨迹动态不同[17]。图9A是聚类为五组的细胞的PCA图。图9B是细胞的EMST表示，显示强轨迹模式，按树维度显示显着的p值。图9C–9F是四种途径基因表达动态的EMST表示。这些图突出显示了细胞和通路轨迹之间EMST表示的结构差异。颜色编码表明，细胞的空间关系既不保留在细胞和通路EMST之间，也不保留在通路EMST之间。这表明通路轨迹模式确实可以与细胞轨迹模式区分开来。
thumbnail 下载：
个人电脑幻灯片
巴新放大图片
断续器原始图像
图 9. 胚胎干细胞单细胞数据中不同的细胞和通路轨迹。
单元格被聚类为由颜色指示的五组。（A）前三个主要组分中观察到胚胎干细胞的基因表达数据。（B）整个转录组的EMST表明细胞水平的强轨迹模式。（C）-（F）胚胎干细胞中轨迹模式的EMST表示，使用（C）肾素，（D）蛋白酶体，（E）不匹配修复和（F）刺猬途径上的基因表达。
https://doi.org/10.1371/journal.pcbi.1009829.g009
在哺乳动物发育过程中，基因表达动态的组织特异性在通路中差异显著-厦门畜牧期刊杂志论文发表
现在我们通过组织特异性来表征通路动力学。如果通路中基因的表达在来自相同组织类型的样品中是同一的，但在来自其他组织类型的样品中是异质的，则该途径是组织特异性的。因此，如果来自该组织的样本在轨迹模式的EMST表示中处于同一邻域，则该通路是组织特异性的。我们利用这一特性来量化数据集中所有组织类型的途径的平均组织特异性。我们检查了哺乳动物发育转录组的集合，涵盖了人类和小鼠从胚胎阶段到成年期的七种组织类型[18]，包括大脑，小脑，心脏，肾脏，肝脏，卵巢和睾丸。
我们测试了轨迹模式的存在，并测量了40个KEGG途径的组织特异性。图10A根据小鼠和人类发育过程中的平均组织特异性对途径进行了排名。钙信号传导是组织特异性最高的途径，平均组织特异性为 0.899。核糖体途径的组织特异性最低，平均组织特异性为0.533。图10B可视化了人与小鼠之间通路组织特异性是如何保守的。40条途径的组织特异性在人和小鼠之间高度相关，斯皮尔曼秩相关系数为0.763（p值= 1.69×10?9).
thumbnail 下载：
个人电脑幻灯片
巴新放大图片
断续器原始图像
图 10. 人类和小鼠发育过程中40种生物学途径基因表达动态的组织特异性。
（A）途径是根据人类和小鼠七种组织类型的转录组数据计算的哺乳动物发育过程中的平均组织特异性排名的。（B）发育中的小鼠和人之间的通路组织特异性是保守的。每个点都是一个途径，轴分别代表其小鼠和人类的组织特异性评分。
https://doi.org/10.1371/journal.pcbi.1009829.g010
具有组织特异性表达动力学最多和最少的途径
图11分别显示了人类和小鼠发育过程中两种组织特异性最强和两种最少的途径的基因表达动态。Wnt和钙信号通路都呈现出强烈的轨迹模式，这些模式也具有高度的组织特异性，如覆盖独特组织类型的轨迹片段所示;错配修复和核糖体途径呈现弱的轨迹模式，同时表现出相对较低的组织特异性。
thumbnail 下载：
个人电脑幻灯片
巴新放大图片
断续器原始图像
图 11. 哺乳动物发育过程中大多数和最少组织特异性途径的基因表达动态。
七种组织类型按颜色编码。人类发育过程中（A）钙信号传导，（B）Wnt信号传导，（C）错配修复和（D）核糖体途径上基因表达的PCA和EMST。（E）–（H） PCA和EMST在小鼠发育过程中相同四个途径的基因表达动态。
https://doi.org/10.1371/journal.pcbi.1009829.g011
通路组织特异性是来自同一组织类型的样品之间通路基因表达动态的均质性程度。在图11中，钙信号传导和Wnt信号通路在哺乳动物发育过程中具有高度的组织特异性：来自相同组织类型的样品在轨迹模式的EMST表示中紧密地结合在一起。大多数高度组织特异性的途径都参与信号转导。
钙信号传导在整个发育过程中很重要，可调节细胞过程，如分裂、迁移、死亡和分化[19]。其机制在调节不同组织类型的钙振荡频率方面具有独特性。研究已经确定了对组织类型具有特异性的多种钙受体亚型[20]。由于这些受体控制钙离子向细胞的运输，受体的多种亚型导致不同组织中独特的下游钙信号。例如，已观察到钙信号传导显示线粒体钙摄取中的组织特异性[21]。
众所周知，Wnt信令是通用的，对发展至关重要。因此，长期以来一直认为该途径导致来自不同组织类型的均匀反应。然而，新出现的证据表明，该途径驱动不同组织类型的独特反应，并被推测为组织特异性[22]。例如，发现Wnt信号通路对收缩压具有独特的易感位点[23]。
发育过程中组织特异性低，修复和核糖体途径不匹配（图11）。许多这样的途径都参与遗传信息处理。
我们的发现揭示了几种信号转导途径的组织特异性高，但在遗传信息处理途径中的组织特异性较低，这部分得到了有限文献的支持。需要进一步的生物学调查来支持或反驳我们关于途径组织特异性的发现，这一方向仍未得到充分探索。
讨论
树维数作为轨迹模式强度的有利统计量
多变量数据的EMST表示既具有计算效率又具有全局最优性，被弹弓[1]等轨迹推理方法广泛使用。我们的前提是，数据中动态模式的存在可以通过其相应EMST的线性程度来表征。我们评估的所有三个EMST统计数据都试图表征EMST的线性程度，实际上，轨迹模式的存在。所有这些都能够区分线性的极端情况。树维度量 Td如果EMST是线性的，则叶子L的数量最小化，如果它是星树，则叶子的数量最大化;树直径 Dm由线性EMST最大化，由星形树最小化。然而，它们的不同之处在于它们所捕获的全球和局部细节。一方面，由局部细节驱动的叶子数量L不会区分线性轨迹和分岔轨迹，其叶子数量与噪声状态相同。另一方面，树直径Dm，由全局最长路径驱动，如果两者具有相同的最长路径长度，则无法区分具有少量长分支的树和具有许多短分支的树。树维度度量值旨在捕获树的方向性/线性度。它提供了更平衡的线性度量，集成了更直观的L和Dm，分别容易受到局部和全局极端的影响。Td尝试结合D的优势m和L，同时减轻其局限性。正如我们在仿真研究中所证明的那样，这使得树维度在性能上比其他替代方案更具优势。
欧几里得最小生成树的使用
我们对EMST的选择取决于许多考虑因素。首先，EMST 是一个全局最优紧致结构，用于汇总欧几里得空间中的多变量数据。其次，它没有可调参数来显著更改相同数据上的表示形式。第三，它在轨迹推理应用中得到了广泛的测试：一些最流行的轨迹推理方法，如TSCAN [5]和Monocle [24]，采用MST进行聚类分析来捕获数据中的轨迹。
值得注意的是，对于EMST被认为不相关的特定应用程序，我们的方法仍然适用于数据集的任何其他类型的树摘要，或由数据集的图形表示产生的树摘要。只要根据某些定义的邻近度量推导出树中的顶点，轨迹测试的上下文就仍然有效。
在细胞簇上找到MST是一个有希望的方向。尽管效率很高，但基于我们早期未发表的工作的轨迹存在测试的性能并不像我们在这里介绍的那样好。性能取决于聚类形状和聚类数。流行的聚类方法在我们的背景下似乎是无效的，这意味着需要制定轨迹友好的聚类策略。
球面多元正态分布与完美线性轨迹之间
我们使用球面 MVN 分布来表示用于为树维度检验统计量生成空分布的数据。球面MVN表示没有任何方向或方向的各向同性图案，因此没有轨迹。另一个极端是完美的线性模式，其中EMST是一个路径图，对应于一维流形。虽然 TDT 旨在优先处理具有内在维数为 1 的轨迹，但它能够捕获从高度线性轨迹到各向同性模式（如球面 MVN 分布所示）的频谱上的情况。因此，如果数据驻留在具有较低内在维数的子空间中，而不是一个高维空间，我们的方法可能会接受轨迹存在。尽管维数可能大于 1，但模式可能比各向同性更接近轨迹。
树维度检验统计量的空分布
树维统计量 S 的对数正态零分布的参数（均值和标准差）取决于样本数量和数据集的维数。我们的方法从零总体中抽取样本，以获得参数的最大似然估计值。如果输入数据集具有较大的样本量和高维数，则采样是计算密集型的，就像单像元数据通常的情况一样。将来可能的工作是通过推导检验统计量 S 的理论渐近零分布来避免采样。
通路轨迹可能不遵循细胞轨迹
由于通路空间由基因子集跨越，因此很明显，来自通路的点云的动力学可能与所有基因跨越的点云的动力学不同。因此，如果使用整个基因集，通路可能会揭示不存在的强动态模式。也有可能并非所有途径都会导致强烈的轨迹模式。
在图9中，通路动力学不同于所有基因的全局动力学。此外，按组织类型进行颜色编码的点突出了不同途径之间空间分布的差异。集中在一个途径的表达动力学亚空间中的一个邻域中的点可能分散在另一个途径的表达动力学的子空间中，从而产生不同的模式。
结论
我们开发了树维数测试来推断多变量数据中的轨迹存在。我们对模拟和生物学数据的研究验证了该方法相对于其他选项的有效性。该方法能够通过超出整体细胞动力学的轨迹强度来确定途径优先级。通过所提出的亚树覆盖方法，我们进一步发现了哺乳动物发育过程中途径动力学的组织特异性的显着差异：几种信号转导途径在基因表达动力学中具有高度的组织特异性，而遗传信息处理途径在组织特异性方面往往较低。我们的工作补充了现有的轨迹推断方法，为潜在轨迹模式的重要性提供统计支持，也打开了一扇窗，通过动力学确定路径的优先级，以揭示生物系统的详细分子机制。
引用
1.Street K， Risso D， Fletcher RB， Das D， Ngai J， Yosef N， et al.弹弓：单细胞转录组学的细胞谱系和伪时间推断。BMC 基因组学。2018;19(1):477.pmid：29914354
查看文章PubMed/NCBI谷歌学术搜索
2.Tran TN，Bader GD. Tempora：使用时间序列单细胞RNA测序数据的细胞轨迹推断。PLoS计算机生物学 2020;16（9）：e1008205.pmid：32903255
查看文章PubMed/NCBI谷歌学术搜索
3.Todorov H，Cannoodt R，Saelens W，Saeys Y. TinGa：快速灵活的轨迹推断与不断增长的神经气体。生物信息学。2020;36（Supplement_1）：i66–i74.下午：32657409
查看文章PubMed/NCBI谷歌学术搜索
4.Saelens W， Cannoodt R， Todorov H， Saeys Y.单细胞轨迹推理方法的比较。纳特生物技术。2019;37:547–554.下午：30936559
查看文章PubMed/NCBI谷歌学术搜索
5.吉志，吉H.单细胞RNA-seq分析中的假时间重建与评估。核酸研究。2016;44（13）：e117–e117.pmid：27179027
查看文章PubMed/NCBI谷歌学术搜索
6.Welch JD， Hartemink AJ， Prins JF.SLICER：从单细胞RNA-seq数据推断分支的非线性细胞轨迹。基因组生物学。2016;7(1):106.pmid：27215581
查看文章PubMed/NCBI谷歌学术搜索
7.查扎尔 F. 米歇尔·拓扑数据分析简介：数据科学家的基本和实践方面。arXiv.2017;第1710.04019页。
8.Wolf FA， Hamey FK， Plass M， Solana J， Dahlin JS， G?ttgens B， et al. PAGA： 图抽象通过单细胞的拓扑保存图来调和聚类与轨迹推理。基因组生物学。2019;20(1):59.pmid：30890159
查看文章PubMed/NCBI谷歌学术搜索
9.Vandaele R， Saeys Y， Bie TD.通过森林表示在图形中挖掘拓扑结构。机器学习研究杂志。2020;21(215):1–68.
查看文章谷歌学术搜索
10.Vandaele R， Rieck B， Saeys Y， De Bie T. 通过图形近似对度量树进行稳定的拓扑特征。模式识别字母。2021;147:85–92.
查看文章谷歌学术搜索
11.March WB， Ram P， Gray AG.快速欧几里得最小生成树：算法、分析和应用程序。在：第16届ACM SIGKDD知识发现和数据挖掘国际会议论文集。KDD'10.美国纽约州纽约：计算机协会;2010. 第603–612页.可从： https://doi.org/10.1145/1835804.1835882.
12.Curtin RR， Edel M， Lozhnikov M， Mentekidis Y， Ghaisas S， Zhang S. mlpack 3： 一个快速、灵活的机器学习库。开源软件学报.2018;3:726.
查看文章谷歌学术搜索
13.Cannoodt R， Saelens W， Todorov H， Saeys Y. 包含轨迹的单细胞组学数据集 （https://doi.org/10.5281/zenodo.1443566）;2018. 可用于： https://doi.org/10.5281/zenodo.1443566.
14.卡特尔 RB.因子数的碎石测试。多元行为研究。1966;1(2):245–276.下午：26828106
查看文章PubMed/NCBI谷歌学术搜索
15.Gorsuch R，Nelson J. CNG碎石测试：确定因素数量的客观程序。在：在多变量实验心理学学会年会上发表;1981.
16.李磊， 董军， 闫丽， 永军， 刘旭， 胡毅， 等.单细胞RNA-seq分析绘制了人类生殖系细胞和性腺生态位相互作用的发展图。细胞干细胞。2017;20(6):858–873.pmid：28457750
查看文章PubMed/NCBI谷歌学术搜索
17.Ochiai H， Hayashi T， Umeda M， Yoshimura M， Harada A， Shimizu Y， et al.小鼠胚胎干细胞中转录爆发的全基因组动力学特性。科学进展.2020;6(25).下午：32596448
查看文章PubMed/NCBI谷歌学术搜索
18.Cardoso-Moreira M， Halbert J， Valloton D， Velten B， Chen C， Shao Y， et al.哺乳动物器官发育中的基因表达。自然界。2019;571(7766):505–509.pmid：31243369
查看文章PubMed/NCBI谷歌学术搜索
19.Brodskiy PA，Zartman JJ.钙作为上皮组织发育中的信号整合器。植物学 2018;15（5）：051001.pmid：29611534
查看文章PubMed/NCBI谷歌学术搜索
20.Bootman MD， Berridge MJ.钙信号传导的基本原理。细胞。1995;83(5):675–678.下午：8521483
查看文章PubMed/NCBI谷歌学术搜索
21.Paillard M， Csordás G， Szanda G， Golenár T， Debattisti V， Bartok A， et al.细胞质Ca（2+）信号的组织特异性线粒体解码由MICU1/2和MCU的化学计量控制。细胞代表 2017;18（10）：2291–2300.pmid：28273446
查看文章PubMed/NCBI谷歌学术搜索
22.索德霍尔姆 S， 坎图 C.WNT/β-连环蛋白依赖性转录：组织特异性业务。WIREs系统生物学和医学。2020;第1511页。pmid：33085215
查看文章PubMed/NCBI谷歌学术搜索-厦门畜牧期刊杂志论文发表
23.赵毅， 布伦科韦 M， 石霞， 舒磊， 李维安 C， 安 IS， 等.整合基因组学分析揭示了组织特异性途径，网络和血压调节的关键调节因子。心血管医学前沿。2019;6:21.pmid：30931314
查看文章PubMed/NCBI谷歌学术搜索
24.Trapnell C， Cacchiarelli D， Grimsby J， Pokharel P， Li S， Morse M， et al.细胞命运决定的动力学和调节因子通过单个细胞的伪时间顺序揭示。纳特生物技术。2014;32(4):381–386.下午：24658644
查看文章PubMed/NCBI谷歌学术搜索