TimeFlow 2：一种用于流式细胞术数据的无监督细胞谱系检测方法

《Cytometry Part A》：TimeFlow 2: An Unsupervised Cell Lineage Detection Method for Flow Cytometry Data

本文推荐了一种名为TimeFlow 2的新型无监督计算方法，用于从单次静态流式细胞术（FC）或质谱流式（MC）数据快照中推断细胞分化轨迹（谱系）。该方法无需预先设定通路数量、细胞类型或时间标签，通过伪时间（pseudotime）分段、高斯混合模型（GMM）聚类、基于最优传输（OT）的路径分组等步骤，自动构建细胞状态级分化路径，并在健康骨髓样本中成功解析了单核细胞、中性粒细胞、红细胞和B细胞等不同成熟阶段的四个独立谱系，其性能在FC数据集上优于其他方法，在MC数据集上保持竞争力。

1 引言

细胞轨迹推断（TI）方法旨在从细胞分化动态过程中推断细胞谱系通路。这些方法基于细胞基因表达谱的相似性，为每个细胞计算一个称为伪时间的数值，并利用该值将细胞从最不分化到最分化进行排序。轨迹推断方法的一个典型应用是造血过程，即血细胞形成的过程。该过程发生在骨髓（BM）中，使造血干细胞和祖细胞（HSPCs）能够产生功能不同的血细胞类型。TI方法在研究HSPCs命运决定、造血基因调控网络、健康和失调模式（如急性髓系白血病（AML））的比较以及不同年龄组基因表达改变的分析中做出了重要贡献。

虽然大多数TI应用使用RNA测序（RNA-seq）数据，但它们在流式细胞术（FC）或质谱流式（MC）数据中的应用也引起了兴趣。轨迹结构取决于该过程中涉及的完全分化细胞类型的数量。细胞可能向单一成熟细胞类型进展（线性轨迹），沿着终止于不同细胞类型的通路发散（树状和多分支轨迹），或在断开的图上演化。不同的TI方法专门处理每种情况。方法也可以根据输入数据的形式进行划分：要么是一个静态细胞快照（假定捕获了过程中的所有过渡状态），要么是在不同实验时间点采集的多个快照（时间序列数据集）。

许多TI方法使用单个静态快照。它们通过在细胞簇上构建分段线性最小生成树（MST）或使用主曲线拟合和主流形算法构建更平滑的MST来检测分支。或者，它们用k近邻图（k-NNG）上的节点表示细胞，并模拟随机游走来构建吸收马尔可夫链，将吸收态视为成熟细胞状态。其他基于图的TI方法使用启发式方法，如随机游走的相关性、图节点中心性度量、最小成本流算法、生长神经气体、社区检测和统计连通性度量、分层混合建模以及用于随机游走聚类的持久同源图。

需要具有时间标签的多个快照的TI方法对细胞群体动力学进行建模，假设细胞以最小成本方式分化，因此解决离散最优传输（OT）问题以推断连续快照中细胞群体之间的概率耦合。从时间序列数据中学习连续细胞动力学的方法使用其他OT框架，例如动态OT公式、Jordan–Kinderlehrer–Otto流方案和薛定谔桥。然而，这些方法没有明确地将细胞与谱系通路关联起来。StationaryOT和MultistageOT提出了用于单个静态快照的OT耦合，但两者都需要手动选择的成熟细胞来推断轨迹。

本研究专注于使用通过FC获得的单个静态快照进行多分支轨迹中的细胞谱系检测。专门为FC/MC数据设计的TI方法包括用于线性轨迹的Wanderlust、用于单分支轨迹的Wishbone、一种使Wanderlust适应树状轨迹但需要手动指定完全分化细胞的方法、CytoTree（在大型数据集中计算伪时间时需要下采样）、t-Space和tviblindi。其他相关方法是TrackSOM和ChronoClust，它们跟踪时间序列数据集中的细胞群体。为了克服诸如先验知识、下采样和对时间序列数据的依赖性等挑战，我们扩展了我们之前的伪时间计算方法TimeFlow，提出了一种用于大型细胞术数据的无监督细胞谱系检测方法。我们的方法将每个细胞与一个通路相关联，允许某些细胞在不同的通路之间共享，反映了它们的共同祖先。我们展示了TimeFlow 2在使用七个先前发表的BM数据集解析造血细胞谱系方面的潜力，并讨论了其在建模标记物动力学方面的有用性。我们还将我们的方法与其他知名的自动检测细胞谱系的TI方法进行了比较，并根据真实门控标签解释了结果的生物学相关性。

