Basic Information 英文标题:Cell2fate infers RNA velocity modules to improve cell fate prediction 中文标题:Cell2fate推断RNA速度模块以改进细胞命运预测 发表日期:03 March 2025 文章类型:Article 所属期刊:Nature Methods 文章作者:Alexander Aivazidis | Omer Ali Bayraktar 文章链接:https://www.nature.com/articles/s41592-025-02608-3 Abstract Para_01 RNA velocity利用包含在拼接和未拼接RNA计数中的时间信息来推断转录动力学。 现有的速度模型通常依赖于粗略的生物物理简化或数值近似来解决潜在的常微分方程(ODEs),这可能在具有挑战性的环境中,如细胞轨迹中复杂或微弱的转录率变化,会损害准确性。 在这里,我们提出了cell2fate,一种基于速度ODE线性化的RNA速度公式,这允许以完全贝叶斯的方式求解一个在生物物理学上更准确的模型。 因此,cell2fate将RNA速度解决方案分解为模块,提供了RNA速度与统计降维之间的生物物理联系。 我们在现实世界设置中全面评估了cell2fate,证明其提高了解释能力和重建复杂动态和稀有及成熟细胞类型中微弱动态信号的能力。 最后,我们将cell2fate应用于发育中的人脑,在此过程中,我们将RNA速度模块空间映射到组织结构上,连接了组织的空间组织与转录的时间动态。 Main Para_01 ‘RNA速度’这一概念涉及从单细胞RNA测序(scRNA-seq)中的剪接和未剪接计数推断转录动力学,显示出显著的潜力。 RNA速度模型的第一批实现经历了概念和技术上的改进,包括参数推断的改进以及使用数值方法来解决基础微分方程。 然而,这些现有的改进存在取舍,要么引入粗略的生物物理近似,要么依赖大量的数值近似。 因此,根本的挑战仍然是定义一个数学上合理的框架,既能捕捉现实的转录动力学,同时保持计算和数值的可行性。 Para_02 为了解决上述限制,我们提出了cell2fate,这是一种基于更现实的生物物理转录动力学模型的全贝叶斯RNA速度模型。 Cell2fate使用线性化将描述复杂转录模式(如转录增强)的微分方程分解为可解析的组件。 通过这样做,该模型同时具有表现力、可解释性和计算效率。 将速度问题分解为组件的方法还提供了RNA速度与使用生物物理解决方案进行降维之间的联系。 Para_03 我们评估并在不同的现实场景中基准测试 cell2fate,展示了它能够捕捉稀有和成熟细胞类型中的复杂动态和微弱动力学信号。 最后,我们展示了 cell2fate 如何与空间转录组学结合,从而将转录动态与其空间组织环境联系起来。 Results The cell2fate model cell2fate 模型
Para_01 Cell2fate 建立在已确立的概念基础上,用于 RNA 速度分析1,2,采用一种动态模型来解释单个基因和细胞的剪接(s)和未剪接(u)读取计数的变化(图1a),这可以通过两个耦合的常微分方程来定义: Fig. 1: Cell2fate model overview.
- 图片说明
◉ Cell2fate 允许推断复杂的和微妙的转录动力学(a)◉ 通过建模基因特异性的转录率(b)◉ 使用更少数量的独立模块,这些模块具有简单的动态特性,同时也在 RNA 速度(c)和计数(d)中产生模块化结构。◉ λ 表示模块激活(ON)或失活(OFF)的速率。◉ 上标 'biol' 表示计数(u 和 s)不包括诸如测序深度等技术变异因素。
Para_02 这里,α、β、γ分别表示不同基因g的转录率、剪接率和降解率。 通过求解u和s的常微分方程(ODEs)并将方程拟合到观察到的计数,可以估计未知的速率参数,然后将这些参数代入方程(2),得到每个细胞中剪接计数的变化率,即(\frac{{ds"}}{{dt"}}),这进而产生了通常所说的‘RNA速度’(参考文献1)。 RNA速度模型的一个重要挑战是,作为分化时间函数的转录率(({\alpha }{g"}(t)))可能是复杂且非线性的,反映了对核内活跃转录因子(TF)丰度的隐含依赖性13,然而转录率函数({\alpha } {g"}(t))的积分需要保持可处理性,以便有效地求解常微分方程。 因此,现有的方法要么假设({\alpha }_{g"}(t))为简化的阶跃函数(参考文献1、2、5、7),要么采用数值近似来求解常微分方程(补充说明和补充图3)。 Para_03 在cell2fate中,我们使用了转录率导数的一个扩展,用可单独积分的基本函数表示,我们称之为以下的线性化: Para_04 我们将每个基函数称为一个模块,并用下标m表示。下标i表示模块的状态(开或关)。 ({\hat{\alpha }}_{{mgi"}}) 是当模块处于开状态时(图1b,右下角),该模块的目标转录率,对于所有基因均取非零值。 ({\lambda }_{{mi"}}) 是达到目标转录率的速度,状态i取决于细胞特异性时间尺度Tc上的开关时间Tm,ON/Tm,OFF(图1b,上部)。 这种基函数的选择使得每个单独的常微分方程以及它们的总和都可以解析求解(方法和补充说明)。 参数Tm,ON/Tm,OFF、({\lambda }{{mi"}}) 和Tc在所有基因中是共享的,这大大减少了与现有模型相比需要估计的参数数量,但仍提供了基因特异性的转录率({\alpha } {{mg"}})。 Para_05 除了计算原因之外,线性化还提供了RNA速度和统计降维之间的生物物理联系。 当将线性化视为混合隶属模型时,这种联系变得明显,在该模型中,转录率、RNA速度以及每个基因的剪接和未剪接计数由M个模块的线性组合控制(方法和图1b)。 然后,混合系数可以类比于因子分析或主成分分析中的基因载荷进行解释。 从机制上讲,模块可以被解释为近似由所有活跃的调控蛋白引起的转录速率变化,这些变化作为一组独立的影响,每个影响在由Tm,ON/Tm,OFF定义的时间窗口内有效。 Para_06 Cell2fate 是一个完全贝叶斯模型,它将原始细胞水平计数作为输入,并包括一系列改进以考虑技术变异来源,包括过离散、拼接和未拼接 RNA 分子的检测敏感度变化、环境 RNA 和已知批次(方法)。 分层先验结构实现了对模型参数的有效正则化,同时在基因之间共享证据强度(方法)。 该模型使用概率编程语言 Pyro14 实现,并使用 Pyro 的随机变分推理框架进行拟合,采用定制的变分分布结构来考虑模型参数之间的依赖关系(方法)。 软件实现附带了指导方针和启发式方法来确定超参数,例如模块的数量(方法),并且基于 scvi-tools15 进行构建,以促进其在现有工作流程中的集成。 Improved predictions of cell fates and transcription rates 改进的细胞命运预测和转录率预测
Para_01 为了评估cell2fate的性能,我们通过评估估计的细胞命运轨迹与先验知识的一致性,将其与现有的RNA速度方法进行了比较。 简而言之,我们考虑了交叉边界方向正确性(CBDir)指标来评估替代模型,从而对已知转换的细胞簇边界处的转移概率一致性进行评分(方法)。 Para_02 我们考虑了涵盖不同模型类别和参数推断方法的十种RNA速度方法(表1和方法章节)。 所有方法都被应用于五个单细胞RNA测序数据集,包括广泛使用的基准数据集,如发育中的小鼠齿状回16和胰腺17(补充图1)。 为了测试不同方法解析复杂转录动力学的能力,我们还检查了小鼠红系成熟18和人类骨髓19:两个在细胞轨迹中具有多个转录爆发的数据集12。 最后,我们考虑了一个小鼠骨髓数据集,该数据集的独特分子标识符(UMI)计数明显较低1,从而评估了这些模型在低覆盖率数据中的表现程度。 Table 1 Summary of RNA velocity model versions and parameter settings used in benchmark 表1 基准测试中使用的RNA速度模型版本和参数设置摘要
Para_03 平均而言,在全部五个数据集中,cell2fate 获得了最佳的 CBDir 分数。 更重要的是,cell2fate 在所有数据集中推断出了细胞命运转换的正确方向性,而所有其他方法(pyroVelocity_model2 除外)在至少一个基准设置中推断出反向动态(对应于负的 CBDir 值)。 检查基准测试结果,我们可以将其性能归因于克服了下面详细说明的两个主要挑战。 Para_04 首先,cell2fate 提供了足够的统计功效来从微妙的转录动力学中识别正确的速度流。 例如,在小鼠齿状回数据集中,其他方法未能解决颗粒神经元后期成熟轨迹的问题,从而错误地推断出成熟细胞转变为未成熟细胞(图2b,蓝色内嵌框和补充图2)。 当将 CellRank 应用于速度估计而不是一致流形逼近和投影(UMAP)时,得到了相同的结果,这表明只有 cell2fate 能够解决朝向成熟颗粒神经元的轨迹问题(补充图3和4)。 Fig. 2: Enhanced performance of cell2fate in benchmark of ten methods across five datasets.
- 图片说明
◉ 十种方法在五个数据集上重建已知轨迹的表现。显示的是CBDir3,正值对应正确的谱系重建。◉ b-d,展示了一些具有特定挑战的数据集的例子,包括成熟细胞类型的弱信号(OPC,少突胶质细胞前体细胞;OL,少突胶质细胞)(b)和复杂的转录率动态(c,d)。◉ d,推断出的两个选定基因的转录率,据推测这两个基因的转录率有逐步变化18。◉ e,来自cell2fate的细胞特异性时间估计。左边,星形胶质细胞的时间远高于神经元谱系。右边,一个连续的时间范围表示单一谱系。◉ f,后验时间的标准差可以作为不确定性度量,用于评估数据集对cell2fate分析的适用性。在一个外周血单核细胞的稳态数据集中,这个标准差在整个过程中接近于1。◉ NK,自然杀伤细胞;Treg,调节性T细胞。
Para_05 其次,cell2fate 正确地重构了复杂的转录动力学。 在小鼠红系成熟的数据集中,该模型解析了正确的细胞轨迹,而其他模型往往表现不佳(图 2c 和补充图 5)。 先前基于对拼接和未拼接计数的视觉检查,通过手动注释的细胞簇进行分析18,提供了证据表明小鼠红系谱系形成包含许多‘多速率动力学’基因,如 Hba-x 和 Nudt4,这些基因在整个细胞成熟轨迹中显示出协调的转录速率变化18。 一致的是,cell2fate 重现了这些多速率动力学基因中的逐步转录速率提升18(图 2d,青色线)。 相比之下,其他方法,例如 pyroVelocity_model2,只能预测一个非零的转录速率,因为它们的基础动力学模型更简单(图 2d,绿色线),我们再次使用 CellRank 证实了这一结果(补充图 6 和 7)。 先前的分析18基于对拼接和未拼接计数的视觉检查,通过手动注释的细胞簇进行了分析,提供了证据表明小鼠红系谱系形成包含许多‘多速率动力学’基因,如 Hba-x 和 Nudt4,这些基因在整个细胞成熟轨迹中显示出协调的转录速率变化18。 Para_06 作为 CBDir 的替代指标,我们还评估了聚类之间估计时间差异与先验知识的一致性(补充表 8),并将时间输出与不同小鼠样本已知的发育年龄进行了比较(补充图 8 和 9);这些替代指标和基准验证了 cell2fate 的性能优势。 值得注意的是,具有较高 cell2fate 估计时间的颗粒神经元亚簇包含更多来自后期发育阶段的细胞(补充图 8b),验证了 cell2fate 的估计。 替代方法以及标准扩散伪时间分析未能解析这一颗粒神经元成熟轨迹(补充图 8c 和 10)。 此外,cell2fate 对先验分布假设的变化具有鲁棒性(补充图 11),并且我们确认了当将模型应用于半合成数据时,该模型能够估算真实的动态速率(补充图 12)。 最后,我们注意到尽管 cell2fate 的计算需求高于一些现有方法,但该模型的计算需求与替代方案一致,因此 cell2fate 可以方便地应用于更大的数据集(补充表 6 和 7 以及补充图 13)。 Para_07 我们注意到cell2fate的细胞特异性时间尺度(图1a、b)提供了两个额外的用例,可以帮助获得更深入的见解。 首先,它有助于识别细胞谱系进展和不同的细胞谱系。 例如,在小鼠齿状回数据集中,颗粒神经元和星形胶质细胞被分配了显著不同的时间点,而少突胶质细胞则占据了中间的时间段(图2e,左),这与这三个细胞类型的谱系起源不同一致16。 相比之下,在小鼠红系成熟数据集中,识别出一个具有单一连接时间范围的单一谱系(图2e,右)。 其次,检查细胞特异性时间的贝叶斯后验可以为跨数据集和数据集内的RNA速度值提供一种原则性的置信度测量。 在上述两个数据集中,个体细胞时间后验分布的变异系数(CV)始终被估计为接近零,表明不确定性较低(图2e,底部)。 相比之下,应用于外周血单核细胞稳态数据集的cell2fate21,其中不期望存在一致的转录动力学12,导致置信度估计的CV接近1,表明高不确定性(图2f)。 后验不确定性还可以用来估计个体转换的置信水平。 为此,cell2fate实现了一个基于从一个簇的细胞特异性时间中的后验样本中有多少高于另一个簇样本的第90百分位数的启发式置信分数。 这个启发式方法与相应转换的CBDir得分相关性良好(ρ = 0.56),表明它是经过良好校准的(补充图14)。 因此,细胞特异性时间的后验可以作为质量控制来评估cell2fate是否在一个给定的数据集中识别出有意义的动力学。 Para_08 总之,我们的基准测试展示了cell2fate增强的统计能力来估算细胞轨迹和解析复杂的转录动态,并且能够量化速度估计中的不确定性。 RNA velocity modules map fine stages of late cell maturation RNA速度模块映射晚期细胞成熟中的精细阶段
Para_01 Cell2fate模块是随着时间顺序激活的基因表达程序。 鉴于它们在转录动力学中的生物物理基础,我们预计RNA速度模块可以比缺乏机械基础的传统降维技术(如矩阵分解或聚类)更细致地表征细胞分化过程中的动态过程。 此外,cell2fate还附带了一系列下游分析和可视化工具,使用户能够探索动态过程并获得生物学见解。 Para_02 为了展示cell2fate工具包,我们考虑了作为基准研究一部分的小鼠大脑单细胞数据集(图2b),该数据集描绘了海马齿状回区域在两个发育阶段的情况16。 除了神经祖细胞早期分化为神经元和星形胶质细胞外,该数据集还涵盖了颗粒神经元晚期成熟轨迹(即,在未成熟神经元阶段之后的后期分化),这是一个关键过程,然而这一过程尚未得到充分理解。 更普遍地讲,尚不清楚这种晚期成熟过程是否会在连续的转录阶段展开。 应用于该数据集的先前RNA速度方法2,3能够区分神经元与星形胶质细胞谱系轨迹;然而,最成熟的颗粒神经元的正确轨迹无法被确定(图2b)。 Para_03 Cell2fate应用于该数据集揭示了16个不同的RNA速度模块(补充图15),涵盖了所有预期的细胞轨迹,其中主要谱系对应于颗粒神经元的分化和成熟,源自神经中间祖细胞(nIPCs)、神经母细胞和未成熟神经元。 放射状胶质样前体细胞主要承诺为星形胶质细胞(图3a),这一点从cell2fate的时间估计和CellRank命运概率中都可以看出(补充图16)。 我们还观察到,齿状回中的另一个神经元群体——苔藓细胞——被分配到了颗粒神经元轨迹的中期阶段。 尽管苔藓细胞被认为具有不同的谱系起源,但它们的转录发育与颗粒神经元高度相似22。 Fig. 3: Module decomposition of cell2fate resolves final stages of granule neuron maturation.
- 图片说明
◉ a, 牙齿回旋数据的细胞命运速度图嵌入。◉ b, 在颗粒神经元谱系中,由选定模块随时间产生的拼接计数丰度。◉ c, 选定示例细胞命运模块的激活,选定MOFA23因子的权重以及Leiden聚类用于比较。◉ d, 每个细胞中的模块可以根据它们达到最大稳态计数的程度以及它们是否在表达上增加或减少来分类为状态,以便于解释和可视化。◉ e, 定义为每个细胞中由模块产生的拼接计数的激活。◉ f, 晚期颗粒神经元成熟模块的状态。◉ g, 模块标记基因,定义为它们的转录速率中有很大一部分可以用该模块解释的基因。◉ h, 模块TF基因,定义为也是已知转录因子的模块标记基因。
Para_04 为了更深入地探索神经元分化的过程,我们使用拟合的cell2fate模型来估计每个细胞在推断时间内的九个颗粒神经元谱系模块的总拼接转录本丰度(图3d顶部)。 这一分析确定了从放射状胶质细胞到nIPCs、神经母细胞和未成熟神经元的早期分化过程中模块的连续诱导(模块1-3)。 令人惊讶的是,cell2fate还恢复了成熟颗粒神经元中的动态过程,这些过程由六个模块(模块4-9)依次激活并在时间上重叠于成熟颗粒神经元之间,从而精细地解剖了这些细胞晚期成熟的各个转录窗口(图3d)。 该模型还正确识别了未成熟颗粒神经元与成熟颗粒神经元之间的时序间隙(图3d),这与先前的预期一致16。 cell2fate可视化工具通过提供基于估计分化时间的动态见解,补充了t分布随机邻域嵌入或UMAP,并且还可以可视化其他元数据,如细胞类型注释或发育年龄(图3d底部两个面板)。 To explore the dynamics of neuronal differentiation in greater depth, we used the fitted cell2fate model to estimate the total spliced transcript abundance for each of the nine granule neuron lineage modules in individual cells across the inferred time (Fig. 3d, top). Para_05 总拼接计数估计也可以用来可视化RNA速度模块在细胞中的动态,例如,在常规的UMAP图上(图3b)。 可以检查每个细胞中不同模块的激活情况,类似于传统矩阵分解中的因子活性。 我们还将这些模块激活估计与传统的因素分析和聚类方法进行了比较。 简而言之,多组学因素分析(MOFA)23产生了捕获互补变异来源的因素,其活动模式在分化轨迹上时间上更加分散。 具体来说,这些因素没有对晚期颗粒神经元成熟进行分层(图3b和补充图17)。 我们还观察到细胞2命运模块与MOFA因素基因载荷之间的总体相关性较低,特别是在晚期神经元成熟模块4-9之间(补充图18)。 同样地,在不同分辨率下对scRNA-seq数据集进行Leiden聚类识别出了与神经元成熟不一致的聚类(图3b和补充图19)。 总体而言,这些观察表明细胞2命运捕捉了与现有分解方法互补的变异方面,并且非常适合于进行颗粒神经元成熟的精细解剖。 Para_06 除了模块的活动外,动力学可以根据细胞内它们是增加还是减少的表达来进一步分类(图3c和方法部分)。这两个量在图3e和f中展示了颗粒神经元分化的情况。 这种可视化中的额外动态信息显示,例如,模块9尚未达到稳定状态,这意味着颗粒神经元的成熟过程超出了该数据集所涵盖的时间范围。 Para_07 最后,我们研究了RNA速度模块能在多大程度上为颗粒神经元分化晚期提供更深入的见解。 我们将基因按其转录速率被每个模块解释的程度进行排序,并使用排名靠前的基因作为‘模块标记’(图3g和补充图20)。 我们确定了Rbm24、Tafa1(Fam19a1)和Sptb(模块4)以及Pakap(Palm2)、Pdzd2和Usp19(模块5)是未成熟颗粒神经元中开启的标记物(图3g和补充图20)。 相比之下,Fst、Rapgef5和Moxd1(模块6)、Prr16、Kdm5d和Rpa3(模块7)以及Rgs2、Nudt13和1700048O20Rik(模块8)提供了颗粒神经元成熟晚期阶段的标记物(图3g)。 据报道,Fam19a1能抑制神经干细胞的维持并促进分化24,25,这与其在成熟颗粒神经元中的表达模式一致。 Rgs2在神经活动过程中动态表达26,并参与突触可塑性27,这与它在模块8中晚期诱导的现象一致。 除这两种基因外,据我们所知,此处报告的标记基因尚未在颗粒神经元或脑发育中进行过功能研究。 我们进一步将每个模块的前30个标记基因与一组30个阿尔茨海默病风险基因28和250个自闭症谱系障碍风险基因29,30进行了交叉分析,突出了晚期颗粒神经元成熟模块8,在该模块中,与阿尔茨海默病相关的基因Trappc6a和与自闭症相关的基因Kdm6a和Nr4a2均属于标记物(补充表9)。 这些结果表明,即使是在广泛表征的细胞谱系如颗粒神经元中,细胞命运也能带来新的生物学见解。 Para_08 此外,我们可以提取作为‘模块转录因子’的顶级模块标记基因(图3h)。 转向这样的模块转录因子,Zmat4(模块6)和Tfam(模块8)在晚期颗粒神经元成熟阶段富集(图3h)。 据报道,Zmat4在出生后第7天的小鼠听觉皮层中相对于成年小鼠上调31,而Tfam敲除会导致未成熟的神经元表型32。 然而,它们在颗粒神经元分化中的作用迄今尚未研究。 我们还发现,根据ProBound算法33预测,最有可能被前20个模块转录因子结合的假定启动子序列的基因,在前300个模块基因中出现得更频繁,而不是那些最不可能被这些转录因子结合的基因(补充图21)。 这些转录因子提供了晚期颗粒神经元分化的潜在候选调节因子。 Para_09 总的来说,我们的研究结果展示了cell2fate模块分解对于单细胞RNA测序数据集的强大解释性和统计功效,能够精细地剖析细胞过程,并表明晚发育颗粒神经元的成熟是由不同的阶段组成的。 Spatial mapping of RNA velocity modules RNA速度模块的空间映射
Para_01 生物过程通常在组织中按时间顺序进行。 例如,细胞分化和迁移经常耦合在一起,细胞在其分化轨迹中与不同的空间信号微环境相关联。 在这里,我们试图通过将RNA速度模块映射到一个人类大脑发育的新生成的空间转录组学数据集中(图4a),将由cell2fate捕获的时间信息与组织的空间结构联系起来。 Fig. 4: Cell2fate interfaces with cell2location to spatially map the cortical neurogenesis process in human brain development.
- 图片说明
◉ 细胞2命运稳态表达模块作为参考谱供细胞2定位推断模块在整个空间位置上的丰度。◉ 左侧,实验设置图示。中间,放射状胶质细胞(绿色)和中间祖细胞(紫色)首先产生深层神经元(蓝色),然后转而产生ULn(黄色)。◉ 右侧,祖细胞(绿色、紫色)位于皮层内部更深的位置,在那里它们产生新生神经元(黄色),这些神经元分化并迁移到外部。◉ 人类大脑数据的细胞2命运速度图UMAP嵌入。◉ 每个细胞2命运模块在推断时间上产生的剪接计数。◉ 三个模块位置的总结空间图,根据它们达到稳态表达的细胞类型命名。◉ 单个选定模块的状态、标记和位置。DL,深层;UL,上层。
Para_02 我们专注于胎儿人类的大脑皮层,在这个区域兴奋性神经元的成熟过程在空间和时间上遵循一个高度固定的轨迹。 被称为放射状胶质细胞和中间祖细胞的神经前体细胞位于大脑皮层的发生区,在那里它们依次产生不同的神经元亚型,随后这些神经元迁移到皮层板的深层和表层。 在早期妊娠期间,位于深层的神经元(DLn)比位于表层的神经元(ULn)出生得更早;因此到中期妊娠时,DLn相对于ULn更为成熟。 因此,皮层兴奋性神经元的成熟状态与其空间位置紧密相关。 Para_03 为了检查人类大脑皮层的细胞分化轨迹,我们最初对一个中期供体进行了单核RNA测序(snRNA-seq)分析(10x版本3.0)。 然后我们遵循标准的snRNA-seq处理工作流程(方法)来聚类细胞,并使用文献中标记物注释细胞类型。 我们注释了不同的神经前体细胞(放射状胶质细胞和中间前体细胞)以及不同成熟阶段的兴奋性神经元群体(图4c)。 如预期所料,成熟的神经元表达了DLn标记物,而新生和未成熟的神经元则显示ULn标记物的富集表达(补充图22)。 我们也注释了抑制性神经元和胶质细胞类型,但将它们排除在随后的兴奋性神经元轨迹分析之外。 Para_04 我们随后将cell2fate应用于这项人类大脑snRNA-seq数据集,并观察到从神经祖细胞到新生、未成熟和成熟神经元的预期兴奋性神经元分化轨迹(图4c)。 RNA速度模块将神经元轨迹细分为更精细的成熟阶段,确定了在整个未成熟和成熟神经元中激活的七个顺序且时间上重叠的模块(图4d和补充图23)。 虽然这些模块包含了一些DLn和ULn细胞类型的标记基因,但它们也包括了许多在成人人类皮层中广泛表达于所有兴奋性神经元中的基因,如PSMC3、KRR1和BMPER(补充图24和25)。 这表明这些模块部分地识别出了DLn和ULn共同的神经元成熟轨迹。 Para_05 与cell2fate相比,其他RNA速度方法如scVelo无法准确识别成熟神经元中的速度流向(补充图26)。 此外,cell2fate的综合测量模型使我们能够考虑拼接和未拼接计数的不同检测概率,并纠正我们的人类大脑snRNA-seq数据集中的批次效应(补充图27),这对于从观测计数估计真实的转录动力学至关重要(方法部分)。 Discussion Para_01 在这里,我们介绍了cell2fate,一种RNA速度的贝叶斯模型,能够推断复杂变化或稀有和成熟细胞类型中微弱信号下的转录动力学。 cell2fate的核心创新在于一个基于线性化的速度问题公式,这使得使用解析可处理的线性化组件来解决更符合生物物理特性的模型成为可能。 这种公式的一个好处是这些线性组件可以被检查为可解释的RNA速度模块。 这为cell2fate与统计降维方法之间提供了直接的生物物理联系。 我们通过描述颗粒神经元晚期成熟轨迹来展示这一特性,这些轨迹用其他方法难以捕捉。 此外,RNA速度模块可用于定位空间转录组学数据中的分化轨迹。 我们在发育中的人脑中举例说明了这一点,在那里,神经元分化的RNA速度模块显示出高度的空间组织性。 Para_02 尽管 cell2fate 提高了生物物理准确性,该模型仍然做出了简化的假设,例如恒定的降解率和没有随机分支。 然而,cell2fate 中提出的概念是通用的,并且可以引出几种扩展,这些扩展可以在未来进一步提高模型的生物物理准确性,而无需诉诸数值近似。 这包括具有细胞特异性剪接和降解率的 RNA 速度模型、谱系分支点上的随机速率以及不同时间点转录率之间的因果关系,相当于动态基因调控网络(补充说明)。 从长远来看,动力学模型还应该包括基于空间转录组学测量的信号分子的细胞间相互作用的效果。 向这一目标迈出的一步将是结合 RNA 速度模块映射与空间细胞间交互工具,如 NCEM39,这可能识别驱动分化过程特定步骤的假定交互。 Methods Para_01 细胞命运模型的完整描述以及与其他方法的比较可以在补充说明中找到。 Generative model for cell2fate 细胞命运生成模型
Para_01 方程(1)到(3)产生了cell2fate所基于的这些新的RNA速度方程: 错误!!! - 待补充 Para_03 这里,({H"}{{ce"}}) 表示将细胞分配到实验批次的一次性热编码分类,(,{l"} {{cg,j"}}) 描述了跨细胞的基因检测效率(例如,测序深度、读取比对和量化)差异,({s"}_{{eg,j"}}) 模型化了每个批次中每个基因的环境RNA('汤')。然后,这些'测量'期望值被用来参数化负二项分布(NB)观察模型中观察到的原始计数值 Xcgj = (Ucg, Scg) 的均值。 这些'测量'期望值被用来参数化负二项分布(NB)观察模型中观察到的原始计数值 Xcgj = (Ucg, Scg) 的均值。 Para_04 这里,agj是每个基因的NB过离散参数,分别适用于拼接和未拼接的计数。所有参数具有分层先验分布(补充说明1.6),并通过使用Pyro概率编程框架中的随机变分推理进行推断(补充说明1.8)。 Downstream analysis of cell2fate 细胞命运分析的下游分析
Computation of the RNA velocity graph RNA速度图的计算
Para_01 我们遵循了scVelo包中提出的程序,用于从RNA速度估算中计算细胞间转换概率图,通过将速度图平均化处理到100个后验速度样本上,使得具有高后验不确定性的噪声基因速度在估计转换概率时权重更小。 可选地,也可以通过使用我们的方法得到的平均速度估计和scVelo Python包中原函数来精确地遵循原始过程。 Computation of module activation 模块激活的计算
Para_01 我们通过将后验参数值代入描述剪接计数时间演变的方程来计算模块激活,我们将其定义为细胞中模块产生的总剪接计数。 然后,可以通过在x轴上绘制每个细胞的后验时间,在y轴上绘制每个模块计算出的模块激活来绘制随时间变化的模块激活图。 Calculation of normalized module activation and state 归一化模块激活和状态的计算
Para_01 模块标准化激活表示由模块产生的计数除以模块的稳态计数(图3d,灰色线),这是通过将时间设置为无穷大来计算捕获剪接计数随时间演变的方程得出的。 如果细胞时间小于模块的开启时间或模块的标准化激活低于0.05(图3d,灰色),则定义模块状态为关闭;如果其标准化激活高于0.95(图3d,亮绿色),则定义为开启。 否则,模块处于诱导或抑制状态,具体取决于推断出的细胞时间是在开关关闭时间之前还是之后(图3d,浅绿色或橙色)。 Benchmarking of RNA velocity methods RNA速度方法的基准测试
Processing of datasets 数据集处理
Para_01 我们使用了论文中所有结果的3,000个变异最大的基因,并且这些基因至少在所有样本中检测到20次。 此外,我们应用了主要scVelo教程中建议的标准预处理流程(针对小鼠胰腺发育),但不包括Pyro-Velocity和cell2fate,因为它们不需要预处理。 这包括总体计数标准化、对数变换以及在30个主成分PCA空间中的30个最近邻近点计算平均表达量('kNN平滑')。 Application of velocity models 速度模型的应用
Para_01 我们遵循基于小鼠胰腺发育数据集的在线教程进行了所有方法的学习,然后使用相同的分析流程生成了本文中的基准测试结果。 我们在表1中总结了所有方法的参数设置,并在可获得的情况下包括了方法版本。 对于cell2fate,我们保持了所有训练和模型参数的默认值。 Definition of ground truth cell state transitions 真实细胞状态转换的定义
Para_01 对于小鼠骨髓数据集,我们定义了以下真实集群转换:'分裂'到'前体细胞'再到'激活'。 对于其余数据集,我们使用了来自UniTVelo RNA速度研究的真实转换,该研究也使用了这些数据集进行基准测试。 所有真实转换均包含在补充表格1-5中。 Benchmarking metric 基准评估指标
Para_01 我们使用UniTVelo Python包提供的函数计算了CBDir。 以下是来自UniTVelo出版物对这一指标的解释:'CBDir通过给定的真实情况下的边界细胞来衡量从源簇到目标簇转换的正确性。 这里的源簇边界指的是该簇中与目标簇相邻的细胞,反之亦然。 使用边界细胞是因为它们反映了生物体在短时间内的发展,并且CBDir是通过这种方式计算出来的。 Para_02 其中 (C_A) 是目标簇A中的细胞集合,N(c)表示指定细胞c的邻近细胞。(v_c) 和 (x_c) 是表示细胞c计算出的速度和位置的低维向量。因此,(x_{c'}-x_c) 是在短时间内空间上的位移。 ‘这段时间内的。 CellRank analysis CellRank分析
Para_01 对于齿状回数据集,我们使用默认参数运行了CellRank,并计算了如果该方法将其识别为宏观状态的星形胶质细胞、少突胶质细胞和颗粒成熟状态的命运概率。 当我们使用默认参数n_states = (4, 12)时,VeloVI和VeloVAE返回了一个错误:'ValueError: 将12个聚类将拆分复共轭特征值。请求一个更多的聚类或多一个聚类' 因此我们将n_states设置为(4, 11),这成功运行了。 对于cellDancer和VeloVAE,我们也得到了以下错误:'[…] 值不等于1(rtol=1e-3)。尝试减小容差作为'tol = …',指定预处理器作为'preconditioner = …'或使用直接求解器作为'solver = 'direct''如果矩阵较小' 即使遵循错误消息中的步骤以及关于相关GitHub问题的建议,也无法解决此错误。 对于红系分化数据集,我们也使用默认参数运行了CellRank,并计算了向前体1状态和红系3状态的命运概率。 Simulation benchmark 仿真基准
Para_01 首先,我们将cell2fate拟合到齿状回区域的数据集。 接着,我们从cell2fate生成模型中生成数据,保持所有参数在其拟合值上,仅将其中一个选定参数乘以0.25、0.5、1、2或4。 我们以这种方式研究了四个生物学和实验参数:剪接率、降解率、转录物的检测概率以及NB过离散参数。 然后,我们将cell2fate拟合到这些模拟数据,并计算推断的剪接率和降解率值与已知真实情况的相关性(因为RNA速度是剪接率和降解率的函数)。 Time and memory requirements 时间和内存需求
Para_01 我们使用Python的'time'包测量了所有方法的时间需求。时间测量的起点设置在数据预处理开始前,终点设置在获得速度值后。 我们通过每10秒运行一次'htop'(针对基于CPU的方法)或'nvidia-smi'命令(针对基于GPU的方法),持续最多2分钟,并记录最高值来测量内存需求。 Comparison of decomposition methods on the dentate gyrus dataset 齿状回数据集上分解方法的比较
Application of cell2fate cell2fate的应用
Para_01 在齿状回数据的基准测试运行之后,我们应用了上述下游分析方法来计算和绘制模块激活和状态。 , Application of multi-omics factor analysis 多组学因素分析的应用
Para_01 我们使用了至少有20次检测到的计数的3000个变异最大的基因作为输入,这与cell2fate分析相同。 我们将分析限制在参与神经元分化轨迹的集群上(径向胶质细胞样、nIPC、神经母细胞、颗粒细胞未成熟、颗粒细胞成熟)。 我们将剪接和未剪接的计数相加,然后使用相应的scanpy函数对数据矩阵进行归一化、对数转换和缩放。 我们运行MOFA使用高斯似然性和十个因素(这与cell2fate在相关集群中找到的因素数量相同)。 进一步的运行选项是spikeslab权重=True,ard因子=True,ard权重=True。 ProBound algorithm ProBound算法
Para_01 为了生成补充图29,我们应用了ProBound算法。 ProBound可以预测大多数转录因子(TF)与用户提供的DNA序列的结合亲和力。 我们使用refdata-cellranger-arc-mm10-2020-A-2.0.0基因组作为参考。 为了获得参考中每个基因的假设启动子序列,我们提取了每个基因转录起始位点上游450 bp和下游149 bp的DNA序列。 对于每个转录因子,我们根据其假设转录起始位点的预测ProBound结合亲和力对基因进行排名。 然后,我们使用每个转录因子的前十个和后十个结合目标来生成补充图29所示的结果。 Spatial integration of RNA velocity RNA速度的空间整合
Human tissue 人体组织
Para_01 甲醛固定的石蜡包埋(FFPE)人胎儿大脑第二孕期样本块由受医学研究委员会和惠康基金资助的人类发育生物学资源(HDBR,http://www.hdbr.org)提供,获得了相应的母体书面同意以及富勒姆研究伦理委员会(REC参考号18/LO/0822)和纽卡斯尔及诺森伯兰一号研究伦理委员会(REC参考号18/NE/0290)的批准。 HDBR 受英国人体组织管理局(http://www.hta.gov.uk)监管,并且依照相关的人体组织管理局行为规范运作。 Developing human brain single-nucleus library preparation and sequencing 开发人类大脑单核库制备和测序
Para_01 单细胞核从冷冻胎儿脑组织中分离出来,根据已发表的方案42进行。 简而言之,组织在匀浆缓冲液(250 mM蔗糖,25 mM氯化钾,5 mM氯化镁,10 mM Tris缓冲液,pH 8.0,1 M 1,4-二硫苏糖醇,0.1% Triton X-100,1×蛋白酶抑制剂,0.4 U l−1 RNasin Plus RNase抑制剂,0.2 U l−1 SUPERase·In RNase抑制剂)中通过Dounce匀浆处理,并用40 µm的细胞过滤器过滤。 通过27%的Percoll密度离心去除滤液中的碎片。 批次中的所有核在滴定封装前以相等浓度混合,使用10x Chromium Single Cell 3′版本3.1试剂盒。 文库根据制造商的方案(CG000204)生成,并在NovaSeq 6000系统(Illumina)上使用NovaSeq S4流式细胞进行28-8-91循环的单端测序。 Developing human brain Visium library preparation and sequencing 开发人类大脑Visium文库制备和测序
Para_01 FFPE组织切片(5微米)根据10x基因组Visium CytAssist用户指南(CG000520)进行了染色和成像。 对于胎儿大脑组织,使用的染色时间为:苏木精3分钟;蓝化1分钟;伊红1分钟。 在探针杂交和连接后,使用Visium CytAssist仪器将分析物从玻璃载玻片转移到具有42.25平方毫米捕获区域的Visium CytAssist空间基因表达载玻片上。 探针延伸和文库构建遵循标准的FFPE Visium工作流程(CG000495),在仪器外部进行。 文库通过Illumina-HTP NovaSeq 6000系统上的配对末端测序,在SP流细胞中进行配对末端双索引测序(读取1为28个循环;i7为10个循环;i5为10个循环;读取2为90个循环)。 使用Loupe浏览器生成JSON文件,使用Space Ranger管道版本2022.0705.1(10x基因组)和GRCh38-2020-A参考来处理FASTQ文件。 Data availability Para_01 所有单细胞和Visium数据的原始UMI计数和元数据以anndata格式可供下载:https://cell2fate.cog.sanger.ac.uk/browser.html。 人类脑单核和Visium数据的FASTQ文件可在ENA数据库中通过访问编号PRJEB79988获取。 Code availability Para_01 cell2fate 软件包可以在 GitHub 仓库安装:https://github.com/BayraktarLab/cell2fate。 它也可以在 Zenodo40 上获取:https://zenodo.org/records/13883214。 使用此存储库中的笔记本可以重现 cell2fate 方法的结果:https://github.com/AlexanderAivazidis/cell2fate_notebooks。 所有方法的基准测试结果以及稳健性分析和与实际发育年龄的比较都是通过此存储库中的笔记本完成的:https://github.com/AlexanderAivazidis/fate_benchmarking。