Wishbone identifies bifurcating developmental trajectories from single-cell data
1、发表时间：2016年4月
2、期刊：Nature biotechnology
3、代码链接：https://github.com/dpeerlab/wishbone/blob/master/notebooks/Wishbone_for_single_cell_RNAseq.ipynb
4、用途：推断细胞分化关系、辅助鉴别细胞类型
5、正文：
摘要：最近的单细胞分析技术为阐明细胞分化路径提供了前所未有的机会。在这里，我们介绍了Wishbone，这是一种用于沿高分辨率的分叉发展轨迹定位单细胞的算法。Wishbone使用多维单细胞数据，例如mass cytometry（大规模流式细胞仪）或RNA-Seq数据作为输入，并根据细胞的发育进程对其进行排序，并通过将每个细胞标记为分叉前或分叉后的两个细胞命运之一来确定分叉点。使用30通道的mass cytometry数据，我们显示了Wishbone准确地恢复了小鼠胸腺T细胞发育的已知阶段，包括分叉点。我们还将该算法应用于小鼠骨髓分化，并证明了其对其他谱系的推广。Wishbone与diffusion maps，SCUBA和Monocle的比较表明，它在排序细胞和识别分支点方面均优于这些方法。
1、引言
多细胞生物是从经历了许多增殖和分化阶段的单细胞发育而来的，从而导致了各种各样的祖细胞和终末细胞类型。尽管这些过程中的许多关键阶段和细胞群已使用荧光激活的细胞分选和遗传扰动进行了表征，但许多发展仍是未知的。新兴的高通量技术，例如单细胞RNA-Seq和mass cytometry，可以在单细胞中同时测量大量参数，并在不干扰的情况下对整个组织进行检查。由于许多组织通过连续和异步发育维持体内平衡，这提供了以高分辨率在几乎所有成熟阶段测量细胞的机会。面临的挑战是设计一种能够根据细胞的成熟度对细胞进行排序并确定产生功能不同细胞的完整Complement的分支点的算法。
最近，一些报告已经证明了根据单细胞的成熟度对其排序的方法[3,4]。但是，这些方法采用的是非分支轨迹，因此很难模拟多个细胞命运。构建分支轨迹的关键挑战是根据细胞的发育成熟度对其排序，确定分支点以及将细胞与其各自的分支相关联。诸如SCUBA[5]之类的方法可以识别数据中的分支以及单元的伪时间顺序，但是在时间分辨率和准确性上却有相当大的损失。
在这里，我们介绍了分支系统的轨迹检测算法Wishbone。 我们使用大规模流式细胞仪数据来测量小鼠胸腺中T细胞的发育，其中淋巴样祖细胞分化为CD8 +细胞毒性或CD4 +辅助T细胞，以证明Wishbone的准确性和鲁棒性。Wishbone 以高精度和发育分辨率恢复了T细胞发育的已知阶段。我们沿统一的分支轨迹从单细胞快照中排序双阴性（DN）1-4，双阳性（DP），CD4 +和CD8 +细胞。结果表明，与竞争方法相比，Wishbone以更高的准确性和分辨率恢复了T细胞发展的已知阶段，产生的轨迹和分支符合目前流行的T细胞分化模型和完整的细胞类型。
我们认为，标志物表达的异质性在很大程度上是由于发育成熟过程中的重编程而不是表达的随机性。此外，我们将Wishbone应用于大规模流式细胞仪[2]，产生的早期和晚期人类骨髓分化数据，以及使用单细胞RNA-Seq [6] 生成的小鼠骨髓分化数据。Wishbone成功地确定了从头开始的骨髓发育中的成熟和分支点，证明了其广泛适用于多种单细胞技术的分支轨迹的系统。
2、结果
2.1 学习分支的分化轨迹
为了根据数据推断分支轨迹，Wishbone对数据进行了以下假设：（i）分化是一个连续过程，（ii）组织的快照能够表征整个分化过程，（iii）细胞在命运的交叉路口只能选择两条分化轨迹中的一条。保守表型标记表达的上升和下降决定了分化，因此标准距离度量（欧几里得度量）无法充分捕获两个细胞之间的成熟度差异（图1a）。与我们之前发表的非分支轨迹检测算法Wanderlust [3]类似，我们使用最近邻图来捕获发育距离并使用最短路径来识别细胞的初始顺序。图中的每个节点代表一个单元格，edge根据表达谱的相似性将细胞与最相似的细胞相连（图1a）。使用最短路径来计算单元格之间的距离，即通过图的neighbors的一系列short steps，其中紧密相关的细胞之间的每一步都可能表示成熟度上的相似性。
Wishbone使用来自输入“分化早期细胞”的最短路径来构建细胞的初始顺序，随后使用选定的一组细胞（称为“waypoint”）进行顺序进行优化。最后，waypoint之间距离的不一致性用于标识所有细胞的branch point和分支关联。最近邻图的质量对于准确排序至关重要，而噪声的主要来源是“short-circuits”的存在，即成熟度相距较远的细胞之间的虚假边缘[3]。值得注意的是，单个short-circuits足以在发育距离较远的细胞之间重定向所有最短路径，从而导致错误的排序。short-circuits在分支数据集中特别普遍，因为分叉点之后的细胞的表型特征可能不够明显（补充图1）。Wishbone通过在使用扩散图（diffusion maps）[7] 生成的缩小尺寸的投影空间中重构graph来克服这些short-circuits。 扩散图考虑了任何一对细胞之间的所有可能路径，以大大减少短路。Wishbone使用top扩散分量（diffusion components）来构造graph，捕获数据中的主要几何结构，同时消除可能由测量噪声引起的微小波动。
该算法使用一组选定的单元格（称为Waypoints）来充当graph不同区域的引导。Waypoints是随机采样的细胞，选择它们代表整个轨迹及其分支上的区域。每个waypoint根据其与所有其他细胞的计算距离（图1b）提供一个视角。细胞在轨迹中的位置是通过平均所有waypoint点的视角来确定的，越近的waypoint的权重越高。因此，更近，更可靠的waypoint主要决定了细胞的位置，同时保留了远处waypoint的影响程度，以得出一致的全局结构（图1b，右下图）。
waypoint也是识别分支点的关键。如果两个waypoint i和t沿同一轨迹，则从早期细胞到t的最短路径与经过i的路径之间的差值接近于零（Lsi+Lit = Lst）（图1c，左图）。另一方面，如果两个waypoint在不同的分支上，则该差值实质上大于零（Lsj+Ljt > Lst）（图1c，中间图）。在存在真实分支的情况下，两个分支的waypoints之间的分歧会累积，以创建两组在每个集合内一致而在各集合之间不一致的waypoint。这些分歧产生了一个结构化的矩阵（图1c，右图）：主干上的航点与所有waypoints的分歧较低，一个分支上的waypoint与同一分支上的其他waypoints一致，并且与不同分支上的所有waypoints都有高度分歧（在线方法）。可以使用聚类方法识别此结构。
具体而言，从频谱聚类技术来看，此矩阵的第二个特征向量总结了给定waypoint的所有分歧，并为waypoint提供了分支关联的定量度量（图1d，左图）。偏离零的程度是细胞成熟度的函数（图1d，左图）。Wishbone恢复沿其发育轨迹的细胞顺序，找到分支点，并将分支点之后的细胞分配给两个分支之一（图1d）。
2.2 小鼠胸腺流式数据集分析
在小鼠胸腺的T细胞发育过程中，CD4 +辅助性T细胞和CD8 +细胞毒性T细胞从淋巴样祖细胞分叉（图2a）[8,9]。我们将Wishbone应用于小鼠胸腺的流式细胞仪数据，并根据其在T细胞发育中的广泛功能选择了表面标记和转录因子（补充表1）。我们从Black6小鼠收集了五个独立的胸腺的数据。每个样品平均收集了23万个细胞。
仅使用表面标记来计算细胞相似性（补充表1）并将DN细胞群体定义为起点，将Wishbone独立应用于每个胸腺数据。沿结果轨迹的标记趋势如图2b所示。Wishbone准确地恢复了T细胞发育的已知阶段（图2b和补充图2），包括分叉成两个单一的阳性谱系（CD4 +和CD8 +）。 具体而言，轨迹开始于DN阶段（CD4-CD8-），过渡到DP阶段（CD4 + CD8 +），最后分支到两个单阳性（SP）种群。 我们注意到Wishbone正确地排列了DN种群：DN2（CD44 + CD25 +），DN3（CD44-CD25 +）和DN4（CD44-CD25-），即使这些细胞很罕见并且构成胸腺细胞的<1％ 。DN1（CD44 + CD25-）细胞极为罕见，我们在轨迹开始时确实观察到了类似于这些细胞的特征（图2b）。
为了进一步测试Wishbone的准确性，我们评估了学习轨迹时未使用的标记的表达趋势：转录因子Runx1和Bcl11b，以及信号分子Notch1（图2c）。 所有这些标记的丰度与其在T细胞发育的DN阶段中的已知作用和时机一致（补充注释1）。
Wishbone准确性的其他证据是标记在整个轨迹过程中变化的紧密性（图2d）。 不仅中位标记水平遵循预期趋势，而且几乎每个单个单元格都正确放置在轨迹中，如跨大部分轨迹的标记的低方差所示。无论标记是否用于学习轨迹（图2d和补充图3），标记的方差都很小，从而增强了Wishbone结果的鲁棒性。
先前关于胸腺发育的研究主要依靠遗传扰动和随后的细胞分选，而这些分选总是消除特定的发育区室。在同时测量30个通道的情况下，我们可以沿着统一的分叉轨迹将单个胸腺的DN，DP，CD4 +和CD8 +细胞放置在一个轨迹上，并准确无误地从胸腺组织直接测量沿轨迹的多个事件的过程。我们使用导数分析在单个参照系中沿轨迹（图2e）对关键事件进行时间计时，发现Wishbone使用从复杂的原始组织收集的细胞来恢复细胞的精确时间顺序和分支以及高分辨率和准确性（补充说明1和图2b，e）。
2.3 Wishbone结果在所有重复试验中保持一致，并且对参数选择具有鲁棒性
我们研究了轨迹和分支对各种自由参数的敏感性：用于图构建的邻居数k，waypoint数nW，来自细胞的waypoint点采样以及使用的扩散分量数。细胞及其分支点的排序对于跨重复的这些不同参数选择非常鲁棒性（补充图5和6，以及补充说明2）。Wishbone的结果在很大程度上排除了用于学习的单个标记（补充图7）。 此外，无论是否使用DN或SP种群的细胞作为输入早期细胞，由wishbone鉴定的分支都保持一致（补充图8）。
2.4 独立细胞类型内标志物水平的成熟度控制
我们在经典的细胞表面标志物中观察到相当多的异质性。我们假设这种变化的至少一部分可能是发育成熟的结果，其中来自不同发育阶段的细胞被汇集到一个单独的门控种群中。使用Wishbone轨迹的精细时间分辨率，我们将以细胞发育进程为条件的标记差异与在门控细胞群中观察到的差异进行了比较。为了进行比较，我们首先使用标准门控方案在两个谱系标记：CD4和CD8的表达上鉴定了两个SP细胞群（图3b）。 接下来，我们将这些门控群的方差与相应标记的方差进行了比较，其条件取决于wishbone轨迹。 在两个SP种群中，沿轨迹的成熟度受到控制时，Wishebone谱系标记CD4和CD8以及共受体CD3的变异比门控种群的种群变异要低得多（补充图9a）。
作为一项附加测试，我们在学习轨迹的同时没有使用CD3作为标记之一就跑了Wishbone。 所确定的轨迹和分支与获得的包括CD3的结果相似，并且在整个轨迹上仅伴随着CD3方差的小幅增加（图3c）。 但是，与门控种群的变异相比，所有受体和共受体分子CD4，CD8和CD3的变异沿轨迹仍显着较低（图3d）。 当从学习中排除CD4或CD8中的任何一个时，这些结果类似地成立（补充图9b，c）。 总的来说，这些结果表明，门控种群中标记表达异质性的很大一部分是比较沿不同发育阶段的细胞而不是标记表达随机性的结果。
2.5 沿着SP轨迹的转录因子动力学
接下来，我们利用一个包括转录因子ThPOK、Gata3和Runx3在内的修正的panel来探索关键转录因子沿着两个SP轨迹的动态变化。ThPOK和Gata3已被证明是CD4 SP的关键，而Runx3被证明是CD8 SP的关键。这些因素沿轨迹的动态如图4a，b所示。
为了将这些转录因子的动态联系起来，我们使用CD69和CD62L来识别成熟的标志，例如谱系承诺和成功的阴性选择（图4c，d）。我们的研究结果表明，这些因素在实现承诺时遵循不同的动态模式。在阳性选择期间，ThPOK和Gata3上调，但Runx3似乎仅在检测到的分支点之后上调（图4c、d、补充注释3和补充图10）。Gata3被证明可以调节ThPOK的表达11，并且可以解释为什么ThPOK在表达变化上落后于Gata3（图4c）。不同的动态可能预示着不同的调控机制，通过这些机制，这些因素实现了血统承诺。需要进一步的实验来阐明这些机制。
我们比较了沿着Wishbone轨迹的标记动力学和从发育中的SP细胞的选通得到的动力学，并比较了沿着Wishbone轨迹每个群体中细胞的顺序。
我们观察到，大多数门控群体中的细胞沿着轨迹分布（图4e和补充图12b），特别是CD4+CD8int（中间型为int）群体的细胞，在这里谱系决定被认为发生15。为了解决Wishbone和gating之间的差异，我们根据其在Wishbone轨迹中的位置将CD4+CD8int细胞分为“早期”和“晚期”两组（图4e，在线方法），并比较两组中已知成熟标记物CD69、CD24和CD62L的表达。与“晚期”细胞相比，“早期”细胞表现出明显的高表达CD69和CD24，CD62L的表达较低（图4f；P<1×10−6，Kolmogorov-Smirnov检验），表明“早期”和“晚期”的细胞分别是不成熟和成熟的。对于CD4和CD8分支中的附加门控群体的类似结果（补充图12）表明，传统门控方案导致在每个门中包含不同成熟阶段的细胞。我们的结论是，Wishbone提供了更可靠的细胞成熟度估计，因此标记动力学沿着SP成熟。
为了了解这种差异的来源，我们比较了免疫基因组计划（ImmGen）分类人群中标记物的平均基因表达量16和Wishbone标记物动力学，并且我们的门控群体中的平均蛋白表达量和ImmGen群体中的平均mRNA表达量是相关的。因此，沿叉形杆轨迹观察到的动力学与选通之间的差异不是特定于数据集的观察（图4g）。在每个门中发育不同的细胞的混合可能会导致在成熟过程中表达改变模式的混杂效应。尽管CD24沿着叉骨轨迹呈持续下降趋势，但在门控和ImmGen群体中表达的变化更大（图4g，左图）。在门控群体中，由于CD69本身用于门控（图4g，右图），所以在门控细胞群中没有观察到CD69的持续上调。最后，尽管Gata3没有用于选通（补充图12g），但Gata3在选通群体中的表达变化并未显示在Wishbone中观察到的动态（图4c（1，5）），进一步证明了Wishbone在更高分辨率下恢复标记动态的能力。
2.6 Wishbone在人髓系发育中的应用
我们在两个人类髓系发育数据集上评估了Wishbone的性能。我们使用了先前发表的流式细胞术数据[2]，这些数据包括适合于恢复早期和晚期髓样分叉的标记物[17]，但不适用于过渡性髓样细胞群的精细分析。
Wishbone能够跟踪造血干细胞和祖细胞（hematopoietic stem and progenitor cells, HSPCs）的单核细胞（CD14+CD11b+CD11c+）和红系细胞（CD235ab+）的分化（图5a和补充图14a），并精确地恢复单核细胞和红系分支（图5a）。此外，沿着轨迹的标记的表达与已知文献18一致（补充图14c）。
Wishbone精确地恢复了从HSPCs开始的轨迹，以及两类单核细胞的分支，经典单核细胞和CD16单核细胞（CD16+CD15+）（图5b）。这是一个更困难的问题，因为除了特征标记CD15和CD16（补充图14d），大多数标记在两个群体之间显示相同的分布。标记沿着轨迹的表达与已知文献一致（补充图14e），检测到的分叉点与CD38的显著下调一致。
2.7 Wishbone扩展到单细胞RNA序列数据
单细胞RNA序列技术可以描绘数千个单细胞，并能够对正在发育的系统进行全基因组鉴定。然而，这些数据带来了挑战，因为许多基因的行为与发育成熟度无关，而是与细胞周期和压力等过程有关。因此，轨迹和分支检测的成功依赖于去除不相关的因素，并保留那些具有差异性的跟踪因素。
我们使用最近发表的单细胞RNA序列数据，从HSPC中选择参与髓系和红系祖细胞分化的细胞（图5c和补充图15a）。我们设计了一个Wishbone的扩展，适用于单细胞RNA序列，它使用扩散图来帮助关注与发育和成熟相关的成分。扩散图捕捉数据中的主要结构和趋势，在流式细胞仪的情况下，不同的扩散成分跟踪组成细胞类型之间的差异（补充图6a）。我们将基因投射到每个扩散组分上，根据基因在这个组分上的表达情况对基因进行排序，然后用这个排序进行基因集富集分析（GSEA）。一些扩散组分被富集用于免疫相关功能（例如，防御反应、抗原处理和吞噬作用），而其他组分被富集用于其他生物过程（例如，细胞周期、核糖体生物发生和代谢过程）（补充图15c，在线方法）。这为保留与差异化过程最相关的组件提供了一种自然的方法。基于类似的推理，Buettner等人21使用潜在变量模型来消除单细胞RNA序列中细胞周期的贡献。
我们构建了Wishbone的邻域图，该图只基于将数据投影到与分化相关的组分上，一旦该图被构建，我们就按照流式细胞仪的方法继续进行Wishbone的操作。Wishbone精确地恢复了从HSPCs开始到两个前体细胞类型和分支联系的轨迹（图5c）。标记趋势显示HSPC标记物CD34沿着轨迹持续减少，髓系标记物Mpo22沿髓系分支的表达增加（补充图15c）。与已知生物学一致，Gata2在Gata1之前沿着红系23上调（图5d）。
2.7 Wishbone在轨迹识别和分支识别方面都优于竞争方法
我们比较了Wishbone与Diffusion maps[7]，SCUBA[5]和Monocle[24]的性能（图6）。尽管我们使用扩散图构建了kNN图，但我们测试了扩散图是否可以概括发展轨迹。请注意，扩散贴图并未明确提供分叉，我们只能评估其概括准确排序的能力。扩散图正确地恢复了T细胞发育的各个已知阶段（补充图16b），尤其是在早期的DN，但是在DP和SP群体中却缺乏足够的分辨率（补充图16a，b）。 此外，尽管扩散图在两个髓样数据集（补充图16c，e）中恢复了正确的顺序，但在单核细胞数据集中，扩散图对成熟细胞之后的前体进行了排序（补充图16d）。 因此，尽管扩散图大大减少了数据中的噪声（补充图1，在线方法），但Wishbone采取的细化细胞顺序的额外步骤对于得出鲁棒的高分辨率轨迹至关重要。
接下来，我们将Wishbone与SCUBA进行了比较。SCUBA具有较大的内存占用空间，因此只能通过从胸腺数据集中对20,000个单元进行二次采样来运行。胸腺的SCUBA轨迹未正确排序各阶段，我们观察到了不同的DN细胞散布在DP细胞之间（图6a）。 SCUBA确实将两个SP种群确定为两个分支，但与Wishbone相比，在分叉点的分辨率降低了（图6a）。而且，不同的随机细胞样本在轨迹和分支方面都导致不一致的结果（补充图17a，c）。大量细胞计数单核细胞和单细胞RNA-Seq髓系数据集中的SCUBA轨迹与已知生物学一致，但产生了大量不连贯的分支（图6c和补充图17d）。 此外，SCUBA无法正确恢复单核-红系数据集中的顺序和分支（补充图17e）。
最后，我们将Wishbone与Monocle24进行了比较，Monocle是专为应用于单细胞RNA-Seq数据而开发的。Monocle不能运行超过1,000个单元，因此我们从每个数据集中对1,000个单元进行了二次采样。 在散布着DN和DP细胞的情况下，Monocle无法恢复胸腺数据中的正确顺序（图6b）。 尽管轨迹确实终止于两个SP种群，但Monocle鉴定出的分支并不对应于T细胞发育的任何特定阶段，并且两个SP种群都被确定为同一分支的一部分（图6b）。 重复的数据二次抽样导致两个SP种群被重复分配到同一分支的结果大体上不一致（补充图18a–c）。Monocle还无法通过错误的细胞排序和缺乏连贯分支的检测来恢复单细胞RNA-Seq髓系数据集中的轨迹和分支（图6d）。单核细胞确实恢复了单核细胞数据集中的排序，但所有髓样数据集中的分支结果均不对应于正确的成熟细胞群（补充图18d，e）。
因此，Wishbone在细胞的精细排序，分支点和分支关联的识别以及复制间一致的鲁棒性方面优于竞争方法。
3.讨论
我们已经开发了一种算法，该算法能够沿着分支的发展轨迹对细胞进行精确且高分辨率的排序（补充图19）。 我们首先使用大规模流式细胞仪的通量收集每个样品中≥200,000个细胞，证明了Wishbone对小鼠胸腺T细胞发育的影响。Wishbone构造了从DN阶段到两个SP谱系成熟的分叉轨迹，提供了事件的顺序和时机，这些事件和事件的时间紧紧概述了该系统的先前研究15。 Wishbone的高分辨率使我们能够识别谱系标记的细微但关键的动力学，例如在DN向DP细胞过渡期间检测稀有CD8 + CD4-中间SP细胞以及向DP末端过渡CD4 + / CD8low中间状态。
选择一个好的标记集是我们获得分辨率的关键。 可以通过先验知识和初步筛选的组合来指导标记的选择。 但是，在骨髓分支中，我们证明了即使只有有限的面板（仅包含少量可区分的骨髓标志物），Wishbone仍能正确排序细胞，识别分叉以及将细胞关联至正确的分支。 尽管不一定可以提供明确的事实依据，但SCUBA和Monocle都无法在这些更具挑战性的数据集中恢复与已知生物学一致的表达趋势和分支。Wishbone只需要几个规范标记就可以正确地识别分叉，并且由于包含了其他标记，因此在过渡人群中获得了越来越精细的分辨率。
单细胞RNA-Seq是大规模流式细胞术的一种有吸引力的替代方法，因为其无偏见的，全基因组的性质提供了数千种基因的测量，并且避免了对先验选择有限标记集的需求。 但是，与发育无关的转录变化会混淆分析，甚至与发育相关的基因的数据也具有很大的噪音，包括脱落效应27。 我们使用扩散图巩固了关键的生物学趋势，并删除了不相关的生物学过程。 我们证明，横臂骨大大优于专门为单细胞RNA-Seq数据开发的方法5。
即使单细胞RNA-Seq的吞吐率不断提高，目前的数据集仍包含数千个细胞，而大规模细胞计数法中却有数十万个细胞。 由于已显示出过渡种群与1 / 10,000个细胞一样稀少3，因此大规模细胞计数的通量更适合于实现更精细的时间分辨率。 我们认为，这两种技术是互补的。 例如，单细胞RNA-Seq可以用于研究较少的发育系统中的无偏标记选择，然后可以使用已鉴定的专家组通过大规模细胞术获得更精细的时间分辨率。