2 材料与方法

2.1 问题定义

我们首先形式化我们的谱系检测方法的目标和设置。给定一个标记物表达矩阵，其中每一行对应一个细胞，每一列对应一个FC标记物，如分化簇（CD）、表面或细胞内标记物，我们的目标是将每个细胞与一个在[0,1]范围内的伪时间值和一个谱系通路相关联。预期输出是一组元组，每个元组对应一个推断的通路。每个元组包括构成特定通路的细胞的索引，以及这些细胞的伪时间。形式上，我们将一个元组表示为一个序列，其中包含细胞的索引、其伪时间以及通路中的细胞总数。这是一种直接的输出来拟合每个通路中标记物表达作为伪时间函数。如前所述，一些细胞可能被分配到一个以上的通路中（例如，尚未承诺向特定谱系发展的多能细胞）。通路的数量不是预先固定的，而是由方法自动确定。

2.2 方法概述

TimeFlow解决了细胞排序任务。它通过计算在k-NNG上从细胞到轨迹根（即指定的具有高CD34表达的干细胞）的密度驱动最短路径，为每个细胞分配一个伪时间值。这样，它保留了局部邻域内细胞转换的连续性。然而，伪时间排序不足以重建谱系。在TimeFlow 2中，我们采用一种方法来近似全局轨迹结构，该方法保留了（a）来自伪时间的排序信息，和（b）来自细胞特征（此处为标记物）的生物学信息。我们通过将伪时间轴划分为几个段并在每个段内定义粗略的细胞状态来解释（a）。这对于理解这些有序的细胞状态如何沿着伪时间分化为不同的分支至关重要。为了满足（b），我们假设向相同谱系承诺的细胞通过相似的路径分化。因此，我们在细胞状态级别构建路径，然后根据它们的相似性将其聚合成生物学相关的组。尽管只能访问单个细胞快照，标准（a）和（b）允许我们同时在伪时间和标记物空间中跟踪细胞状态路径的演化。更具体地说，TimeFlow 2采取以下步骤：

1. 1.
   伪时间分割：将伪时间轴分割成等宽、部分重叠的段。
2. 2.
   每段细胞聚类：在每个段内对细胞进行聚类。与几种其他TI方法强制整个快照中的细胞进行聚类不同，这里我们将聚类限制在每个段内以获得有序的簇或细胞状态。
3. 3.
   路径构建：基于最大细胞重叠连接连续段之间的细胞簇，以构建细胞状态级别的路径。
4. 4.
   路径分组：基于OT的相似性度量对路径进行分组，以构建轨迹主干。
5. 5.
   轨迹细化：细化轨迹主干。

图1提供了TimeFlow 2的示意图概述。

2.3 方法细节

2.3.1 伪时间分割

我们将伪时间轴划分为等宽的段，使得每对连续段按其宽度的百分比重叠。我们允许段之间的细胞重叠，以保留转换的连续性，避免严格的边界，并随后指导细胞状态连接。因此，一个段的下界位于其前一个段的上界之前。我们用索引表示每个伪时间段，其中总段数。设是所有伪时间段的细胞集合。伪时间段的默认值为，段宽度重叠的默认值为35%。

2.3.2 每段细胞聚类

我们使用高斯混合模型（GMM）对每个段内的细胞群体进行建模，假设每个高斯分量描述一个粗略的细胞状态。在细胞术中，GMM及其变体在自动门控/聚类免疫细胞类型方面特别有用。在为每个段估计GMM之后，我们获得一个混合模型，其中表示段索引，表示段中的总簇数，是具有均值向量和协方差矩阵的高斯分布，是每个高斯分量的权重，满足。我们注意到每个分量都有自己的完整协方差矩阵，并由其相关的概率度量来表征。我们设置每个段的最大簇数，并使用贝叶斯信息准则（BIC）来定义每个段的高斯分量数量。

