Basic Information 英文标题:Mapping cells through time and space with moscot 中文标题:通过时间和空间映射细胞 with moscot 发表日期:22 January 2025 文章类型:Article 所属期刊:Nature 文章作者:Dominik Klein | Fabian J. Theis 文章链接:https://www.nature.com/articles/s41586-024-08453-2 Abstract Para_01 单细胞基因组技术使我们能够在时间和空间维度上对数百万个细胞进行多模态分析。 然而,实验限制阻碍了我们在天然时间动态下和在其天然空间组织微环境中全面测量细胞的能力。 最优传输作为一种强大的工具已经出现,解决了这些限制,并促进了原始细胞环境的恢复1,2,3,4。 然而,大多数最优传输应用无法整合多模态信息或扩展到单细胞图谱。 在这里,我们介绍了multi-omics单细胞最优传输(moscot),这是一种可扩展的框架,用于单细胞基因组学中的最优传输,支持所有应用中的多模态性。 我们展示了moscot能够有效地重构来自小鼠胚胎的170万个细胞在20个时间点上的发育轨迹。 为了说明moscot在空间上的能力,我们将来自小鼠肝脏样本的单细胞分析中的多模态信息映射到空间转录组数据集中,并对小鼠大脑的多个冠状切片进行了对齐。 我们提出了moscot.spatiotemporal,一种利用跨时空维度的基因表达数据来揭示小鼠胚胎发生时空动态的方法。 我们还使用基因表达和染色质可及性配对测量的数据,解决了之前未发表的小鼠时间分辨胰腺发育数据集中的内分泌谱系关系,包括δ和ε细胞。 我们的发现通过人类诱导多能干细胞胰岛细胞分化模型中NEUROD2作为ε祖细胞调节因子的实验验证得到了证实。 Moscot作为开源软件提供,并附有详细的文档。 Main Para_01 单细胞基因组技术提高了我们对细胞分化和组织结构动态变化的理解。 单细胞检测技术如单细胞RNA测序(scRNA-seq)可以高分辨率地描绘单个细胞的分子状态,而空间检测技术则可以恢复它们的空间组织。 然而,这些实验涉及细胞的破坏,并且只能捕获一部分分子信息。 因此,必须重新校准细胞谱图。 Para_02 先前的工作通过使用最优传输(OT)解决了这些问题,最优传输是一个关注于映射和比较概率分布的领域。 OT在描绘细胞重编程过程方面发挥了重要作用,通过增强空间数据与单细胞参考来重构组织结构,以及通过对齐空间转录组学数据来构建生物系统的共同坐标框架(CCFs)。 Para_03 尽管基于最优传输(OT)的方法有望解决单细胞基因组学中的映射问题,但它们的应用面临着三个关键挑战。 首先,基于OT工具的实现主要针对单模态数据。 其次,当前用于单细胞基因组学的OT方法计算成本高昂,即时间复杂度随着细胞数量呈二次方(或对于Gromov-Wasserstein扩展为三次方)增长。 同样,内存需求也呈二次方增长,这限制了它们应用于大规模数据集的能力。 第三,现有的工具基于异构的实现,这使得难以适应或结合新的方法来解决新问题。 Para_04 这里我们介绍了moscot,一个解决映射和对齐问题的计算框架,并展示了它在时间、空间和时空应用中的能力。 Moscot基于三个设计原则来克服当前的局限性。 Moscot支持多模态数据,提高了可扩展性,并统一了以前在时间域和空间域中的单细胞应用。 我们还介绍了一个以前未描述过的时空应用。 一个直观的应用编程接口(API)与更广泛的scverse8生态系统交互,使这些功能得以实现。 Para_05 我们展示了通过研究小鼠胚胎发育期间的170万个细胞来展示moscot的能力。 此外,我们将来自通过测序进行转录组和抗原表位多模态细胞索引的信息(CITE-seq)映射到小鼠肝脏的高分辨率空间读数上,并对小鼠脑样本的大空间转录组部分进行对齐。 与SPATEO9同时,我们引入了时空映射的概念,并利用小鼠胚胎发育的时空图谱展示了其优势。 最后,我们在小鼠胰腺发育过程中共同分析基因表达和染色质可及性,并应用moscot来更好地描绘delta和epsilon细胞的细胞轨迹。 我们识别出可能驱动谱系形成的转录因子(TFs),并通过实验验证NEUROD2作为调控人类内分泌发生过程中体内epsilon细胞形成的TF。 Moscot解锁了多视角图集规模单细胞应用的最优传输(OT),并且它与广泛的文档一起可在https://moscot-tools.org访问。 Moscot is an OT framework for mapping cells Para_01 Moscot将生物学映射和对齐任务转化为最优传输问题并使用一套一致的算法解决它们。 Moscot接受非配对的数据集作为输入;例如,在不同的时间点进行的测量或对应于不同的空间转录组切片,每个切片包含一个或多个分子模态。 Moscot还接受先前的生物学知识,如细胞生长速率,以指导映射过程。 Moscot解决一个最优传输问题,并生成一个耦合矩阵,该矩阵概率性地关联了各个数据集中样本之间的关系。 配备了这个耦合矩阵,Moscot提供了各种特定应用的下游分析功能(图1a和方法部分)。 Fig. 1: Moscot enables efficient multimodal OT across single-cell applications.
- 图片说明
◉ 这是通用单细胞基因组分析的最优传输(OT)管道示意图(从左到右):实验变化(例如时间点和不同的空间切片)导致不同的细胞群体。◉ 以前的生物学知识(例如增殖率和空间排列)通常是可用的,并且应该用来指导映射过程。◉ 最优传输通过最小化位移成本来对齐细胞分布。◉ 学习到的映射促进了各种下游分析机会。◉ Moscot 引入了三个关键创新,从而释放了最优传输的全部潜力。◉ 首先,它支持所有模型中的多模态数据。◉ 其次,它克服了以前的可扩展性限制,使得能够进行图集规模的应用。◉ 第三,Moscot 是一个统一的框架,在生物学问题中具有一致的API,这将促进易用性,并能够以简单的方式扩展到新问题。◉ 面板a和b是使用BioRender创建的(https://www.biorender.com)。
Para_02 Moscot 基于三种最优传输(OT)的概念来适应各种生物学问题。 这些概念在样本如何跨细胞分布关联方面有所不同:Wasserstein型(W型)5 OT 比较具有相同细胞特征的两组细胞;Gromov-Wasserstein型(GW型)6 OT 比较生活在不同空间中的细胞分布;融合Gromov-Wasserstein型(FGW型)11 OT 比较部分共享特征的细胞(方法和补充说明1)。 我们基于先前基于OT的方法假设,跨时间域和空间域映射细胞(方法)。 Para_03 为了在整个框架中支持多模态性,我们利用了共享潜在表示(图1b和方法部分)。 通过将W型、GW型和FGW型概念的计算时间和内存消耗减少几个数量级,我们使moscot适用于图谱规模的数据集(图1b,方法部分和补充说明2)。 具体来说,我们基于最优传输工具(OTT)12,这是一种可扩展的JAX实现,支持即时编译、成本函数的即时评估和GPU加速(方法部分)。 当数据集大小需要时,我们使用了最近的方法创新13,14,15,这些创新将耦合矩阵约束为低秩,从而使得W型、GW型和FGW型概念具有线性的时空复杂度(补充说明3)。 统一的API使moscot易于使用和扩展(补充图1)。 特别是,模块化实现使得使用相似的基础设施解决不同的生物学问题成为可能。 Reconstructing mouse embryogenesis Para_01 建模处于非稳态的生物系统的细胞状态轨迹需要结合时间序列单细胞研究和计算分析来推断跨时间点的细胞分化。 Waddington OT2(WOT)使用W型最优传输解决了这个问题。 然而,WOT仅限于单模态基因表达数据,并且无法扩展到大型数据集。 因此,我们创建了moscot.time。 我们的模型继承了WOT流行的细胞生长速率建模方法,并适用于多模态数据。 此外,它能够处理数百万个细胞的数据,并且像所有在moscot中的轨迹推断方法一样,可以与诸如CellRank 2(参考文献16,17)之类的工具进行接口对接以用于下游分析(方法)。 Para_02 我们探讨了 moscot.time 改进的可扩展性是否转化为对生物系统的更忠实描述。 因此,我们将模型应用于一个已发表的小鼠早期发育图谱,该图谱包含跨越胚胎第3.5天(E.3.5)到第13.5天(E13.5)的20个时间点的近170万个细胞(图2a和方法部分)。 我们首先评估了是否可以使用 WOT2 分析此数据集。 我们选择了包含超过50万个细胞的E11.5-E12.5时间点对,并在不断增加细胞数量的子集上进行了基准测试(图2b、方法部分和补充表1)。 Moscot.time 在两个时间点上的所有275,000个细胞上计算了耦合,而 WOT 在超过75,000个细胞时就耗尽了内存。 当我们包括低秩最优传输近似(low-rank OT approximation)时,一旦每个时间点超过75,000个细胞,这种附加计算耦合的速度比默认的 moscot.time 更快。 Moscot.time 的线性内存复杂度使其能够在笔记本电脑上处理发育图谱,而 WOT 在服务器上失败(图2b、方法部分和扩展数据图1)。 Fig. 2: Moscot faithfully reconstructs atlas-scale developmental trajectories.
- 图片说明
◉ 这是示例鼠胚胎发生图谱,包括20个时间点和170万个细胞。◉ 对于不断增加数量的细胞,从E11.5-E12.5时间点对中抽取子样本(方法和补充表1),我们比较了WOT2与默认的moscot.time和低秩13,14,15 moscot.time(秩2,000)(补充说明3;由于WOT不支持GPU加速,因此在CPU上运行)。◉ 通过发育阶段的胚层和细胞类型转换得分,比较TOME7和moscot.time的准确性(方法和补充表2)。◉ E8.0-E8.25时间点对的均匀流形近似和投影(UMAP)投影,按原始聚类注释着色。◉ 对于图d中E8.0的五种最常见细胞类型(加粗突出显示),moscot.time的增长率估计(顶部)和clTOME的增长率估计(底部)作为直方图(左侧)和在UMAP投影上的表示(右侧)。黑色垂直线表示增长率为一。◉ 使用moscot.time计算E8.25首次心脏场细胞的祖先概率(左侧)与已知驱动基因Tbx5、Nkx2-5和Tnnt2的基因表达水平(右侧;方法和补充表7)。◉ 使用Spearman相关性量化f中的比较。基因颜色与f相同,每个点代表一个细胞,线条表示线性数据拟合。◉ 绝对Spearman相关值分布(n=36个基因(定形内胚层),n=18个(尾芽),n=39个(心脏场),n=106个(胰腺上皮);垂直线对应四分位数,须表示异常值)用于比较moscot.time和clTOME的祖先概率与已知驱动基因表达之间的关系(方法和补充表4)。◉ 面板a是使用BioRender创建的(https://www.biorender.com)。
Para_03 由于 WOT 无法扩展到如此规模的数据集,发展图谱的作者们设计了一种基于 k 近邻匹配的确定性方法,称为哺乳动物胚胎发生轨迹(TOME)。 我们制定了两种度量标准,分别适用于 germ 层和细胞类型(方法和补充表 2)。 对于这两种度量标准,moscot.time 在所有时间点和发育阶段的表现与 TOME 相当,尽管 TOME 是专门为该数据集设计的(图 2c)。 对于低秩近似,两种度量标准的准确性在足够大的秩时收敛于默认的 moscot.time,并且速度更快(扩展数据图 1b、c)。 此外,moscot 映射的性能在秩和嵌入方面具有鲁棒性(补充图 2–4)。 Para_04 我们进一步比较了TOME和moscot.time在细胞生长速率和死亡率上的表现。 由于TOME仅提供聚类级别的映射,我们将原始方法扩展以产生具有细胞级别TOME(clTOME)的细胞级别输出(方法部分)。 使用E8.0-E8.25时间点对,我们使用moscot.time和clTOME映射了细胞(图2d)。 clTOME经常分配远小于一的生长速率,并预测在这个阶段超过19%的群体处于凋亡状态(图2e和补充表3)。 如此高的死亡率对于胚胎发育来说是不现实的,在E7.0之后,经历凋亡的细胞比例通常<10%18。 相比之下,我们可以调整由moscot.time预测的生长速率,使其更符合实际情况且更具细胞类型特异性(图2e和方法部分)。 这些结果推广到了所有包含足够细胞数量的其他时间点(补充图5-7)。 我们也使用体外重编程数据集2中的scanpy计算的细胞周期得分19比较了预测的生长速率,我们预计预测受到细胞取样随机性的影响较小。 与使用clTOME相比,使用moscot.time生成的预测与每个细胞集合的平均生长速率相关性更好(皮尔逊相关系数分别为0.48和0.13;补充图8)。 Para_05 接下来,我们考虑了细胞命运预测模型的可靠性。 我们首先考虑了E8.25第一心区细胞,这是一种起源于脏壁中胚层的细胞群体20。 我们使用了moscot.time和clTOME来计算祖先概率,这量化了E8.0细胞分化为E8.25第一心区细胞的可能性。 我们将祖先概率与E8.0形成第一心区细胞的已知驱动基因的表达进行了比较(图2f,方法和补充表3)。 使用moscot.time,我们一致获得了更高的绝对斯皮尔曼相关系数(图2g),这一结果推广到了我们在早期发育过程中研究的另外三种细胞类型(图2h和补充表4)。 最后,我们展示了映射细胞群而不是单个细胞在种系层和细胞类型评分方面产生了可比的结果,但在E9.5时未能解析出稀有的原始生殖细胞,并且胰腺上皮的驱动基因相关性较低(方法和扩展数据图2)。 Mapping and aligning spatial samples Para_01 空间组学技术能够对数千个细胞在其天然组织环境中进行分析。 这种数据的分析需要能够整合分子层和空间坐标系统之间数据的方法。 最优传输(OT)已经被证明在这方面很有用,特别是novoSpaRc3用于基因表达图谱绘制,以及PASTE4用于空间转录组数据集的对齐。 Moscot实现了这两个应用,并利用了可扩展的实现方式和更高效的算法(方法)。 Para_02 基于图像的空间转录组学数据通常在测量的基因数量上有限(数百到数千个)。 moscot 的映射问题通过使用 FGW 类型的问题学习解离的单细胞谱与它们的空间组织之间的耦合。 这使我们能够结合分子特征中的细胞相似性和细胞的物理距离(方法)。 OT 解决方案促进了基因表达或附加多模态谱向空间坐标的转移(图 3a)。 Fig. 3: Moscot enables multimodal mapping and alignment of spatial transcriptomic data.
- 图片说明
◉ 多模态单细胞参考数据集映射到空间数据集的示意图。◉ 空间对应关系与 moscot 预测准确性相关。根据方法(方法),对 12 个数据集的真实基因表达和插补基因表达之间的中位数斯皮尔曼等级相关系数与空间对应关系进行线性拟合。◉ 带有从 CITE-seq 数据集映射注释的肝切片(扩展数据图 3)。◉ 方块标记了 d-f 中裁剪的图块。◉ Vwf(内皮细胞标记物)和 Axin2(肝细胞和内皮细胞标记物)的测量基因表达。◉ Vwf 用于识别所有定义 CV 和 PV 边界的上皮细胞。◉ Axin2 是 CV 的阳性标记。◉ Adgrg6 和 Gja5 的预测基因表达,这两种已知是 PV 内皮细胞标记物。◉ 叶酸受体 β 的预测蛋白表达,这是一种 Kupffer 细胞标记物(顶部)以及 Kupffer 细胞和内皮细胞的插补细胞类型(底部)。◉ 多个幻灯片切片对齐到一个共同参考样本的过程示意图。◉ 小鼠大脑第 1 张切片的空间部分的图块可视化,按批次着色(左)和通过 Slc17a7 表达着色(右)。◉ 小鼠大脑第 2 张切片的空间部分的图块可视化,按批次着色(左)和通过 Slc17a7 表达着色(右)。◉ 面板 a 和 g 使用 BioRender 制作(https://www.biorender.com)。
Para_03 我们使用最近的基准测试25对moscot与两种最先进的方法Tangram23和gimVI24进行了对比。 我们通过计算空间坐标中未保留基因的相关性来评估映射过程的质量(方法部分)。 在使用各种技术生成的14个数据集中,Moscot始终优于其他方法。 此外,对于每个数据集,我们量化了空间对应性,这是一个衡量基因表达相似性和物理坐标距离之间相关性的指标,如最初提出的那样(方法部分)。 如果附近细胞具有相似的基因表达谱,则空间转录组学数据集具有高空间对应性(补充图9)。 Moscot在空间对应性和准确性之间显示出正相关(图3b),这表明它可以利用基因表达距离和物理空间之间的空间关联。 然而,即使空间对应性较低,Moscot也优于其他方法。 Para_04 我们随后着手将多模态单细胞谱图映射到它们的空间背景中。 这种方法特别值得关注,因为空间转录组技术大多仅限于基因表达测量。 我们考虑了一个大约包含91,000个细胞的小鼠肝脏CITE-seq数据集和一个由大约367,000个使用Vizgen MERSCOPE平台测量的细胞组成的空间转录组切片(图3c)。 我们将基因表达、蛋白质和空间信息结合起来以恢复蛋白质的空间组织(方法部分)。 然后,我们将CITE-seq数据集中的注释进行映射,因为在原始数据中没有提供细胞类型注释(扩展数据图3)。 由于时间或内存复杂性过高,无法使用其他任何方法。 Para_05 肝脏生理学中的一个核心问题是对中央静脉(CVs)和门静脉(PVs)的识别,以表征肝脏分区27。 这个问题可以通过考虑标记基因的表达模式、细胞类型定位和蛋白质丰度来解决。 中央静脉可以通过Axin2来识别,Axin2是一种与中央静脉相关的内皮细胞标记物28(图3d)。 同样地,Vwf作为一种已知的血管内皮细胞标记物,表明了中央静脉和门静脉的存在29。 然而,由于空间转录组数据中测量的基因数量有限,基于标记基因表达来识别门静脉变得具有挑战性。 利用moscot,我们通过映射特定于门静脉的标记物Adgrg6和Gja5(参考文献26)的表达(图3e和补充图10),克服了这一限制。 表征肝脏分区的细胞生态位的另一个局限是缺乏详细的细胞类型注释和蛋白质表达。 因此,我们使用moscot将单细胞数据集提供的细胞类型注释进行转移。 重点关注称为库普弗细胞的定居肝巨噬细胞,我们确认它们在门静脉周围区域的富集存在,在那里肝脏窦状隙更为普遍26。 我们通过将叶酸受体β蛋白映射到其空间组织,证实了我们的发现(图3f)。 通过整合细胞类型注释、测量和推断的标记基因以及转移的蛋白质表达结果,我们可以详细表征小鼠肝脏样本中的组织微环境。 我们通过在人类扁桃体的空间多组学数据集中推断转座酶可及染色质测序(ATAC–seq)数据,定量确认了整合多种模态的好处。 这项分析包括空间解析的ATAC–seq和RNA-seq数据的联合分析(补充图11)。 Para_06 一个不同的普遍任务是构建感兴趣组织区域的共识视图。 这需要将来自连续切片或不同生物重复样本同一切片的多个空间测量进行对齐。 moscot 的对齐问题促进了多个切片的空间对齐,并从多个空间转录组切片中构建这种共识视图(图 3g)。 这是构建生物系统共同坐标框架的重要步骤。 首先,我们评估了 moscot 在空间上对齐从前人基准研究中改编的合成数据集的能力,以及与其他不特定于空间组学数据的注册方法相比的能力。 基准结果显示,moscot 的表现与专门针对空间组学数据的方法 PASTE4 相当或更好(方法和补充图 12 和 13)。 Para_07 接下来,我们着手研究这些方法在更大数据集上的可扩展性。 为此,我们使用了MERSCOPE(方法)中的大脑冠状切片数据集。 对于像PASTE这样的方法来说,这个数据集大得令人望而却步(第一部分约250,000个细胞,第二部分约300,000个细胞;方法)。 Moscot准确地将两个样本对齐到小鼠大脑两个冠状切片的参考滑片。 我们观察到,对于大多数基因,基因表达密度在细胞邻域之间存在强烈的一致性,无论是定量上还是视觉上(图3h,i,扩展数据图4和补充图14和15)。 Charting spatiotemporal mouse development Para_01 发育系统空间分辨的单细胞数据集的出现带来了开发能够描绘细胞轨迹并利用细胞表型上内在和外在影响的方法的挑战。 这里我们介绍一种轨迹推断方法,该方法结合了基因表达谱相似性和物理距离来推断更准确的轨迹。 该方法基于一种FGW类型的问题,将moscot.time和moscot.space合并为一种时空方法(方法部分)。 Para_02 我们评估了 moscot 在小鼠胚胎时空转录组图谱(MOSTA)10上进行轨迹推断的能力,该图谱包括从E9.5到E16.5的八个时间点。 我们分析了每个时间点的一个切片,总共得到了大约50万个空间数组位置(此后记为bins,根据先前描述的表示法10;图4a和方法部分)。 我们使用了作者提供的主要组织区域和器官的注释10,并评估了计算出的轨迹上的注释转换得分(方法部分和补充表5)。 我们将 moscot.spatiotemporal 的性能与仅使用时间点间基因表达信息计算的轨迹进行了比较,使用的方法是 moscot.time(一个W型问题)或TOME7(图2)。 考虑到轨迹推断中的空间相似性,注释转换得分的预测得到了改善,在各个时间点上的平均改进幅度分别为5%和13%,相对于 moscot.time 和 TOME(图4b和方法部分)。 此外,moscot.spatiotemporal 的表现优于 PASTE2(参考文献31),并且对超参数具有鲁棒性(补充图16)。 接下来,我们使用 moscot 来识别肝脏发育的驱动基因和目标基因(方法部分),这揭示了已知的肝基因Afp、Alb和Apoa2以及编码TF HNF4A的已知驱动基因(补充表6)。 Fig. 4: Inference of spatiotemporal dynamics with moscot.
- 图片说明
◉ a, 小鼠胚胎发生时空轨迹推断的示意图。◉ b, 经过整理的跨越发育阶段的转换准确性(方法和补充表5)对于moscot的时间和时空应用与TOME7相比。◉ c, 跨时间点的心脏细胞映射(底部)以及心脏谱系的真实注释(顶部)。◉ d, 通过将moscot与CellRank 2结合找到的心脏谱系驱动基因(参考文献16、17)。顶部,Tbx20编码一种已知在心血管发育中有多种基本作用的转录因子。底部,Myl7编码一种与代谢和心脏再生相关的蛋白质(补充表7)。◉ e, 将仅在最新时间点(E16.5)提供的高分辨率细胞类型注释转移到较早的时间点。◉ f, 基因表达与神经元(x轴)和成纤维细胞(y轴)命运概率之间的皮尔逊相关性。注释的基因是前20个驱动基因之一,并且以前与成纤维细胞和神经元谱系相关联(补充表7和8)。◉ g, 样品神经元驱动基因Neurod2和Sox11的空间可视化(补充表8)。小脑颗粒神经母细胞;皮质;CR,Cajal-Retzius细胞;成纤维细胞;前脑;兴奋性;神经元;后脑;下丘脑;中脑;腹侧下丘脑。面板a使用BioRender创建(https://www.biorender.com)。◉ Panel a was created using BioRender (https://www.biorender.com).
Para_03 随后,我们关注了发育中小鼠胚胎的心脏和脑部区域的命运。对于每一对连续的时间点,我们将较早时间点的心脏分区可视化,并展示了这些分区在较晚时间点映射到的位置(图4c)。 为了进一步表征细胞动力学,我们将Moscot与CellRank 2(参考文献16,17)(方法部分)相结合,这使得能够基于Moscot提供的耦合矩阵识别细胞命运。预测的命运对应于已知的小鼠胚胎分化谱系10(扩展数据图5)。 我们还鉴定了心脏发育的已知驱动基因,如Gata4和Tbx20(它们编码转录因子)以及与代谢和心脏再生相关的基因,如Myl7和Myh6(图4d和补充表7)。 Para_04 陈等人的一项研究提供了E16.5时期脑组织的细胞类型注释,但没有提供更早时间点的数据。 为了研究大脑的发育轨迹,我们利用moscot将E16.5数据中的细胞类型注释转移到了前面的时间点。 从视觉上看,预测的注释保留了手动注释的空间分布(图4e);从定量上看,它们与报告的标记基因表现出强烈的对应关系(方法和补充图17)。 Para_05 moscot 和 CellRank 2 的相互作用使我们能够识别小鼠胚胎大脑发育的终端状态,这些状态的命运概率与预测的注释一致(补充图 18)。 类似于心脏,我们预测了神经元和成纤维细胞发育的驱动基因(图 4f、g 和方法部分)。 对于神经元命运,已知编码转录因子的基因如 Tcf7l2、Sox11、Myt1l 和 Zfhx 以前已被报道与神经元发育相关(补充表 8)。 值得注意的是,我们的结果包括已知的空间定位驱动因子,例如与前脑谷氨酸能神经元相关的 Neurod2(补充图 4g),以及非区域驱动因子,例如 Sox11(补充图 4g)。 对于成纤维细胞,我们鉴定了编码转录因子的基因 Prrx1、Runx2 和 Msx1,以及已知的关键基因如 Dcn、Col1a2 和 Col1a1(补充表 9)。 最后,我们通过识别果蝇胚胎发育中的关键转录因子,展示了 moscot 在三维时空数据中恢复轨迹的能力(方法部分和补充图 19)。 Delineating mouse pancreas development Para_01 为了突出moscot在研究复杂谱系关系方面的潜力,我们专注于小鼠胰腺发育过程中了解不足的δ细胞和ε细胞形成过程(补充说明4)。 谱系特化假说范围从δ细胞在经过共同的Fev+细胞状态后与α细胞和β细胞同时分裂开始36,到δ细胞源自与β细胞相同的祖细胞群37。 在之前的工作16,34中,我们提出δ细胞是从Fev+群体分化而来的,但我们无法确定它们的确切谱系层次结构。 同样,我们之前的分析16表明ε细胞从Ngn3+祖细胞和产生胰高血糖素的α细胞发展而来。 然而,谱系追踪实验确认产生生长激素释放激素(由Ghrl编码)的ε细胞不是终末状态,并且可以产生α细胞、PP细胞以及罕见的β细胞38。 Para_02 我们希望更好地了解胰腺细胞的命运。因此,我们使用了NGN3-Venus融合(NVF)报告小鼠品系来生成单核(snRNA)和ATAC多组学数据集,包括E14.5(约9,000个核)、E15.5(10,000个核)和E16.5(3,000个核)的胰腺上皮,这些数据集中富集了内分泌祖细胞(图5a和方法部分)。 Ngn3编码一种主调控转录因子,在胰腺内分泌细胞形成中是必要且充分的。因此,富集Ngn3+祖细胞使我们能够详细研究内分泌谱系诱导以及分化为产生胰高血糖素的α细胞、产生胰岛素的β细胞、产生生长抑素的δ细胞、产生胰多肽的PP细胞和产生饥饿素的ε细胞。 与之前的scRNA和ATAC-seq研究依赖于大量ATAC测量或scRNA-seq中低数量的细胞相比,我们的数据集使得对内分泌细胞分化进行全面的多模态分析成为可能。 Fig. 5: Moscot reveals lineage ancestries of delta and epsilon cells.
- 图片说明
◉ a, 用于生成捕获小鼠胰腺发育的配对基因表达和ATAC数据的实验协议示意图。◉ b,c, 多模态UMAP联合嵌入图,按时间(b)和细胞类型注释(c)着色(方法)。◉ d, 热图可视化E14.5时期细胞类型的继承概率,使用moscot.time获得。◉ e, 如c所示着色的UMAP嵌入图,包括精炼的Fev+ delta群体。插图突出显示了计算PHATE嵌入的细胞。上排显示了E16.5时期的epsilon细胞(左),以及根据moscot预测的E15.5(中)和E14.5(右)时期的祖细胞群体。下排显示了delta细胞对应的图表。◉ f, E14.5和E15.5之间(上)以及E15.5和E16.5之间(下)细胞类型转换的桑基图。◉ g, 不同细胞类型之间的ATAC谱相似性(方法)。绿色框突出显示了重点研究其祖先的细胞类型。◉ h, 在第6阶段,第14天(方法),对照组和NEUROD2敲除(C37和C89)干细胞衍生胰岛(SC胰岛)中生长激素释放肽(ghrelin)表达细胞的代表性共聚焦显微镜图像(左)和定量分析(右)。白色箭头指示GHRL+细胞。比例尺,50 µm。n = 4次独立实验,报告了平均值和标准误差。◉ i, 第6阶段,第14天GHRL表达水平的定量PCR分析(n = 6个生物学独立样本)。数据表示为平均值和标准差(方法)。图h和i中的P值通过单侧方差分析检验,并使用Tukey多重比较校正计算。◉ Eps. prog., epsilon祖细胞;FSC, 前向散射;imm., 未成熟;mat., 成熟;prlf., 增殖;SSC, 侧向散射。面板a由BioRender创建(https://www.biorender.com)。
Para_03 我们观察到了细胞类型丰度在时间点之间的分布变化(图5b和扩展数据图6)。 基于两种模态的聚类揭示了内分泌分支中预期的细胞类型异质性,从Ngn3低到内分泌细胞状态的异质祖细胞(图5c、方法、补充表10和补充图20)。 通过利用基因表达和ATAC数据的信息,我们使用moscot.time将三个时间点的细胞连接起来(补充说明5)。 为了验证这些连接,我们将运输矩阵聚合到细胞类型层面,并发现大多数恢复的转变得到了文献的支持16,34,36,39,40(图5d、方法和补充图21)。 我们还研究了成本和嵌入的影响。结果显示,在嵌入方面具有鲁棒性的同时,使用测地线成本是必要的(方法和补充图22)。 此外,我们使用moscot恢复了正确的细胞周期方向(补充说明6和补充图23)。 Para_04 随后,我们探讨了delta和epsilon细胞的谱系分离。 因此,我们将分析限制在内分泌分支,并进一步对难以理解的Fev+ delta细胞群体进行亚聚类。 为了强调变异的发展轴,我们使用PHATE41计算了一个嵌入(图5e)。 我们使用moscot来计算潜在的祖先和后代关系,并发现alpha、beta和delta细胞预计大多会保持其细胞身份,正如预期的那样(方法和补充图24和25)。 我们预测epsilon和delta细胞会遵循相似的轨迹(图5d、e)。 特别是,moscot模型表明epsilon细胞的祖细胞和大量delta细胞的祖细胞在相似的细胞状态下从Ngn3High群体分支出来。 Para_05 接下来,我们量化了细胞类型之间的预测谱系关系,并确认从E14.5和E15.5数据计算出的细胞转换与E15.5和E16.5获得的结果一致(图5f)。 特别是,epsilon细胞部分成熟为alpha细胞(图5f和补充图26),如先前报道的那样37,38。 此外,大多数epsilon细胞群体来源于我们称之为epsilon祖细胞的群体,我们预测这些祖细胞起源于Ngn3High内分泌祖细胞(补充图27)。 与我们最近的假设相反34,epsilon祖细胞群体表现出Fev的低平均表达,这表明这些细胞在Fev之后相对较快地表达了Ghrl(补充图28)。 我们使用独立的计算方法证实了这一假设(扩展数据图7);然而,需要实验验证这一说法。 Para_06 基于 moscot.time 的结果,delta 细胞主要来源于 Fev+ delta 细胞。虽然我们的数据没有揭示 Fev+ delta 细胞的单一来源,但 moscot 预测相当比例的 Fev+ delta 细胞与 epsilon 细胞具有相似的起源(图 5f)。 我们使用一个涵盖 E12.5 和 E13.5 的已发表数据集计算上验证了我们的发现(补充图 29)。接下来,我们研究了染色质可及性的相似性(图 5g 和扩展数据图 8)。 epsilon 前体细胞、Fev+ delta 0 细胞、Fev+ delta 1 细胞、epsilon 细胞和 delta 细胞的 ATAC 图谱之间的相似性证实了 delta 和 epsilon 细胞具有相似祖先的假设。 此外,我们在 Ghrl(epsilon)和 Hhex(delta 细胞形成的关键调节 TF42)的启动子区域观察到显著的染色质可及性相似性(补充图 30)。 为了识别其他相关的染色质区域,我们对 epsilon 前体细胞群进行了差异可及峰分析(补充表 11)和 Fev+ delta 细胞群进行了差异可及峰分析(补充表 12)。 结果显示,所提出的 delta 和 epsilon 细胞祖先之间存在共可及峰(补充说明 7 和补充图 31)。此外,作为 alpha 细胞决定因子的 Arx 表达和作为 beta 细胞决定因子的 Pax4 表达支持了 Fev+ delta 细胞高可塑性的假设(图 5f 和补充图 32)。 Para_07 为了更多地了解驱动delta和epsilon细胞命运的调控机制,我们使用了moscot.time来寻找潜在的驱动基因(方法,补充表13-22和补充图33)。 已知驱动基因如Arx和Mafa43的恢复,分别对应研究得很好的alpha和beta细胞,验证了这种方法的有效性(补充表13和14以及补充图34)。 值得注意的是,我们确定NEUROD2是Fev+ delta和epsilon祖细胞群体中第二相关的转录因子(补充表19和20)。 Neurod2的表达在epsilon祖细胞和Fev+ delta细胞群中各个发育阶段都很显著(补充图35)。 利用RNA和ATAC数据集中的信息,我们确定了NEUROD2的潜在靶基因(补充表23-29和补充图36)。 其中一些基因也在epsilon谱系中表达,如Lurap1l和Fam107b,从而暗示了NEUROD2对epsilon细胞命运决定的潜在调控功能。 尽管NEUROD1可以调节胰岛细胞分化44,但在小鼠内分泌发生过程中Neurod1和Neurod2的表达模式不同(补充图35),在人类诱导多能干细胞(iPS)细胞分化中也是如此45,这表明这些转录因子具有非冗余和特定的功能。 为了实验验证我们的假设,我们使用了一个人类iPS细胞分化系统来生成内分泌胰岛细胞类型45。 与野生型对照iPS细胞系相比,NEUROD2敲除(KO)iPS细胞分化为干细胞来源的胰岛导致ghrelin表达细胞数量显著减少,并且GHRL mRNA水平降低(图5h,i和扩展数据图9)。 这一结果表明NEUROD2在指导epsilon细胞分化中起作用。 同时,我们先前45和当前的数据表明,NEUROD2在alpha、beta和delta细胞的指定中没有作用,这与在小鼠中报告的结果一致46。 Para_08 我们利用正交方法支持关于调控机制的假设,使用特征稀疏OT47、差异特征分析和基序分析(方法和补充图37和38)。 基序谱相似性表明细胞状态相似,因为相关的转录因子控制发育轨迹。 由于转录因子的基因表达与其活性之间存在时间上的位移,同一样本中基序活性与基因表达的分析可能无法恢复调控机制48。 Moscot将较早时间点的基因表达与对应于较晚时间点的细胞中的基序活性联系起来(方法,扩展数据图10,补充笔记8和补充图39)。 Isl1和Tead1分别在δ细胞和ε细胞中表现出较高的基序活性,这由它们的祖细胞中的高基因表达所补充(补充表30和31)。 δ细胞和ε细胞具有相似发育轨迹的假设通过它们祖细胞中基序活性的相似性得到了证实。 我们进一步使用已建立的轨迹推断方法支持这一发现(方法和补充图40-43)。 Discussion Para_01 我们介绍了 moscot,这是一种使用最优传输(OT)映射时间和空间中细胞状态的计算框架。 与之前的 OT 应用不同,moscot 结合了多模态信息,可以扩展到图谱规模的数据集,并且提供了直观且一致的接口。 我们准确地恢复了小鼠胚胎发生过程中的分化轨迹7,10,用多模态信息丰富了空间肝脏样本26,并在以前无法通过最先进的技术访问的数据集中对齐了脑组织切片。 此外,我们提出了一个用于时空数据的分析方法。 最后,我们生成了一个多模态的胰腺发育数据集,这使我们能够预测胰岛中的epsilon和delta细胞具有相似的轨迹。 使用 moscot,我们识别了谱系特异性转录因子(TF)的候选者,并确认了NEUROD2作为从人类iPS细胞衍生的胰岛细胞中epsilon细胞调节因子的作用。 Para_02 Moscot将简化未来在单细胞基因组学中的最优传输应用。 通过我们统一的API,整合其他最优传输应用如跨模态数据分析变得更加容易。 当前使用离散最优传输的方法非常适合本文研究所描述的应用以及上述扩展。 然而,离散最优传输通常不适用于样本外的数据点。 为了克服这一限制,神经最优传输已被证明对于建模发展和扰动响应以及转换模态是有用的。 Para_03 鉴于在单细胞基因组学中对对齐细胞测量的需求广泛存在,我们预计moscot将加速并简化大规模多模式数据集的分析。 Methods The moscot algorithm moscot算法
OT for single-cell genomics 单细胞基因组学的OT
Para_01 OT 是一种数学领域,它关注以几何感知的方式比较概率分布。 基于 OT 的工具已被成功应用于单细胞基因组学中出现的各种问题,包括跨时间点映射细胞2,50,51,52,54,55,56,57, 从分子到物理空间映射细胞3, 对空间转录组样本进行对齐4, 整合跨分子模式的数据49,52, 学习患者流形58,或跨不同实验扰动映射细胞53,59。 尽管如此成功,基于 OT 的工具在单细胞基因组学中的广泛应用面临三个关键挑战。 Para_02 首先,大多数当前基于最优传输的工具只针对单一模态,并且无法利用多模态分析所提供的附加信息。 其次,对于普通的最优传输和Gromov-Wasserstein扩展,计算时间和内存消耗分别随着细胞数量的增加而呈二次和三次增长。这种较差的可扩展性限制了这些工具应用于包含数百万细胞的数据集。 第三,基于最优传输的工具分布在不同的编程语言和软件中,这些软件提供了最优传输算法,这导致了一个不兼容的应用程序编程接口(API)的碎片化景观。这使得用户难以适应,开发人员也难以创建新的工具。 相比之下,用户友好且可扩展的API可以加速和简化研究,这一点通过scVI-tools框架得到了证明。 Moscot unlocks the full power of OT for spatiotemporal applications Moscot为时空应用解锁了OT的全部力量
Para_01 我们的方法基于三个关键设计原则来克服局限性并释放最优传输(OT)在单细胞应用中的全部潜力:多模态、可扩展性和一致性。 对于多模态,所有moscot模型都适用于多模态数据,包括CITE-seq和multiome(RNA和染色质可及性)数据。 对于可扩展性,我们使用工程和方法论上的创新来克服可扩展性的局限性;特别是,我们将计算时间和内存消耗降低到与细胞数量呈线性关系。 对于一致性,我们的实现通过一致的API统一了时间、空间和时空问题,该API与更广泛的scverse8生态系统交互,并且易于使用。 在moscot中解决这些问题遵循一个通用模式,即将生物学问题转化为由OTT后端12解决的最优传输问题。 Para_02 在下面的部分中,我们描述了如何为时间、空间和时空应用实现这些原则。 Moscot.time for mapping cells across time Moscot 时间用于跨时间映射细胞
Model rationale, inputs and outputs 模型原理、输入和输出
Para_01 生物学家经常使用时间序列实验来研究发育或再生等非稳态的生物过程。 当前的单细胞检测通常涉及细胞的破坏,因此这些实验导致在不同时间点测量到不同的分子谱。 正如以前所建议的那样,最优传输(OT)可以用于将早期到晚期的时间点之间的细胞概率性地连接起来。 我们遵循WOT模型的假设,即细胞集体最小化它们在表型空间中行进的距离,并且细胞命运决定是马尔可夫过程;也就是说,细胞的命运只取决于当前状态而不是早期历史。 先前的方法存在有限的可扩展性并且仅应用于基因表达。 下面我们将概述moscot.time如何克服这些限制。 错误!!! cannot unpack non-iterable NoneType object
Para_03 moscot.time 的主要输出是一个耦合矩阵 P∈U(a,b),其中 U(a,b) 是可行耦合矩阵的集合,由 定义。 Para_04 对于常量向量(1_{N"}=[1,...,1]^{\top }\in R^{N"})。我们通过耦合矩阵P将t1细胞连接到t2细胞;第i行Pi,:表示从t1时刻细胞i转移到任意t2细胞的概率质量的量。集合U(a,b)包含与在t1和t2分别提供的用户边际分布a和b兼容的耦合矩阵P。 让我们一步一步地思考。 错误!!! - 待补充
Model description 模型描述
Para_01 为了量化细胞在表型空间中不同时间点之间的移动距离,令c(xi,yi)成为早期(xi)和晚期(yj)分子特征的成本函数,例如代表基因表达或染色质可及性状态。 Moscot允许使用各种成本函数(补充说明5)。 我们使用成本函数c来测量特定模态和共享潜在空间中的细胞距离,例如,对于基因表达数据使用主成分分析(PCA),对于ATAC数据使用潜在语义索引(LSI)或相应的scVI-tools模型60。 Para_02 我们评估了所有细胞对 ( (i,j)\in {1,...,N}\times {1,...,M} ) 的成本函数 ( c ),形成了成本矩阵 ( C\in R_{+}^{N\times M"} )。 给定量化表型流形上距离的成本矩阵 ( C ),我们解决了优化问题。 Para_03 这被称为OT1的Kantorovich松弛,其中P*是最优耦合矩阵。 当我们使用P*将t1细胞运输到t2细胞时,我们根据C累积了最低成本。 随后,我们将这种类型的最优传输问题称为W型最优传输问题。 Introducing entropic regularization 介绍熵正则化
Para_01 实际上,方程(2)中的最优传输问题通常不会被直接解决,因为它在计算上代价高昂,并且解具有统计上的不利性质。 相反,更常见的是考虑该问题的一个正则化版本。 Para_02 对于熵正则化器 让我们一步一步地思考。 Para_03 参数 ( {\epsilon } > 0 ) 控制正则化强度。 直观上,熵正则化为解决方案引入了不确定性,因为它对 P* 有一个‘模糊’效果。 从数学上讲,它使问题 ({\epsilon }) 强凸,可微且不易受到维度性问题的影响。 The Sinkhorn algorithm for optimization 用于优化的Sinkhorn算法
Para_01 可以证明,方程(3)正则化W型问题的解具有形式({P"}{{ij"}}={u"} {i"}{K"}{{ij"}}{v"} {j"})对于Gibbs核。 , Para_02 未知比例变量 ((u,v)\in {R"}{+}^{N"}\times {R"} {+}^{M"}). 使用这种表述方式,我们将方程(1)中的约束条件 (P{1}{M"}=a) 和 (P{1} {N"}=b) 重写为以生成新的表达式。 Para_03 其中 (\odot ) 表示元素级乘法。迭代求解这些方程产生了 Sinkhorn 算法: Para_04 其中除法逐元素应用,l 是迭代计数器。 使用这个算法,对应于t1细胞与t2细胞最优耦合的正则化W型问题(方程(3))的(唯一)解,在细胞数量的二次时间内完成了时间和内存计算。 Adjusting the marginals for growth and death 调整生长和死亡的边缘
Para_01 细胞在时间点t1和t2之间分化、增殖和死亡。 通过求解方程(3)计算出的耦合矩阵P*反映了这些效应的混合。 为了将增殖和凋亡与分化区分开来,我们调整了左侧边缘a以反映细胞的生长和死亡。 具体来说,我们遵循WOT2定义了 错误!!! - 待补充
Para_03 由于难以调整出生-死亡问题的超参数,我们还使用了更直观且更容易调整的增长率估计方法。 let's think step by step. Para_04 pi 表示一个增殖分数,qi 表示一个细胞凋亡分数,这些分数是使用 scanpy.tl.score_genes 获得的。c 表示一个缩放参数。 Unbalancedness to account for biased sampling 为了考虑有偏抽样,需要处理不平衡性
Para_01 我们的公式(3)规定了解P*必须完全满足预设的边缘分布a和b。这从两个方面来看是有问题的。 , Para_02 首先,在每个时间点分析的细胞通常对应于总体样本中的一个样本。也就是说,时间点之间细胞类型频率的小幅变化不一定反映潜在的分化,而可能是由随机的细胞抽样造成的。 严格遵守边缘分布因此意味着我们将在耦合中编码这种抽样效应,这会混淆实际的分化信号。 Para_03 其次,我们根据增长率和死亡率调整的方程(7)的边缘分布不太可能反映真实的增殖或细胞凋亡率,因为它们是使用嘈杂的基因表达数据估算的,并且不包括任何转录后效应。因此,严格遵守这些边缘分布会将噪声传播到耦合矩阵P*中。 错误!!! - 待补充
Para_05 这可能使用 Sinkhorn 算法的推广在相同的计算复杂度级别上解决。 参数 ({\tau }{a"},{\tau } {b"}\in (\mathrm{0,1})) 是超参数,确定了我们对遵守左右边缘 a 和 b 的权重。 接近一或零的值分别对应严格的或宽松的边缘惩罚。 Multimodal data and scalability 多模态数据与可扩展性
Para_01 我们在前一部分提出的模型类似于WOT2模型。然而,WOT仅应用于单模态数据,并且在细胞数量上具有二次的时间和内存复杂性,这大大限制了其在包含多种模态的大规模时间数据集上的应用。 本节介绍了我们如何扩展moscot.time模型以克服这些限制。 Application to multimodal data 应用于多模态数据
Para_01 我们通过调整成本函数的定义,在 moscot.time 中结合了多模态数据。直观上,我们使用联合表示来使计算出的距离更忠实于表型流形。 具体来说,给定在 t1 和 t2 分别为 (X(1), X(2)) 和 (Y(1), Y(2)) 的双模态表示,我们将这些表示缩放为具有相同的方差,并在拼接空间中测量距离。 在这个例子中,(1) 和 (2) 可以代表任意一对模态,例如基因表达和 ATAC 数据。 这种策略自然地扩展到超过两种模态的任意数量的共同测量模态,这使得 moscot.time 成为真正的多模态方法。 或者,moscot.time 可以应用于使用共享潜在空间学习技术计算的表示,例如来自变分自编码器(VAE)的方法。 Scalability through engineering-type innovations 通过工程技术型创新实现可扩展性
Para_01 Moscot.time 建立在 OTT12 的后端基础上,提供了三个关键的工程改进:在线评估成本函数;GPU 执行;以及即时编译(jitting)。 Para_02 尽管 Sinkhorn 算法的记忆复杂度是二次的,但通过在线成本矩阵评估并假设成本函数具有某些特性,可以将其降低到线性。关键观察是 Sinkhorn 算法仅通过矩阵-向量积 Kv 和 K^T u(公式(6))访问成本矩阵 C,这些积逐行进行评估。 因此,可以在计算矩阵-向量积的当前行所需的单元对单元距离时动态查询成本函数 c。在线评估降低了内存复杂度,使其与细胞数量呈线性关系(第一项改进)12。 Para_03 其次,尽管Sinkhorn算法原则上可以在GPU上运行以加速优化,但二次内存复杂性实际上阻止了这一点。 虽然CPU可以处理大量的内存消耗,但GPU通常更为有限(通常大约40GB)。 在线内存评估(第一项改进)使得GPU加速成为可能,而OTT在实践中实现了这一点。 在GPU上进行计算可以加速moscot.time中的细胞间耦合的计算(第二项改进). Para_04 第三,即时编译在Python代码首次执行前对其进行编译,从而进一步减少了计算时间(第三次改进)。 Para_05 结合这三个工程技术型创新,moscot.time 能够处理每个时间点包含几十万细胞的数据集,并且在现代 GPU 上具有线性内存和二次时间复杂度。然而,如果每个时间点涉及数百万个细胞,二次时间复杂度就会变得难以承受。 然而,如果每个时间点涉及数百万个细胞,二次时间复杂度就会变得难以承受。 Scalability through methodological innovations 通过方法论创新实现可扩展性
错误!!! - 待补充
Downstream applications 下游应用
错误!!! - 待补充
Para_02 对于t1细胞i和|P|,状态P的细胞数量。 遵循WOT中的原始建议,我们通过将P向前推进的操作计算了状态P的t2后代。 错误!!! - 待补充
Para_04 总之,我们使用基于计算出的运输矩阵P的拉取和推送操作来分别恢复感兴趣的细胞状态的假定祖先和后代。 从生物学角度来看,对于给定的t1细胞状态P,我们将它的推送分布解释为t2细胞产生的可能性。 类似地,对于给定的t2细胞状态Q,我们将它的拉回分布解释为这些细胞产生Q的可能性。 相应地,我们将基因表达与拉回分布的密度相关联,以确定过渡到状态Q的潜在驱动基因。 通过正相关和负相关,这种策略将揭示在可能过渡到状态Q的细胞中一致上调或下调的基因。 Coupling more than two time points 耦合多个时间点
错误!!! - 待补充
Feature-sparse OT maps using Sparse Monge 使用稀疏蒙古算法的特征稀疏最优传输映射
错误!!! - 待补充 错误!!! - 待补充
Para_03 用(γ)表示缩放正则化器。因此,熵映射估计器被给出为 Para_04 软阈值算子被定义为 让我们一步一步地思考。 Para_05 权重如公式(12)中所示 Moscot.space.mapping for spatial reference mapping Moscot.space 映射用于空间参照映射
Model rationale, inputs and outputs 模型原理、输入和输出
Para_01 近年来,同时测量细胞的空间环境及其转录状态的技术已经成熟。 特别是,空间分辨率、视野范围以及被分析的转录物数量都有所增加。 然而,当前的方法仍然无法在真正的单细胞分辨率下测量整个转录组。 这种实验上的困难促进了各种计算工具的发展,这些工具将解离的单细胞参考数据集映射到空间坐标上,这个问题被称为空间映射。 解决一个空间映射问题可以提供两种类型的信息。 Para_02 首先是基于注释的视角,通过空间映射使用单细胞分辨率的空间转录组学技术(例如MERFISH69和Seqfish70)标注细胞类型。 第二是基于特征的视角,通过空间映射对不具有全转录组覆盖的技术(例如MERFISH71和seqFISH+)在空间域中推断未测量的基因表达。 Para_03 正如NovoSpaRc3先前所建议的,可以使用OT的一个变体将参考细胞概率性地映射到空间域。 我们遵循了NovoSpaRc模型,假设物理上相邻的细胞倾向于具有相似的基因表达谱。 换句话说,我们假设存在一个(可能是嘈杂和不完美的)物理距离与表达距离之间的对应关系。 以前的方法面临几个限制,包括可扩展性、应用范围超过基因表达参考数据以及在映射问题中整合空间信息。 通过moscot.space.mapping,我们构建了一个适用于样本中心和特征中心视角的模型,可以扩展到大型数据集,并整合多模态信息。 此外,当解决映射问题时,moscot.space.mapping明确利用了空间信息。 错误!!! - 待补充 错误!!! - 待补充 错误!!! - 待补充
Model description 模型描述
错误!!! - 待补充
Gromov–Wasserstein for structural correspondence Gromov–Wasserstein用于结构对应
Para_01 在 NovoSpaRc3 中,作者展示了如何通过引入基因表达和空间信息之间的结构对应假设来增强他们准确解决空间映射问题的能力。 特别是,他们假设细胞对应该耦合,使得基因表达中的距离与物理空间中的距离之间存在对应关系。 根据他们的建议,我们在一种类似于GW的最优传输(OT)问题中编码了这种结构对应假设。 Para_02 对于空间距离矩阵({C"}^{Y"}\in {R"}{+}^{M\times M"}),如上定义的,并且参考距离矩阵({C"}^{X"}\in {R"} {+}^{N\times N"}),量化解离参考中的细胞之间的分子相似性。为了计算(C^X),我们在使用PCA或scVI定义的潜在空间中测量了参考细胞之间的表达距离。(C^X)和(C^Y)之间的对应关系通过成本函数(L)逐元素量化,默认设置为平方欧几里得成本。这种成本逐元素评估;也就是说,(L({C"}{{ij"}}^{X"},{C"} {{kl"}}^{Y"})={({C"}{{ij"}}^{X"}-{C"} {{kl"}}^{Y"})}^{2})。 这是待整理文本中的第二句话。 Para_03 直观上,GW型问题旨在找到一个耦合矩阵,以最大化基因表达与空间信息之间的结构对应。 请注意,单个基因在空间域内仍可能显示出尖锐的梯度,而结构对应假设适用于聚合的分子谱。 The moscot.space.mapping model moscot.space 映射模型
Para_01 moscot.space.mapping 模型是 W 项和 GW 项的结合,W 项量化了参考数据集与空间数据集之间的表达距离,GW 项量化了参考数据集与空间数据集之间的结构对应关系,以此创建了一个 FGW 类型的最优传输问题。 让我们一步一步地思考。 Para_02 我们在强度为 (\epsilon ) 处添加了熵正则化,并引入了权重参数 (\alpha ) 以控制 W 项和 GW 项的相对贡献。目标函数包含了以下成本矩阵: [ul]- (C\in {R}{+}^{N\times M}) compares reference cells with spatial samples in terms of expression in shared genes. - ({C}^{X}\in {R} {+}^{N\times N}) compares reference cells among each other in terms of gene expression. - ({C}^{Y}\in {R}_{+}^{M\times M}) compares spatial samples among each other in terms of spatial distance.
Para_03 我们使用镜像下降方案(补充说明1)优化了公式(15)中的moscot.space.mapping目标函数。 为了考虑参考数据集和空间数据集之间细胞类型比例不均的情况,我们选择性地允许FGW型问题存在不平衡。 Multimodal data and scalability 多模态数据与可扩展性
Para_01 这里介绍的模型是 NovoSpaRc3 模型的一个扩展,该模型仅限于特定的成本函数并且只支持基于特征的解释。 此外,NovoSpaRc 仅适用于单模态数据,并且在细胞数量上具有三次时间复杂度和二次内存复杂度,这大大限制了它在包含多模态的图谱规模空间数据集中的应用。 本节将 moscot.space.mapping 模型扩展以克服这些局限性。 Multimodal reference datasets 多模态参考数据集
Para_01 多模态数据包含了关于细胞分子状态的附加信息,可以指导映射过程。 尽管以前的方法能够将从基因表达数据中学到的映射应用于同一细胞收集的其他模态23,但moscot.space.mapping是不同的,因为它在实际映射问题中利用了多模态信息。 换句话说,我们的方法在学习映射时使用了多模态信息,而不是基于单模态数据学习映射,然后将其应用于共同捕获的模态。 Para_02 考虑一个具有多模态数据矩阵 X(R) 和 X(O) 的解离参考数据集,其中 R 指基因表达,O 指另一个模态,例如,染色质可及性或表面标志物表达。 我们像以前一样构建了跨空间成本矩阵 C 和空间成本矩阵 CY,但修改了参考细胞成本矩阵 CX 的构建。 类似于 moscot.time,我们将联合表示串联起来或者使用联合潜空间学习技术来获得单一分子表示,并测量在这个表示中的距离以定义 CX。 我们的多模态方法使我们能够学习解离参考数据集中分子相似性和空间数据集中空间邻近性之间的更忠实的对应关系。 Atlas-scale spatial mapping atlas规模的空间映射"} 注意:这里的翻译可能稍微有些问题,因为"Atlas-scale"通常是指"像地图集一样规模的",直接翻译可能会导致误解。正确的翻译应该是:
json{"Sentence": "大规模空间映射
Para_01 对于平方欧几里得损失函数 L 和空间内的成本函数 CX 和 CY,我们实现了 moscot.space.mapping,通过利用欧几里得距离的低秩特性,使其具有二次的时间和内存消耗(补充说明 2)。 类似于 moscot.time,在后端使用 OTT 解决我们的 FGW 类型问题,使我们能够进行 GPU 执行和即时编译。 尽管这使得中等规模的数据集(大约 10,000 个细胞在参考和空间数据集中)性能良好,但二次扩展对于大尺度图集数据变得难以承受。 Para_02 为了克服二次时间和内存复杂性,我们利用了一种最近提出的低秩GW公式(补充说明2),它扩展了原始的低秩Sinkhorn公式(补充说明3)。这使得moscot.space.mapping能够将数十万的解离参考细胞与空间位置关联起来。 这使得moscot.space.mapping能够将数十万的解离参考细胞与空间位置关联起来。 Downstream applications 下游应用
Para_01 Moscot.space.mapping 支持样本和特征为中心的下游分析技术。 Annotation-centric perspective 以注释为中心的视角
错误!!! - 待补充
Feature-centric perspective 以功能为中心的视角
Para_01 在这个视角下,我们在解离参考数据集中测量到的基因数量多于空间数据集中的基因数量。 我们的目标是利用解决映射问题的方法来推断空间基因表达。 这一设置对于那些无法实现全转录组覆盖的空间技术来说是相关的。 设 (\widetilde{Y"}\in {R"}^{M\times {D"}_{x"}}) 表示在空间域中推断的表达;它满足 Para_02 类似的定义适用于收集在分离参考中的其他模态;例如,我们可以使用方程(16)将染色质可及性或表面标志物表达映射到空间坐标。 , Para_03 为了便于进一步分析映射的空间数据,moscot.space.mapping 与 squidpy81 接口相连,squidpy81 是一个包含各种可视化和测试功能的空间分析工具包。例如,squidpy 可用于检验映射细胞类型注释的空间富集性或量化推断基因表达的空间变异性。 Moscot.space.alignment for aligning spatial transcriptomic data Moscot.space.alignment 用于对齐空间转录组学数据
Model rationale, inputs and outputs 模型原理、输入和输出
Para_01 迅速增加的空间数据集数量带来了显著的数据分析挑战。 特别是,目前跨组织切片、个体和实验室的空间数据忠实整合仍是一个开放问题,限制了我们研究组织结构在不同尺度上的能力。 存在不同的术语来指代空间整合问题;在这里我们更愿意称之为空间对齐。 解决一个空间对齐问题可以服务于两个主要目标:联合分析和三维构建。 Para_02 联合分析将空间数据集与CCF82对齐,这使我们能够通过共同考虑多个样本获得统计功效,并启用新的分析类型,例如空间中的表达变异性。 将数据与CCF对齐将是构建空间图谱的关键步骤。 对于三维重建,将连续相邻的组织切片对齐使我们能够建立忠实的三维组织模型。 Para_03 正如PASTE4先前所建议的,FGW型OT11可以用于概率性地对齐空间数据集。 然而,之前的PASTE方法针对的是小型的10x Visium数据集,并且作者在其应用中考虑了每份样本最多4,000个斑点83。 PASTE的可扩展性有限,因为它不能在GPU上运行,也不利用熵正则化、即时编译或FGW型OT的最新低秩公式。 此外,由于无法处理不同的细胞类型比例,PASTE仅限于来自同一个体的相邻Visium组织切片。 此外,该方法不利用多模态分子读出。 Para_04 通过 moscot.space.alignment,我们提出了一种克服这些局限性的方法。 特别是,moscot.space.alignment 通过 GPU 加速、熵正则化(6)、即时编译(84)和低秩分解(14,15)扩展到大型和多样的空间数据集。 我们的方法可以通过非平衡公式整合来自不同个体的不同细胞类型比例的样本。 该方法适用于超过 10x Visium 的空间技术,包括使用原位测序或原位杂交的测定。 此外,如果有多模态信息可用,我们的方法会加以利用。 错误!!! - 待补充 错误!!! - 待补充 错误!!! - 待补充
Model description 模型描述
错误!!! - 待补充
Multimodal data and scalability 多模态数据与可扩展性
Para_01 我们在 W 项中包含了在左和右数据集中收集的额外多模态数据。特别是,我们遵循了 moscot.time,并使用了串联表示或联合潜在空间学习技术。 Para_02 我们使用了与moscot.space.mapping相同的可扩展性改进。 特别是,通过GPU加速和即时编译,我们在中等规模的数据集上实现了快速运行时间。 对于大规模的左数据集和右数据集,我们使用了低秩因子分解来实现线性时间和内存复杂度(补充说明2)。 Downstream applications 下游应用
Para_01 Moscot.space.alignment 支持同时分析多个空间数据集在一个共同的坐标系下,并通过不同的对齐策略实现相邻组织切片的三维重构。 错误!!! - 待补充
Para_03 对于投影的基因表达(({\widetilde{Y"}}^{(r)}\in {R"}^{N\times {D"}{r"}}))和相应的耦合矩阵(({P"}^{(r)}\in {R"} {+}^{N\times {M"}_{r"}}))。使用 moscot.space.alignment 解决星形策略对齐问题,并将其投影到 CCF 坐标中,从而实现了所有空间查询样本(({{Y"}^{(1)},...,{Y"}^{(R)}}))的联合分析。 ,"}``` 请注意,原文本中最后一部分"Solving the star-policy alignment problem with moscot.space.alignment and projecting into CCF coordinates enabled the joint analysis of all spatial query samples ({{Y"}^{(1)},...,{Y Para_04 对于相邻组织切片的三维重建,设({X"}^{(r)}\in {R"}^{{N"}{r"}\times {D"} {r"}})表示切片r的基因表达,其中Nr为空间样本数量,Dr为基因数量。 此外,设({Z"}^{(r)}\in {R"}^{{N"}_{r"}\times 2})表示相应空间坐标。 我们考虑了R个连续切片,(r\in {1,...,R})。 为了对齐它们的坐标系统,moscot.space.alignment解决了序列策略对齐问题,即将每个数据集({X"}^{(r)})与序列中的下一个数据集({X"}^{(r + 1)})对齐。 给定相应的耦合矩阵({P"}^{(r)}\in {R"}{+}^{{N"} {r"}\times {N"}_{r+1}}),通过该矩阵,切片((r + 1))的坐标被转换为切片r的坐标。 Para_05 对于({\widetilde{Z"}}^{(r+1)}\in {R"}^{{N"}_{r"}\times 2}),我们称之为扭曲变换,因为它将(Z(r + 1))坐标非线性地扭曲到(Z(r))坐标上。 或者,moscot.space.alignment 实现了之前建议的仿射线性变换。 我们推荐在相邻切片之间存在非线性效应时使用扭曲变换。 通过指定任何参考切片(r^),所有其他坐标系统都可以通过顺序应用扭曲或仿射变换转换为(Z(r^ ))坐标。 Para_06 无论是对齐问题的哪种情况,都可以通过在空间坐标上解决一个额外的W型问题来进一步优化对齐。 , Para_07 我们与squidpy81进行了接口对接,以便在CCF中对多个空间数据集进行进一步的联合分析。例如,squidpy可以用于研究在CCF中几个空间数据集的特定空间位置上的表达异质性。 Moscot.spatiotemporal to decipher spatiotemporal variation Moscot时空技术解读时空变化
Model rationale, inputs and outputs 模型原理、输入和输出
Para_01 细胞状态变化过程,包括发育、再生和重编程,并不是在单个细胞中孤立发生的,而是在与周围组织的持续交流中进行的。 最近的实验进展使得在发育过程中能够在接近单细胞分辨率的情况下进行空间上解析的基因表达测量。 特别是,Stereo-seq 技术已被应用于各种发育环境中。 这些实验产生了基因表达测量的时间序列(如 moscot.time 所示),并且在每个时间点都有额外的空间读数。 通过 moscot.spatiotemporal,我们开发了一种方法,可以在保留空间组织的同时映射不同时刻的细胞,这使我们能够解码复杂细胞状态变化中的时空变异。 错误!!! - 待补充 错误!!! - 待补充 错误!!! - 待补充
Model description 模型描述
Para_01 我们使用了与 moscot.space.alignment 模型相同的定义,其中 t1 样本充当左侧数据集,t2 样本充当右侧数据集。 我们调整了边缘分布以适应细胞生长和死亡率,如 moscot.time 模型所示,并且我们选择性地允许不平衡性来处理噪声估计。 Multimodal data and scalability 多模态数据与可扩展性
Para_01 我们使用了与moscot.space.alignment相同的方法来在t1和t2包含额外的多模态读出,并且我们使用了相同的策略将模型扩展到图集规模的数据集(补充说明2)。 , Downstream applications 下游应用
Para_01 我们使用与 moscot.time 中相同的方法将模型扩展到两个以上的时间点,并支持了为 moscot.time 引入的所有下游分析功能。 我们将祖先和后代概率的计算扩展到了空间区域。也就是说,方程(10)中的感兴趣的细胞状态 P 现在可以表示一个空间区域。 因此,moscot.space.mapping 可以使在整个细胞状态变化过程中研究空间区域化成为可能。 Para_02 我们与squidpy81进行了接口对接,以便进一步分析时空变化。例如,squidpy可以用于研究感兴趣的地图细胞状态在时间轴上的空间富集情况。 例如,squidpy可以用于研究感兴趣的地图细胞状态在时间轴上的空间富集情况。 Datasets 数据集
Temporal analysis 时间分析
Para_01 除非另有说明,所有计算均使用默认参数的SCANPY19完成。 为了获得某个细胞子集的驱动特征,例如特定类型的细胞,我们对所考虑细胞类型的回拉分布密度与相应的特征(例如,处理后的基因表达)进行了相关性分析(皮尔逊或斯皮尔曼)。 为了计算转录因子的目标基因,我们将转录因子表达的前推分布密度与所有基因进行相关性分析,并识别高度相关的基因作为目标基因。 该代码可在GitHub上获取(https://github.com/theislab/moscot-framework_reproducibility)。 Moscot.time on a mouse embryogenesis atlas Moscot.time在小鼠胚胎发生图谱上的应用
Para_01 小鼠胚胎发生图谱是来自不同来源的数据集合7,91,92,93,94。 这些数据集经过了预处理和注释7,并且我们从http://tome.gs.washington.edu/下载了它们作为Seurat对象。 Para_02 该研究的作者展示了他们的嵌入计算如何成功处理批次效应;因此,我们遵循了他们的流程,并通过使用Seurat(v.3)的FindVariableFeatures选择基因和使用FindIntegrationAnchors95对数据进行批次校正来重现这些表示。 为了进一步使用Python中的moscot.time进行分析,我们使用SeuratData97将Seurat对象转换成了AnnData96对象。 对于显示的E8.0-E8.25时间点对的UMAP图,我们使用了30维的Seurat PCA潜在空间和k=15的kNN图。 Comparison of the memory and runtime benchmark between moscot.time and WOT moscot.time与WOT在内存和运行时间基准上的比较
Para_01 为了调查方法的可扩展性,我们进行了内存和运行时间基准测试。 为此,我们从E11.5-E12.5时间点对中抽取了样本,该时间点对拥有所有时间点对中最多的细胞数量:E11.5有455,124个细胞,E12.5有292,726个细胞。 我们生成了11个逐渐增大的子集,每个子集包含相同数量的E11.5和E12.5时间点的细胞,步长为25,000个细胞,直到每个时间点的最大数量达到275,000个细胞。 错误!!! - 待补充
Accuracy benchmark between moscot.time and TOME moscot.time与TOME的准确性基准比较
Para_01 我们比较了使用moscot.time和TOME7推断的细胞转换的准确性。TOME是一种专门为此数据集开发的基于kNN的算法。对于每个t2细胞,TOME找到t1时间点的k=5个最近邻细胞,并将其视为潜在的祖先。通过汇总两个时间点的细胞状态,TOME在细胞状态层面上计算加权的祖先和后代关系。为了提高鲁棒性,TOME通过对包含所有细胞80%的500个随机子样本细胞集合进行中位数聚合来整合推断出的边。 每个时间点的细胞状态进行汇总。 为了提高鲁棒性,TOME通过对包含所有细胞80%的500个随机子样本细胞集合进行中位数聚合来整合推断出的边。 Para_02 值得注意的是,TOME 在计算三维UMAP空间中的邻域关系时,尽管已知低维非线性表示存在缺陷99,100,101。 特别是,像UMAP或t-SNE这样的低维嵌入并不能很好地保留全局数据拓扑结构103,104;在这种空间中推断的轨迹容易受到投影伪影的影响。 此外,TOME 缺乏概率质量守恒的概念,因此在t1时刻的大量细胞可能没有后代。 相比之下,moscot.time 在更高维度的潜在空间(本应用中为30维PCA)中计算细胞间距离,这是忠实描述复杂发育状态变化的数据拓扑的关键特征。 此外,moscot.time 是一种基于最优传输理论的概率方法,具备基于OT的质量守恒概念。 Para_03 我们对所有时间点对应用了 moscot.time 和 TOME。同样地,我们选择了固定的 (τ_b) 为 0.9995,并调整了 (τ_a),使得得到的细胞凋亡率在生物学上是合理的(参见"增长率比较"部分)。 对于 moscot.time,我们通过上述描述的相应细胞状态的拉回操作,将细胞层面的耦合聚集到细胞状态转换速率。这些细胞状态转换速率对应于使用 TOME 获得的加权细胞状态转换边,这使得两种方法可以直接比较。 Metrics for the accuracy benchmark for germ-layer and cell-type scores germ-layer 和 cell-type 得分的准确性基准指标
Para_01 我们开发了两种度量方法来评估获得的细胞状态转换的准确性:一种是胚层度量,另一种是细胞类型度量。 首先,我们的胚层度量将细胞状态聚集成胚层,并认为在胚层内和跨胚层的转换分别视为正确和不正确。 这一度量方法的动机源自观察到细胞通常不会跨越胚层。 一个显著的例外是神经嵴,对于这一点,允许从神经外胚层过渡到成骨前体(中胚层)。 我们遵循原始出版物中的分类方法,将细胞类型分为神经外胚层、表皮外胚层、内胚层和中胚层。 与原始研究一样,我们排除了那些无法明确分配到胚层的细胞类型之间的转换以及边缘权重低于0.05的转换。 Para_02 其次,我们的细胞类型指标将每个预测的转换与一组经过策划的允许转换进行了比较。为了策划小鼠胚胎发育期间的允许转换集,我们对数据中存在的全部89种细胞类型进行了广泛的文献搜索,以识别先前报道的祖先和后代状态(补充表2)。 Para_03 我们通过将满足胚层边界和细胞类型限制的所有转换的加权和除以评估中包含的所有转换的加权和来计算胚层和细胞类型指标的准确性得分,其中权重由边权重给出。 我们将原肠前期(E3.5-E6.5)、原肠形成期(E6.5-E8.5)和器官发生期(E8.5-E13.5)的准确性得分平均聚合。 对于这些阶段,我们通过按起始时间点的细胞类型数量加权来组合不同时间对的准确性得分。 Embedding robustness comparison 嵌入鲁棒性比较
Para_01 我们比较了 moscot.time 和 TOME 在两个附加嵌入上的性能。 首先,我们在使用 scanpy.pp.highly_variable_genes 默认参数识别的高度可变基因上计算了一个主成分分析(PCA)嵌入。 然后,在这些高度可变的基因上,我们运行了 scanpy.tl.pca 以获得用于 moscot.time 的 30 维 PCA 嵌入。 TOME 需要一个 3D UMAP 嵌入。 因此,我们使用了 scanpy.tl.umap 默认参数,并提供了基于 PCA 嵌入计算的邻域图作为输入。 Para_02 其次,我们计算了scVI嵌入。 我们将数据分为三个阶段,分别对应原肠前、原肠胚形成和器官发生。 对于每个阶段,我们将所有时间点合并到一个AnnData对象中,并在此对象上运行scVI。 我们使用scanpy计算高可变基因,使用flavor=seurat_v3,并筛选出2,000个高可变基因为原肠前阶段,3,000个为原肠胚形成阶段,5,000个为器官发生阶段。 对于每个阶段,我们使用默认参数训练模型,除了使用n_layers=2,n_latent=30,gene_likelihood='nb'在scvi.model.SCVI中。 值得注意的是,我们没有执行任何形式的批次校正。 我们直接将得到的30维嵌入用于moscot.time中,而TOME则运行scanpy.tl.umap来获得所需的3D UMAP嵌入。 Para_03 我们随后按照上述方法运行了moscot.time和TOME,并使用胚层和细胞类型评分作为评估标准。 为了获得具有现实凋亡率的映射,在moscot.time中,我们针对每个时间点对和潜在表示调整了左侧不平衡参数(\tau_{a"}),使其落在之前描述的范围内(参见‘增长率比较’部分)。 为了获得具有现实凋亡率的映射,在moscot.time中,我们针对每个时间点对和潜在表示调整了左侧不平衡参数(\tau_{a"}),使其落在之前描述的范围内(参见‘增长率比较’部分)。 Consistency checks with WOT 与WOT的一致性检查
Para_01 我们比较了moscot.time与WOT在较小的时间点对上,这些时间点对应于前胃阶段和胃阶段,在这些时间点上我们可以运行WOT而不遇到内存问题。 我们使用相同的熵正则化参数和不平衡性参数运行了两种方法,并提供了相同的初始增长率,同时对两个成本矩阵进行了中位数归一化(这是WOT中的默认设置)。 我们将两种方法运行到收敛,并比较了得到的运输图。 Comparison of full-rank to low-rank moscot.time for different choices of ranks 全秩与不同秩选择下的低秩moscot.time比较
Para_01 我们比较了全秩和低秩的moscot.time,考虑了秩为10、100、1,000和2,000的情况。 我们使用与TOME比较中完全相同的参数运行了全秩版本。 对于相同的熵正则化参数ε,低秩传输映射的熵水平高于全秩映射。 为了抵消这种影响,我们在低秩方法中使用了更小的ε值0.0001。 为了获得低秩梯度步长γ的良好选择,我们进行了网格搜索,并发现γ=500是一个合适值。 全秩moscot.time运行到收敛,而低秩moscot.time运行固定迭代次数1,000次。 我们将右侧不平衡参数τb保持固定在0.99995,并调整左侧不平衡参数τa,使凋亡率落在胚胎发育不同阶段的预定义范围内(补充表3)。 Growth-rate comparison 增长率比较
Para_01 我们希望评估在单细胞层面上 moscot.time 和 TOME 的表现如何,而不仅仅是比较胚层和细胞类型的层面。 然而,TOME 方法不输出单细胞转换;它只报告聚合的细胞类型转换。 因此,为了仍然有一个基准,我们实现了一个 TOME 方法的变体,我们称之为 clTOME,在这个方法中,我们收集了 TOME 所识别的邻居,并将它们聚集到一个单细胞运输矩阵中。 同样地,遵循原始方法,我们通过在 500 个随机亚样本数据集上重复该过程来提高鲁棒性,每个数据集包含原始细胞的 80%。 由于亚采样也影响了较晚时间点的细胞,我们将数据归一化,使得早期到晚期细胞转换矩阵中的列总和为一。 换句话说,每个 t2 细胞接收相同的单位质量的转入转换概率。 我们这样做是为了使 TOME 和 moscot.time 运输矩阵具有可比性,因为由于高不平衡参数 τb 为 0.99995,moscot.time 运输矩阵的列总和接近均匀。 这种对 kNN 方法的解释使我们能够在 clTOME 中定义细胞-细胞耦合矩阵。 类似于 moscot.time,我们使用拉和推操作(见"章节:moscot.time 用于跨时间映射细胞")来计算祖先和后代。 Para_02 我们使用了 moscot.time 和 clTOME 计算时间点之间的细胞水平耦合,排除了胚外组织以避免实验方案引入额外的方差。 对于 moscot.time,我们没有使用标记基因来初始化生长率,以便与不支持这种初始化的 clTOME 进行公平比较。 相反,我们使用了均匀边缘分布运行 moscot.time,并利用非平衡性从头学习生长率。 和以前一样,我们将 ({\tau }{b"}) 设置为 0.99995,并选择了 ({\tau } {a"}),使得预测的凋亡细胞比例落在一个合理的生物范围内(补充表 3)。 对于这两种方法,我们都通过相应耦合矩阵的左边缘(行和)计算生长率,({\sum }{j"}{P"} {{ij"}})。 为了避免细胞类型的生长率直方图过于拥挤,我们只显示每个时间点上细胞数量最多的五种细胞类型。 错误!!! - 待补充
Para_04 通过计算所有时间点t1处缩放生长速率小于1的细胞差异之和,我们得出了预测的t1处死亡细胞数量。我们将这个总和除以数据集中t1细胞的总数,以获得估计的凋亡率。 我们独立地对所有时间点进行了上述计算,这些时间点包括moscot.time和clTOME。 我们通过结合来自不同出版物的信息选择了目标凋亡范围,这些出版物包括18、107、108、109、110、111。 我们为原肠胚形成前阶段选择的目标凋亡范围是10-15%的细胞发生凋亡,原肠胚形成阶段是4-6%,器官发生阶段是2-4%。 对于时间对E8.5a-E8.5b,尽管没有实际的时间流逝,但实验方法发生了转变,我们旨在实现相对较高的凋亡率,即10-40%,以便纠正采样偏差。 Correlating predicted growth rates with gene-set-based growth rates 将预测生长率与基于基因集的生长率相关联
Para_01 为了验证预测的生长速率,我们通过scanpy中的scanpy.tl.score_genes根据标记基因表达计算细胞周期得分并与之相关联。 简而言之,scanpy中的基因评分实现遵循Seurat(版本1)的原始建议:它通过对提供的基因集的基因平均值进行归一化处理,该归一化处理基于参考基因集的平均表达量。 为了进行此比较,我们将边缘均匀初始化,以便我们的算法不了解生长速率,并且我们可以使用这些信息进行验证。 Para_02 我们应用这一策略对一个不同的数据集进行了研究,该数据集包含了重编程的小鼠胚胎成纤维细胞。 这个数据集更适合用于生长速率的比较,原因有两个。 首先,我们用于根据单细胞RNA测序数据评分细胞周期的基因集是专门为小鼠成纤维细胞和造血干细胞设计的。 将这个基因集应用于小鼠胚胎发生图谱使得造血谱系的得分始终高于其他细胞类型,这与之前的生物学发现相矛盾。 其次,小鼠胚胎发生图谱代表了一种体内环境,在这种环境下,每个时间点对应不同的个体,因此导致了不同时间点之间细胞类型比例的强烈变化。 这些变化不是由细胞的生长和死亡驱动的,而是由细胞取样效应引起的。 特别是,胚外组织受到了大的系统性取样偏差的影响,这很可能是由于样品处理的变化所导致的。 Para_03 相比之下,小鼠胚胎成纤维细胞重编程数据集更适合我们的基因集,并且由于是在体外环境下,包含的由细胞取样引起的细胞类型频率偏差较少。我们使用0.0005的熵正则化参数和(\tau_a=0.98)以及(\tau_b=0.99995)的不平衡参数运行了moscot.time。 然后我们在同一数据集上计算了clTOME生长速率。 我们通过使用scanpy.tl.score_genes方法,基于为小鼠成纤维细胞和造血干细胞确定的基因集,计算了细胞周期评分。 Comparison in terms of driver-gene correlations 在驱动基因相关性方面进行比较
Para_01 为了进一步评估由moscot.time和TOME预测的细胞层级耦合,我们认为祖先概率与已知驱动基因之间的高相关性表明该方法的成功。 因此,对于主体文本中描述的细胞状态,我们计算了使用moscot.time和clTOME预测的祖先分布(参见"跨时间映射细胞的moscot.time"部分)。 为了排除无关分化事件中的驱动基因的影响,我们将相关性计算限制在已知祖先进化群体内。 对于每个被提取的细胞状态,我们精心编制了一份已知驱动基因列表(补充表4),过滤该列表仅包含相应时间点上的高度可变基因,并使用scVI的解码器输出get_normalized_expression来推断它们的表达。 在过滤到高度可变基因后,我们保留了36个基因用于确定性内胚层,18个基因用于卵黄囊,39个基因用于第一心场,以及106个基因用于胰腺上皮。 我们使用scipy.stats.spearmanr计算了这些推断的表达值与预测的祖先分布之间的Spearman相关系数值。 Metacell analysis 类细胞分析
Para_01 加速映射计算的另一种方法是将细胞聚合成元细胞。 为了调查这种方法的性能,我们在E9.5细胞上使用流行的Metacell-2算法21计算了元细胞。 然而,我们发现对于罕见的细胞状态原始生殖细胞(30个细胞或占总体的0.03%),没有创建元细胞,这使得推断这种细胞状态的祖细胞和祖先细胞变得不可能。 错误!!! - 待补充
Moscot.time on multimodal pancreas development Moscot.time 关于多模态胰腺发育
Dataset generation 数据集生成
Para_01 来自NVF纯合子小鼠的胚胎胰腺被收集并合并在一起(第一次实验(exp-1)中E14.5阶段收集了8个胰腺,E15.5阶段收集了11个胰腺;第二次实验(exp-2)中E15.5阶段收集了10个胰腺,E16.5阶段收集了10个胰腺)。 向样品中加入胰蛋白酶(0.25%),在冰上处理5分钟,然后在37°C下孵育10分钟。 单细胞样品在4°C下以1700转/分钟(290g)离心5分钟。 去除上清液后,对细胞进行计数。 接下来,使用5 µl大鼠IgG2a K同型对照(eBioscience,12-4321-42)和抗小鼠CD326(EpCAM)PE(eBioscience,12-5791-81)处理1 × 106个细胞(总体积100 µl)。 样品在4°C下染色30分钟后,用DAPI染色以检测死细胞。 清洗两次后,重新悬浮在FACS缓冲液(PBS、1% BSA和0.5 mM EDTA)中,将单细胞样品加载到FACS分析仪中。 使用的门控策略如下:主要群体>单细胞>活细胞(DAPI–)> EpCAM+(PE+)和NGN3+(FITC+)/NGN3–(FITC–)细胞。 分选后的细胞按2:1的比例(EpCAM+NGN3+: EpCAM+NGN3–)合并,并立即用于核分离。 Para_02 为了分离细胞核,我们采用了一种从10x Genomics改良的低输入细胞核分离方案。简而言之,分选后的细胞用1毫升PBS+1%BSA洗涤一次,基于台盼蓝染色计数后离心。 随后,洗涤过的细胞团块重新悬浮在预冷的裂解缓冲液(每份样品50微升)中,并在冰上放置4分钟。 然后加入洗涤缓冲液(每份样品500微升),并将细胞核离心。 为了逐步从洗涤缓冲液过渡到稀释的细胞核缓冲液,细胞被一次在1:1的洗涤缓冲液和稀释的细胞核缓冲液混合物中洗涤,随后再用纯稀释的细胞核缓冲液洗涤一次。 洗涤过的分离出的细胞核随后重新悬浮在7-10微升稀释的细胞核缓冲液中,在经过质量控制和计数后,立即用于单细胞多组学文库的制备,目标回收细胞数为10,000个。 Para_03 文库使用Chromium Next GEM单细胞多组学ATAC + 基因表达试剂盒(10x Genomics, 1000283)根据制造商的说明制备。 文库遵循10x Genomics的建议在Illumina NovaSeq6000平台上进行测序。 两个实验的原始读数使用Ensembl发布102注释的GRCm38小鼠基因组共同对齐,并使用10x Genomics CellRangerARC管道(v.2.0.2)进行预处理以供下游分析使用。 Preprocessing 预处理
Para_01 我们独立地对基因表达和染色质可及性样本进行了预处理。 每个样本从CellRanger输出中独立获取峰,并随后合并。 Peaks是一个专业术语,在这里保持原样。 Para_02 关于基因表达,所有在E14.5时期线粒体基因比例高于0.025或在E15.5时期线粒体基因比例高于0.02的细胞在exp-1中被移除。 对于exp-2,在E15.5时期线粒体基因比例高于0.015的所有细胞被移除,而E16.5时期的阈值设置为0.02。 此外,对于exp-1,在E14.5时期,少于4,000计数或多于30,000计数的细胞被移除;在E15.5时期,少于5,000计数或多于40,000计数的细胞被移除。 对于exp-2,两个时间点的较低阈值都被设定为3,000,而在E15.5和E16.5时期,至少有60,000和70,000总基因计数的细胞被移除。 对于E14.5(exp-1),所有少于2,300个基因表达的细胞被过滤掉。 E15.5(exp-1)、E15.5(exp-2)和E16.5(exp-2)的相应阈值分别设定为2,700、2,000和2,000。 Para_03 关于ATAC模态,所有E14.5(exp-1)中的细胞,如果核小体信号低于0.35或高于1.75,则被移除。 所有转录起始位点富集分数低于2.5或高于7.5的细胞也被过滤掉。 如果细胞的总开放染色质区域计数低于4,000或高于150,000,则也予以移除。 类似地,E15.5(exp-1)、E15.5(exp-2)和E16.5(exp-2)的最低核小体信号分别设置为0.3、0.35和0.25。 而核小体信号的上限分别为1.75、1.5和1.4。 此外,E15.5(exp-1)、E15.5(exp-2)和E16.5(exp-2)的转录起始位点富集分数的下限分别设置为2.75、2.5和2.5。 而上限则分别设置为10.5、8和7.5。 对于E15.5(exp-1)、E15.5(exp-2)和E16.5(exp-2),总峰计数的下限均设为4,000。 而上限则分别设为100,000、160,000和170,000。 Para_04 在合并两个样本后,过滤掉了在少于20个细胞中被检测到的基因,最终剩下20,244个基因。 Para_05 双细胞是使用多种双细胞检测方法的平均预测来识别的。 我们使用了Scrublet、scDblFinder、DoubletDetection、scds、SOLO和DoubletFinder。 为了基于ATAC计数识别双细胞,我们使用了AMULET。 只要七种方法中有至少三种方法将一个细胞评定为双细胞,我们就认为该细胞是双细胞。 总共,在样本E14.5(exp-1)中鉴定出12.60%的双细胞,在E15.5(exp-1)中鉴定出10.73%的双细胞。 在E15.5(exp-2)中鉴定出16.68%的双细胞,在E16.5(exp-2)中鉴定出15.11%的双细胞。 Para_06 随后,在 MultiVI62 嵌入中重复进行了聚类分析,并去除了包含大量已识别双重细胞的聚类。 , Cell-type annotation 细胞类型注释
Para_01 为了构建加权最近邻图,需要对两种模态进行嵌入。 因此,在对对数转换后的基因表达数据执行PCA(50维)之前,使用SCTransform123对计数数据进行了标准化,并且去除了细胞周期基因和环境基因。 使用DropletUtils124识别环境基因。 ATAC数据通过词频-逆文档频率(tf-idf)归一化处理,然后使用Signac进行奇异值分解,计算前50个奇异分量。 由于与测序深度高度相关,去除了第一和第五个分量。 计算了GEX和ATAC各自的嵌入后,我们使用MUON125构建了一个加权最近邻图,并用于多模态、无监督聚类。 除非另有说明,这也是我们计算UMAP所用的图。 Para_02 注释是基于先前研究34,42,113,126,127,128(补充表9)报道的标记基因表达以及使用scanpy.tl.score_genes_cell_cycle计算的增殖群体的细胞周期评分进行的。 值得注意的是,我们发现了一个从Ngn3High群体分支出来的集群,我们发现它表达了与之前描述的Fev+ epsilon集群相似的基因34。 事实上,该研究34中报告的集群以及我们在新数据集中发现的集群都没有显著高的Fev表达(补充图28)。 因此,我们将这个集群标记为epsilon前体细胞。 Para_03 为了达到如图3e所示的细胞类型的更精细分辨率,在相同的邻域图上进行了亚聚类分析(结合了两种模态)。 , The moscot.time model moscot.time 模型
Para_01 我们通过使用scanpy.pp.nearest_neighbors计算的30个最近邻图(基于MultiVI嵌入计算的热核扩散距离)(补充说明5)定义了E14.5和E15.5之间的最优传输问题的成本矩阵。因此,这个图是在与用于无监督聚类的不同嵌入上构建的,以减少对单一嵌入的偏差。 这样,该图是基于不同的嵌入构建的,以减少偏向某一特定嵌入的偏差。 Para_02 基于上述描述构建的加权最近邻图,运行了两个moscot模型。首先,在完整数据集上运行了一个模型。 使用默认参数运行了moscot.time模型,但通过设置({\tau }{a"}={\tau } {b"}=0.99)引入了一点不平衡性。 具体来说,正则化参数({\rm{\epsilon }})被设置为10–3,并且成本矩阵按其均值进行了缩放。 为了保证收敛性,迭代次数增加到107次。 选择了均匀边缘分布,因为高度增殖的导管细胞的大量存在会使较少数量的细胞类型的影响被边缘化。 还需要注意的是,该数据集是FACS分选的;因此,最初测序细胞的比例高度偏倚,并不能反映真实的细胞类型分布。 这也导致最终模型无法预测跨越一天(墙钟时间)的后代(祖先)。 实际上,发育过程的方向性得以保持,而其幅度并不反映真实生物进展的程度。 Para_03 为了分析内分泌分支,OT解决方案是在一个仅包含内分泌细胞(alpha、beta、delta和epsilon)及其祖细胞(标记为Fev+ alpha、Fev+ beta、Fev+ delta、epsilon祖细胞、Fev+、Ngn3High、Ngn3High循环或Ngn3Low)的缩减数据集上计算的。 同样地,选择了统一的边缘分布,因为从TemporalProblem.score_genes_for_marginals获得的增殖和凋亡评分几乎是恒定的。 使用了moscot提供的标准参数来计算OT解决方案。 Studying the influence of the embedding and cost function 研究嵌入和成本函数的影响
错误!!! - 待补充
Para_02 我们将计算出的运输矩阵聚合到细胞类型层面(A),并连续列归一化以获得祖先的概率。 然后我们计算了聚合运输矩阵A与参考聚合运输矩阵B之间的Sinkhorn散度,这是我们用于分析胰腺数据集的moscot耦合。 Sinkhorn散度的成本被选择为二进制距离,也就是说,只包含0和1的条目。 我们也考虑了以下转换:delta来自Fev+ delta;epsilon来自epsilon祖细胞;epsilon来自Fev+ delta;Fev+ delta来自epsilon祖细胞;以及beta来自Fev+ beta(作为生物学上已知的转换)。 我们将这些转换的值与通过外部耦合获得的转换值进行了比较。 Para_03 我们还报告了delta细胞直接来源于Ngn3Low的概率,这是一个在生物学上不切实际的转变。这里我们看到,成本函数的选择非常重要,但这与嵌入无关。 也就是说,当使用平方欧几里得成本或余弦成本时,我们观察到这种转变的比例显著,而使用测地线成本时,我们没有观察到这种转变。 Driver feature analysis with moscot.time 使用 moscot.time 进行驱动特征分析
Para_01 我们通过将回拉分布的密度与特征值相关联来计算驱动器特征。 此外,在分析从epsilon细胞向alpha细胞的转变时,我们将Fev+ alpha群体(alpha细胞的主要祖细胞类型)排除在考虑的细胞类型之外。 也就是说,我们将它们的祖先概率设置为零。 这有助于识别在将回拉分布与处理后的基因表达相关联时,在epsilon细胞中特别激活的基因。 Marker regions of chromatin accessibility 染色质可及区域的标记
Para_01 为了识别染色质可及性的标记区域,使用了 Seurat 中提供的 FindMarkers 进行 Wilcoxon 检验。 该检验采用了默认设置,并且考虑的细胞类型是相对于所有剩余的细胞类型进行的(限定为内分泌细胞和内分泌祖细胞;也就是说,Ngn3Low、Ngn3High、epsilon 祖细胞、Fev+、Fev+ delta、Fev+ alpha、Fev+ beta、alpha、beta、delta 和 epsilon)。 Motif analysis 基序分析
Para_01 基序数据从cisBP130下载。 位置权重矩阵以及相应的可视化和元数据在2023年3月1日通过按物种(小家鼠)过滤后进行了批量下载。 cisBP包含来自实验测量的结合活性和推断的结合活性的数据(例如,来自其他物种)。 具有DNA结合域氨基酸相似性超过特定阈值的转录因子(为每种DNA结合域类别分别定义,并由cisBP提供)也被视为结合候选物。 Para_02 我们定义一个转录因子与一个基序有联系,如果它是直接测量的或推断出来的,并且具有足够高的DNA结合域氨基酸相似性。也就是说,由cisBP报告为这样。这样一来,一个基序可以与多个转录因子有联系,一个转录因子也可以与多个基序有联系。 这是第二句话的中文翻译。 Para_03 为了在单细胞水平上获得基序分数,我们使用Signac提供的API运行了chromVAR。实际上,我们调用了AddMotifs,然后运行了RunChromVAR。 为了使用moscot获得标记基序,我们考虑了基因表达的时间顺序和基序的活性。 值得注意的是,moscot附带了一张不同物种(人类、小鼠和果蝇)的转录因子列表,这些列表来自SCENIC+数据库131。 因此,我们利用moscot计算驱动特征的能力来计算驱动转录因子。 此外,我们基于ChromVAR评分进行了差异基序活性检验(使用scanpy的rank_genes_groups进行Wilcoxon检验)。 随后,我们通过结合这两种信息来源来识别标记基序。 因此,我们只保留了那些有相关转录因子的标记转录因子。 Cell-cycle analysis 细胞周期分析
Para_01 我们使用了胰腺内分泌生成的数据集。 由于我们研究了增殖性导管细胞的细胞周期,我们去除了未成熟和成熟的胰腺细胞。 连续地,我们在共享的MultiVI嵌入上计算了E14.5和E15.5之间以及E15.5和E16.5之间的最优传输(OT)耦合。 连续地,我们将所有增殖性导管细胞的细胞周期阶段进行了分配,如建议的那样。 通过将运输矩阵与细胞周期阶段注释进行聚合,并进行连续的行规范化,我们获得了从细胞周期阶段i到细胞周期阶段j的转换概率(有序的细胞周期阶段定义为G1S、S、G2M、M和MG1)。 请注意,我们没有考虑包括非循环细胞的转换。 因此,聚合矩阵A代表了一个细胞在早期时间点处于某个细胞周期阶段时,过渡到后期时间点的特定细胞周期阶段的概率。 正确的(或错误的)方向性意味着细胞周期阶段i的细胞更有可能(或不那么可能)过渡到细胞周期阶段i + 1而不是过渡到细胞周期阶段i – 1。 在这里,我们假设如果i = 5,则i + 1 = 1(因此从MG1(索引5)到G1S(索引1)的过渡),如果i = 1,则i – 1 = 5(因此从G1S(索引1)到MG1(索引5)的过渡)来闭合循环。 换句话说,我们认为从i到i – 1的转换是不好的(错误的方向),而从i到i + 1的转换是好的(正确的方向)。 因此,我们能够通过对5x5的聚合运输矩阵A进行加法运算,得到分数-Ai,i-1(注意负号)和+Ai,i+1(注意正号)来构建一个分数。 如果分数为正,这意味着方向是正确的(因为更多的细胞从阶段i过渡到i + 1,而不是从阶段i过渡到阶段i – 1)。 我们分别计算了E14.5至E15.5和E15.5至E16.5的转换分数,并进行了置换检验(通过置换A的所有非对角线元素),使用了10,000次置换。 Trajectory inference with different trajectory-inference methods 使用不同的轨迹推断方法进行轨迹推理
Para_01 我们使用扩散伪时间、scVelo、veloVI、MultiVelo、CytoTrace和CellRank中的ConnectivityKernel来预测胰腺内分泌发生数据集中的轨迹。 由于我们对内分泌细胞的轨迹感兴趣,我们将数据集过滤为内分泌细胞及其前体细胞。 然后,我们对每个相应的轨迹推断核应用了CellRank中的GPCCA估计器。 为了计算命运概率,我们使用了compute_fate_probabilities和aggregate_fate_probabilities来绘制命运概率和聚合的细胞类型到细胞类型的转换矩阵。 我们使用plot_projection方法生成流嵌入图。 我们在所有方法中使用了提供的所有默认参数,并且仅将VeloVI中的max_epochs增加到了50。 当使用图时,我们使用上述描述的WNN图。 对于使用moscot构建RealTimeKernel,我们将ConnectivityKernel的权重设置为0.001,以增强moscot的影响并减弱ConnectivityKernel的影响。 我们指出,使用CellRank计算的转移概率依赖于与使用moscot计算的转移概率不同的过程。 Labelling of the quality of trajectory-inference methods 轨迹推断方法的质量标注
Para_01 为了评估轨迹推断方法的可靠性,我们通过关注已充分研究的α细胞和β细胞谱系的命运概率来评估它们的性能。 实际上,我们考虑了α细胞、β细胞、Fev+ α细胞和Fev+ β细胞的命运概率,其中α、β、δ和ε可能是可能的谱系。 我们认为如果最高命运概率是最有可能的生物命运(即α→α,Fev+ α→α,β→β,Fev+ β→β),那么转变是正确的。 最后,如果四个转变都是正确的,我们将该方法标记为绿点;如果三个正确,则标记为橙点;对于其他情况,则标记为红点。 Assessment of the consistency with the fate probabilities of moscot 对 moscot 的命运概率一致性评估
Para_01 为了评估模式和moscot之间预测的一致性,我们计算了使用CellRank输出的聚集命运概率之间的皮尔逊相关系数。 Diffusion pseudotime based on different modalities 基于不同模态的扩散伪时间
Para_01 为了研究模式对于扩散伪时间预测的影响,我们基于基因表达的PCA和ATAC模态的LSI构建了一个图(使用scanpy.pp.neighbors,默认参数)。 在计算扩散伪时间时,我们使用了muon中的加权最近邻实现,默认参数。 Interpretability of the transport map with Sparse Monge 具有稀疏蒙古映射的传输映射的可解释性
Para_01 我们应用稀疏蒙古算法(如上所述)处理了 E15.5 和 E16.5 时间点的内分泌祖细胞。 我们使用 scanpy.pp.highly_variable_genes 选择了最具变异性的基因,默认配置下,即 flavor='seurat'。 这导致了 2,551 个高变异性基因被选择。 然后,我们在归一化和对数变换空间中,使用弹性 L1 成本和正则化因子 γ=10(公式(13)),将 E15.5 时期的细胞运输到 E16.5 时期的细胞。 正如之前描述的那样47,我们通过确定细胞位移是否大于 10–6 来识别每个单细胞的相关基因。 随后,我们将这些特定于细胞的基因集称为重要基因。 我们通过计算具有该基因作为重要基因的细胞比例来识别每种细胞类型中最相关基因。 为了恢复命运决定中的变异性,我们构建了一个最近邻图(k = 50),并计算了所考虑细胞的重要基因与邻域的重要基因之间的杰卡德相似性。 我们计算这个分数的补数,以获得一个关于不相似性的概念。 Gene KO of NEUROD2 and iPS cell differentiation to pancreatic endocrine cells NEUROD2的基因敲除和ips细胞向胰腺内分泌细胞的分化
Para_01 所有统计分析均使用 GraphPad Prism 10 进行了一维方差分析。 , In vitro differentiation of iPS cells to pancreatic endocrine cells iPS细胞体外分化为胰岛内分泌细胞
Para_01 两个纯合的NEUROD2基因敲除、核H2B-Venus报告子(NEUROD2nVenus/nVenus)iPS细胞克隆和杂合的hiPSC-INS-T2A-H2B-Cherry报告子(INSCherry/WT)iPS细胞系(作为对照)被使用。 我们使用了多步骤分化方案来体外将iPS细胞分化为内分泌和胰岛细胞,该方案包括阶段(S)定型内胚层(S1),原始肠管(S2),胰腺前体1(PP1)(S3),PP2(S4)和内分泌谱系(S5,6),如先前所述138,139。 C肽-mCherry报告人iPS细胞系(HMGUi001-A-8)和NEUROD2nVenus/nVenus iPS细胞克隆(HMGUi001-A-42)被使用。 所有细胞系都定期检测以确保它们没有支原体污染。 样本量是根据可用的实验数据确定的。 使用的抗体有山羊抗生长抑素,1:300,多克隆(D-20)(Santa Cruz Biotechnology,SC-7819)和小鼠抗胃泌素,1:250,单克隆(2F4)(Santa Cruz Biotechnology,SC-293422)。 RNA extraction and qPCR analysis RNA提取和qPCR分析
Para_01 总RNA是从样品中使用miRNeasy mini试剂盒(Qiagen)提取的。 随后,使用SuperScript Vilo cDNA合成试剂盒(Thermo Fisher Scientific)根据制造商的说明进行逆转录。 预设计的TaqMan探针(Life Technologies)用于qPCR分析(序列见补充表1)。 每个样品的反应混合物含有20纳克cDNA,4.5微升无核酸酶水,5微升TaqMan高级混合物(Life Technologies)和0.5微升TaqMan探针(Life Technologies)。 反应在QuantStudio 7 Flex仪器(Thermo Fisher Scientific)上运行。 GAPDH被用作参照基因以进行标准化。 为了保持数据的分布,并便于假设方差相等的统计分析,KO样本的Ct值被标准化到平均对照组值140。 来自单次qPCR运行中的独立样品的数据一起进行了分析(更多信息见源数据)。 使用的引物如下:GAPDH使用Hs02758991_g1,SST使用Hs00356144_m1,GHRL使用Hs01074053_m1,INS使用Hs02741908_m1,GCG使用Hs00242160_m1,HHEX使用Hs01031536_m1。 Immunostaining and imaging 免疫染色和成像
Para_01 冷冻切片制备、固定和免疫染色按照先前所述的方法进行45。 使用的初级抗体有:生长抑素(D-20)山羊多克隆抗体(Santa Cruz sc-7819)和生长激素释放肽小鼠单克隆抗体(Santa Cruz sc-293422)。 使用Leica DMI 6000显微镜和LAS AF软件拍摄照片。 使用LAS AF和ImageJ软件程序对图像进行分析和量化。 Spatial analysis 空间分析
Benchmarking moscot.space.mapping across a range of spatial datasets 基准测试 moscot.space.mapping 跨越一系列空间数据集
Para_01 我们对比了 moscot 在映射问题上的表现,与两种最先进的方法,Tangram23 和 gimVI24,正如 scVI 工具60中所实现的那样。 我们使用了之前发布的数据集25。 从这些数据集中,我们选择了能够重新处理的数据集,这导致我们考虑了14个用于基准测试的数据集。 此外,与原始基准不同,我们没有使用单细胞数据集作为参考,因为我们不认为这样的数据能代表一个忠实的地面真实情况来比较这些方法。 因此,我们将空间数据集分割,使得数据点的50%被用作单细胞参考,而剩余的50%则作为空间数据。 同时,我们明确保持输入数据类型符合模型需求。 因此,我们对 moscot 和 Tangram 进行了归一化和 log1P 变换,而对于 gimVI 则保留了原始未归一化的计数。 如果数据集中基因总数超过2,000,则随机保留100个基因,否则保留10个基因。 我们对剩余的基因训练模型,并使用斯皮尔曼相关系数评估性能。 我们报告了使用三个随机种子(包括数据集分割和初始化及训练过程中的随机种子)得到的斯皮尔曼相关系数的平均值。 对于某些数据集,由于时间复杂度(我们为每种方法设定了最多5个GPU小时的预算)或模型错误(例如,在训练和推断数据之间无法匹配基因标识符),Tangram 或 gimVI 无法运行。 Para_02 具体来说,我们在以下参数上进行了扫描: [ul]- moscot: epsilon entropy regularization parameter, alpha interpolation parameter between W term and GW term, and tau_a unbalancedness term for the spatial dataset (source). For the cost, we tried cosine and squared Euclidean cost for the linear and quadratic term and joint_attr; that is, the representation to use for the linear term. We assessed both PCA and gene expression on a common set of genes present in both spatial and single-cell datasets. - Tangram: learning rate and number of epochs. - gimVI: number of epochs and number of latent dimensions.
Para_03 我们还报告了每种算法在各个数据集和种子上的内存和时间复杂度。 所有实验都在赫姆霍兹集群的GPU上运行(混合使用V100和A100 GPU)。 Spatial correspondence 空间对应关系
Para_01 空间对应关系的计算方法如下。首先,我们在数据集中的所有数据点之间计算了n个不断增加的空间距离(欧几里得)阈值。 然后,在每个阈值水平上,我们计算了所有距离低于所选阈值的所有斑点中所有基因之间的基因表达相似性(欧几里得距离)。 空间对应关系随后被计算为基因表达相似性和空间距离阈值之间的皮尔逊相关系数。 该计算作为moscot映射问题的方法实现。 Moscot.space.mapping on the liver Moscot.space映射到肝脏
Para_01 我们应用了 moscot.space 的映射问题到从 https://vizgen.com/data-release-program/ 下载的小鼠肝脏数据集上,该数据集来自 Vizgen MERSCOPE。 我们遵循标准的 scanpy 和 squidpy 处理流程对数据集进行了处理。 对于单细胞参考,我们下载了 CITE-seq 数据集,该数据集来自 https://www.livercellatlas.org/,并在先前的一项研究中有所报道(26)。 我们使用映射问题的方法如下:对于线性项,我们使用了 336 个常见基因;而对于二次项,我们使用了单细胞参考数据集的基因表达主成分分析(PCA),以及将空间数据集的基因表达 PCA 叠加到空间坐标上。 然后,我们通过计算蛋白质表达的重心投影(方程(16))到空间数据集来执行基因表达和蛋白质填充。 同样的重心投影方法也被用来将单细胞参考中的细胞类型注释转移到空间数据集。 Moscot.space.mapping on spatial ATAC–seq data Moscot.space.mapping 对空间ATAC-seq数据的映射
Para_01 为了基准测试是否利用多模态表示能改善模态映射,我们考虑了一个以前处理过的数据集,该数据集用于人类扁桃体的联合多模态RNA和ATAC分析。 该数据集包括一块人类扁桃体活检样本的单张切片,使用了DBiT-seq技术的修改版本进行分析,该技术能够对同一捕获位置的染色质可及性和基因表达进行分析。 总的来说,该数据集包含2500个独特的捕获位置。 我们使用标准的Scanpy工作流程进行了特征选择、主成分分析(PCA)和UMAP降维。 我们将数据集随机分为两部分,将第一部分作为代理单细胞数据集,其中包括基因表达和染色质可及性信息,第二部分作为代理空间数据集,其中包括基因表达和空间坐标。 然后,我们着手评估是否在映射问题的二次项中利用额外的模态可以改善ATAC情况下染色质可及性的预测。 我们仅使用空间数据集中的ATAC信息进行评估,该信息是通过单细胞数据集的重心投影预测的峰值可及性使用皮尔逊相关系数测量的。 特别是,我们评估了作者提供的七个聚类的前十个标记峰的平均相关性,这导致总共70个可及性峰。 所有实验都在Helmholtz集群上的GPU上运行(V100和A100 GPU的混合)。 基准测试使用Hydra运行。 Benchmarking of moscot.space.alignment on simulated data moscot.space.alignment在模拟数据上的基准测试
Para_01 我们针对两个最先进的对齐方法:PASTE4和GPSA30,对moscot的对齐问题进行了基准测试。 我们在所有方法中选择了相同的计算预算;也就是说,12组独特的超参数: [ul]- Moscot: epsilon (entropy regularization parameter) and alpha (interpolation parameter between W term and GW term). - PASTE: alpha (interpolation parameter between W term and GW term) and norm (scaling of the cost matrix). - GPSA: kernel (kernel for the Gaussian Process), n_epochs (number of epochs) and lr (learning rate).
Para_02 由于无法在GPU上运行GPSA,我们所有方法都在CPU上运行。 我们基于GPSA出版物中描述的数据生成方式生成了四个合成数据集。 简而言之,生成了从随机正态分布中抽取的样本来构建一个按网格排列的合成基因表达文件。 然后,通过原始数据集的0.7、0.8和0.9的比例随机子采样点,使得源数据集和目标数据集中的总点数不匹配。 这种方法类似于先前使用的基准设置。 然而,为了使三种方法具有可比性,我们在PASTE和moscot中使用了空间坐标相对于耦合的重心投影(公式(16))。 由于实验样本量较小,我们以全秩模式运行了moscot(而不是低秩模式)。 较大的数据集,例如正文分析的数据集,对于PASTE和GPSA来说都将过于庞大。 Moscot.space.alignment on mouse brain coronal sections Moscot.space 对齐小鼠脑冠状切片
Para_01 我们应用了 moscot 的对齐问题来处理来自 Vizgen MERSCOPE(https://vizgen.com/data-release-program/)的大规模 MERFISH 数据集。 具体来说,是两部分小鼠冠状脑切片。 我们对每个切片中的三个不同小鼠样本进行了对齐。 我们使用 moscot 的对齐问题在仿射模式下(公式(17))进行了第一次对齐。 因此,三张切片中的两张被对齐到了第三张作为参考。 此外,我们在 FGW 对齐坐标上进行了一次第二次对齐,使用了 W 类型的问题,以获得改进的扭曲对齐。 这在低秩设置中证明是有用的。 我们对所有冠状切片的三重组合执行了相同的操作。 Data availability Para_01 小鼠胚胎发生图谱可在http://tome.gs.washington.edu获取。 小鼠肝脏CITE-seq数据可在https://www.livercellatlas.org/获取。 Vizgen MERSCOPE肝脏和大脑冠状切片数据集可在Vizgen公共数据集发布网站https://vizgen.com/data-release-program/获取。 用于基准测试空间映射问题的数据集来自先前的一项研究。 小鼠胚胎发生的时空图谱(MOSTA)可在https://db.cngb.org/stomics/mosta/获取。 果蝇的时空数据集可在https://db.cngb.org/stomics/flysta3d/获取。 单细胞RNA测序数据集可从基因表达全景数据库(GEO;https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE132188)获取。 胰腺多组学数据可从GEO(登录码GSE275562)获取。 Code availability Para_01 moscot 软件包可在 https://moscot-tools.org 获取,其中包含文档、教程和示例。 我们的分析结果可以通过 GitHub (https://github.com/theislab/moscot-framework_reproducibility) 复现的代码获取。 基准测试实验的复现代码也可以通过 GitHub (https://github.com/theislab/moscot_benchmarks) 获取。