2.3.3 路径构建

在此步骤中，我们连接簇以定义细胞状态转换，并在细胞状态级别构建路径。从最终段开始，我们根据它们共享的最大细胞数将段中的每个簇连接到前一个段中的一个簇。我们向后遍历伪时间段，以确保每个簇都连接到一个更早的簇，形成回溯路径。考虑到使用GMM进行聚类，每个路径由跨伪时间段的高斯分量组成，并用路径索引、路径总数、段中的完整GMM模型、路径遍历的特定高斯分量以及路径遍历的总段数来表示。为了计算路径总数，我们计算所有段中的总簇数，并减去第一段中的簇数。

2.3.4 路径分组

为了定义轨迹主干，我们将终止于最终段的路径聚合成生物学信息丰富的组。我们提出一个问题：将一条路径转换为另一条路径需要多少努力。我们通过考虑一个成本函数来量化这一点，该函数对路径在对应段的状态之间的距离求和。我们将总成本解释为这些路径之间的相似度。小成本表示路径的细胞组成相似，因此谱系一致，而大成本意味着它们的细胞组成不同，可能谱系分化。为了测量路径相似性，我们依赖最优传输（OT）。OT解决了两个概率度量（或分布）之间的质量传输问题，通过基于传输成本函数的传输计划将它们耦合起来。当传输成本函数由范数给定时，OT定义了概率度量之间的距离。这被称为p-瓦瑟斯坦距离，它描述了将一个概率度量的单位质量传输到另一个概率度量的单位质量的最小成本。我们在支持信息第S1节中提供了关于OT的更多细节，并请读者参阅文献了解OT的理论基础，以及文献了解基于计算的方法。我们将两条路径之间的总路径距离表示为它们对应高斯分量之间的2-瓦瑟斯坦距离（）之和。对于两条等长路径和，距离由公式给出，其中和是两条长度为的路径。因此，两条等长路径之间的成本遍历它们的所有段，计算它们对应高斯分量之间的。重要的是，为了计算这些段内距离，我们依赖一个封闭形式的解决方案，使我们免于以复杂度计算完整的OT计划。在两种边缘分布都遵循（单/多元）高斯分布的特殊情况下（如此处的情况），可以根据分布的均值向量和协方差矩阵明确写出，如文献中所述。对于和，2-瓦瑟斯坦距离是公式。

成对路径成本作为无监督路径分组的输入，而无需访问实际的组数。从没有组开始，我们遍历路径成本，如果它们的成本低于某个阈值，则合并两条路径。如果一对路径满足此条件，但其路径均不属于任何现有组，则创建一个新组。否则，如果任何路径已分配给一个组，则剩余路径将合并到该组中。因此，与其他路径的成本超过阈值的路径将不会分配给任何现有组，而是形成自己的路径组。在上述过程结束时，属于同一组的路径接收一个标签（例如，整数）。这些标签是人为的，不包含生物学信息。它们用于指示给定阈值的组成员资格。形成自己单独组的路径也会收到一个标签。考虑到为基于的成本选择截止阈值并不明显，我们对一系列候选阈值重复相同的过程。该路径分组过程的评估基于模型本身，因为没有可用的真实标签。我们选择最大化轮廓系数（SI）的阈值作为最佳阈值。为了计算SI，我们使用scikit-learn的实现，并输入人工标签和成本。较高的距离阈值比较低的严格阈值对路径合并更宽容。据我们理解，上述通路分组策略尚未在TI的背景下使用。SI不用于聚类验证，而是用于评估基因表达数据任意组的一个例子见文献。

2.3.5 轨迹细化

在下一步中，我们将所有不终止于最终段的细胞状态路径整合到轨迹主干中。这些细微的路径可能是虚假的，或者描述了向新的终末状态的成熟。因此，它们要么被先前发现的通路吸收，要么在主干中被视为新的候选通路。为了识别此类路径，我们搜索在任何段（最终段除外）中没有前向连接的簇。为简单起见，我们将这些终端簇称为“叶簇”。此外，我们列出属于任何主要通路的所有簇，并将其称为“已探索簇”。我们区分以下以叶簇结尾的路径情况：

1. 1.
   对应于未发现通路的路径，其成熟细胞集中在早于最终段的段中。
2. 2.
   对应于先前发现通路的路径。这涉及汇聚到与已探索簇高度相似的叶簇的路径，原因可能是不同的分化速度或偏离主要通路的随机偏差。考虑到BIC倾向于支持具有大量簇的模型（尤其是在后面的段中），这些偏差是极有可能的。此外，细胞快照包括处于其自身个体成熟不同状态的细胞，因此轨迹中预计会出现变异性和噪声。

为了区分标记潜在新终末细胞状态的叶簇和可吸收簇，我们使用段内距离，并检查每个叶簇的以下哪种情况为真：

1. 1.
   前N个（默认为4）最近簇中没有一个属于已探索簇列表。
2. 2.
   前1个较近簇是已探索簇，但与叶簇不共享高度相关的特征。
3. 3.
   前1个较近簇是已探索簇，与叶簇共享高度相关的特征。

一方面，情况1和2促使将叶簇的路径作为新的主要通路纳入。另一方面，情况3促使叶簇的路径被一个主要通路吸收（在考虑用新通路更新轨迹之后）。在情况1中，叶簇可能代表先前未见的成熟状态。我们要求所有前4个最近簇都是未探索的，以降低轨迹过度分支的风险，尤其是在早期段中。在情况2中，我们将与其已探索的前1个最近簇没有充分相关性的叶簇视为潜在的新终末状态，以防止轨迹分支不足。相反，属于情况3的叶簇，以及不满足情况1条件的叶簇，将被现有通路吸收。我们将叶簇的路径视为时间序列，其中时间步数等于路径遍历的段数。每个时间步由一个维特征向量表征，特征值基于落入其簇的细胞的平均标记物表达计算。我们使用每个叶簇路径的时间序列表示来计算其与前1个最近簇时间序列的皮尔逊相关性，并应用0.85的严格阈值来区分情况2和3。如果所有段间相关性的第10个百分位数高于0.85，表明即使最弱的相关性也保持高度相似性，则可以放宽阈值。叶簇基于快速动态时间规整（DTW）被现有通路吸收，DTW找到两个多元时间序列之间的最优匹配，而不受其长度影响。

在我们方法的最后一步，我们检查每个通路中终末细胞之间的层次关系，以降低虚假通路的风险。我们检索每个通路中具有最大伪时间值的细胞（默认为5%）作为其最成熟细胞的代表，并根据它们所属的通路（通路标签）对它们进行标记。我们计算这些终末细胞的标记物汇总统计量（均值、中位数、第25和第75分位数），遵循细胞术中特征工程的其他示例，并进行层次聚类。随之而来的是，选择切割树层次结构（树状图）的阈值决定了聚类的粒度。过低的阈值可能导致状态过度碎片化，而过高的阈值可能使轨迹过于简化。使用通路标签（不同于门控标签），我们计算不同阈值下的轮廓系数，并选择具有最大值的阈值（如第2.3.4节所述）。我们注意到，我们的方法允许用户根据他们所需的分辨率来改变此阈值。可选地，用户可以查阅每个通路标记物与伪时间的相关散点图，并尝试不同的阈值，如第3.5节所述。

2.3.6 使用FlowSOM的TimeFlow 2变体

TimeFlow 2是模块化的，允许替换其组件以获得更大的灵活性。作为GMM聚类的替代方案，我们考虑FlowSOM，这是一种成熟的非参数方法，用于根据基准测试对细胞术数据进行门控和可视化。与自组织映射（SOMs）类似，FlowSOM训练一个SOM节点网格，并通过将每个数据点（细胞）分配给其最近节点来形成簇。然后，它更新该SOM节点及其邻居的位置以保留数据的拓扑结构。为了便于与用Python 3.11.3编写的TimeFlow 2集成，我们选择了FlowSOM的Python版本。我们将此变体称为TimeFlow 2 (FlowSOM)，以区别于TimeFlow 2 (GMM)。鉴于FlowSOM不假设簇服从高斯分布，我们无法应用公式（3）。我们通过求解熵正则化OT问题，使用Sinkhorn算法快速、可扩展地近似OT计划（支持信息第S1节），来计算SOM节点的段内距离。我们将TimeFlow 2 (FlowSOM)的默认超参数设置为10个段，35%的段重叠和5×5 SOM网格（25个簇），以匹配TimeFlow 2 (GMM)的默认设置，并对OT计划使用1e-1的熵正则化项和2048的批量大小。我们在第3.3节讨论了两个变体对超参数设置的敏感性。

2.4 数据集

我们将TimeFlow 2应用于七个来自无血液病患者的BM样本，如原始出版物中所述。我们重复使用了文献中的P1/2/3-BM FC数据集，这些数据集范围从499,544到597,613个细胞，以20个CD标记物为特征（支持信息第S2节）。在不同成熟阶段的细胞中，我们通过门控识别出15%–20%的成熟中性粒细胞，3%–6%的成熟单核细胞，0.3%–6%的成熟红细胞和0.7%–3%的成熟B细胞。所有细胞计数在表S1中给出。此外，我们考虑了其他四个大型、公开可用的MC BM数据集。我们使用了在文献中发现的Kimmey 6814（347,167个细胞）和Kimmey 6796（647,730个细胞）BM数据集，这些数据集由作者基于SPADE簇进行了注释。两个数据集包括33个细胞群体，其中包括以下细胞类型：中性粒细胞（23%）、CD4记忆T细胞（8%–9%）、CD8记忆T细胞（8%–9%）、自然杀伤细胞（NKs，7%–10%）、B细胞（4%–6%）、自然杀伤T细胞（NKT，1.1%–3%）、单核细胞（1.1%–2%）、血小板（0.9%–1.1%）、常规树突状细胞（cDCs，0.9%–1.6%）、浆细胞样树突状细胞（pDCs，0.7%）、嗜碱性粒细胞（Ba，0.5%）、多染性幼红细胞（0.2%–0.6%）和浆细胞（0.1%–0.4%）。表S2和S3分别总结了Kimmey 6814/6796数据集的细胞计数。支持信息第S2节列出了我们为这些数据集的分析保留的32个抗体靶点。我们还使用了Levine 13和Levine 32数据集，它们分别用13和32个CD标记物测量了167,044和265,627个细胞的细胞表达。Levine 13包含25个细胞群体，包括成熟B细胞（4%）、单核细胞（4%）、成熟CD4 T细胞（8%）、幼红细胞（7%）、成熟CD8 T细胞（4%）、NK细胞（2%）、巨核细胞（2%）、浆细胞（0.2%）、pDCs（0.1%）、血小板（0.003%，5个细胞），而Levine 32包含15个群体，其中CD4 T细胞（9%）、CD8 T细胞（7%）、单核细胞（7%）、成熟B细胞（6%）、CD16阴性NK细胞（1%）、CD16阳性NK细胞（0.8%）、嗜碱性粒细胞（0.4%）、pDCs（0.4%）和浆细胞（0.1%）。我们注意到，Levine 13中50%的细胞和Levine 32中60%的细胞在原始出版物中未被表征并标记为“未分配”。支持信息第S2节提供了有关这些数据集所用panel的更多细节，表S4和S5提供了细胞分数。遵循标准预处理实践，我们使用反双曲正弦（ArcSinh）和系数5转换所有MC数据集，并进一步将其缩放到[0,1]。

2.5 评估指标

为了评估TimeFlow 2和其他TI方法推断的谱系的生物学相关性，我们使用了作者通过手动或自动门控提供的细胞类型标签。具体来说，我们根据原始出版物和先前对造血文献的引用定义了参考谱系。例如，在P1/2/3-BM数据集中，我们预期有以下四个谱系：

1. 1.
   单核细胞：Immature Mono → Intermediate Mono → Mature Mono。
2. 2.
   中性粒细胞：Immature Neu → Intermediate I Neu → Intermediate II Neu → Mature Neu。
3. 3.
   红细胞：Immature Ery → Intermediate I Ery → Intermediate II Ery → Mature Ery。
4. 4.
   B细胞：Immature B-cells → Intermediate