4.1 Pseudotemporal ordering¶
Motivation¶
单细胞测序技术可以对生物组织进行高分辨率的测量,有助于破译和理解细胞异质性以及生物过程的动态变化。相关研究包括量化细胞命运以及鉴定驱动该过程的基因。
然而,经典的单细胞 RNA 测序 (scRNA-seq) 协议在测序过程中会破坏细胞,因此无法追踪细胞的发育过程,例如无法追踪基因表达谱随时间的变化。尽管最近的技术进步允许顺序记录转录组,但它们在实验上具有挑战性,并且目前无法扩展到更大的数据集上。因此,我们只能根据测得的静态数据来估计潜在的动态过程。
尽管传统上样品只取自单个实验时间点,但仍可以观察到多种细胞类型。这种多样性源于生物过程的异步性,因此可以观察到发育过程的不同阶段。 领域内新兴的“轨迹推断 (TI)” 旨在通过对观察到的细胞状态进行排序来重建发育图景。这项任务是通过将离散注释映射到连续域(即所谓的“伪时间”)来实现的,从而使细胞状态沿着发育方向排列。
伪时间根据细胞在发育过程中的相应阶段对其进行相对排序。发育较早期的细胞分配较小的值,成熟细胞分配较大的值。例如,在研究骨髓样本时,造血干细胞会分配较低的伪时间,而红细胞会分配较高的伪时间。对于单细胞 RNA 测序数据,分配伪时间是基于细胞的转录组谱。此外,构建伪时间通常还需要指定一个初始细胞,即整个过程的起点。
Pseudotime construction¶
伪时间构建通常遵循以下通用工作流程:首先,将超高维的单细胞数据投影到低维空间。这种做法的合理性在于观察到动态过程发生在低维流形上 {cite}pt:Wagner2016
。实际应用中,伪时间方法可能依赖于主成分分析(例如 Palantir {cite}Setty2019
)或扩散成分(例如扩散伪时间 (DPT) {cite}Haghverdi2016
)。接下来,伪时间将基于以下原则之一进行构建。
聚类方法:首先将观察结果聚类,然后识别这些簇之间的连接。可以通过对簇排序来构建伪时间。我们将这种方法称为“聚类方法”。常用的聚类算法包括 k-means {cite}
Lloyd1982
, {cite}MacQueen1967
、Leiden {cite}pt:Traag2019
或层次聚类 {cite}Mueller2011
。簇的连接可以基于相似性,也可以通过构建最小生成树 (MST) {cite}Pettie2002
来实现。图方法:首先找到低维表示的观察结果之间的连接。此过程定义了一个图,基于该图定义簇以及排序。例如,PAGA {cite}
Wolf2019
将图划分为 Leiden 簇并估计它们之间的连接。直观上,这种方法在较低分辨率下分析数据的同时保留了数据的全局拓扑结构。因此,提高了计算效率。流形学习方法:与“聚类方法”类似,但通过使用主曲线或图来估计潜在轨迹来定义簇之间的连接。主曲线在高维空间中找到连接细胞观察结果的一维曲线。Slingshot {cite}
Street2018
是这种方法的典型代表。概率框架:将转移概率分配给有序的细胞对。每个转移概率量化参考细胞是另一个细胞祖先的可能性。这些概率定义了用于定义伪时间的随机过程。例如,DPT 被定义为随机游走连续状态之间的差异。相比之下,Palantir {cite}
Setty2019
将轨迹本身建模为马尔可夫链。虽然这两种方法都依赖于概率框架,但它们都需要指定一个根细胞。伪时间本身相对于该细胞进行计算。
伪时间推断 (TI) 是一个研究成熟的领域,提供了丰富的相关方法。为了将合适的方法应用于单细胞数据集的分析,首先需要理解生物过程本身。这种理解尤其包括过程的性质,例如它是线性的、循环的还是分支的。类似地,同一个数据集中存在正交过程也会限制可应用的 TI 方法。dynguidelines {cite}Deconinck2021
提供了算法及其特性的详尽概述,可以帮助识别合适的工具。
下游任务与展望¶
尽管伪时间和轨迹推断 (TI) 已经能够提供宝贵的见解,它们通常只是更精细分析的垫脚石。例如,识别终状态是一个经典的生物学问题,可以通过 TI 进行研究。类似地,谱系分叉和驱动细胞命运决定的基因也可以基于 TI 和伪时间进行识别。能够回答哪些问题以及如何找到答案通常取决于具体的方法。例如,Palantir 通过识别其构建的马尔可夫链的吸收状态来识别终状态。
尽管轨迹推断 (TI) 的成功已经得到充分验证,并且许多方法被提出,但随着测序技术的进步,我们可以利用新的信息来源。例如,ATAC-seq {cite}Buenrostro2015
、CITE-seq {cite}pt:Stoeckius2017
和 DOGMA-seq {cite}pt:Mimitou2021
等技术可以测量除转录组之外的其他调控模式。谱系追踪 {cite}Weinreb2020
和代谢标记 {cite}Erhard2019
, {cite}Battich2020
, {cite}Qiu2020
, {cite}Erhard2022
甚至可以提供给定细胞的 (可能) 未来状态。因此,未来的 TI 工具将能够包含更多信息来更准确、更稳健地估计轨迹和伪时间,并回答新的问题。例如,RNA 速度 {cite}LaManno2018
, {cite}Bergen2020
, {cite}Bergen2021
是一种利用未剪切和剪切 mRNA 来推断经典静态快照数据之外的定向动态信息的技术之一。
成人骨髓的伪时间推断¶
为了展示如何构建伪时间以及比较不同伪时间的结果,我们接下来将分析来自成人骨髓的单细胞 RNA-seq 数据集 {cite}Setty2019
。
Environment setup¶
from pathlib import Path
import scanpy as sc
General settings¶
DATA_DIR = Path("../../data/")
DATA_DIR.mkdir(parents=True, exist_ok=True)
FILE_NAME = DATA_DIR / "bone_marrow.h5ad"
Data loading¶
adata = sc.read(
filename=FILE_NAME,
backup_url="https://figshare.com/ndownloader/files/35826944",
)
adata
AnnData object with n_obs × n_vars = 5780 × 27876 obs: 'clusters', 'palantir_pseudotime', 'palantir_diff_potential' var: 'palantir' uns: 'clusters_colors', 'palantir_branch_probs_cell_types' obsm: 'MAGIC_imputed_data', 'X_tsne', 'palantir_branch_probs' layers: 'spliced', 'unspliced'
为了构建伪时间,数据需要进行预处理。首先,我们过滤只在少数细胞中表达的基因(这里,至少要在 20 个细胞中表达)。值得注意的是,稍后伪时间的构建对阈值的具体选择并不敏感。
在进行第一次基因过滤之后,细胞大小被归一化,计数经 log1p 转换以减少离群值的影响。和常规操作一样,我们还会识别并注释高变异基因。最后,我们基于一个最近邻图来定义伪时间。主成分的数量则根据解释方差来选择。
sc.pp.filter_genes(adata, min_counts=20)
sc.pp.normalize_total(adata)
sc.pp.log1p(adata)
sc.pp.highly_variable_genes(adata)
sc.tl.pca(adata)
sc.pp.neighbors(adata, n_pcs=10)
二维的 t-SNE 图像以细胞类型注释进行染色,结果显示细胞类型良好地聚类在一起。此外,发育层次结构也清晰可见。
sc.pl.scatter(adata, basis="tsne", color="clusters")
Pseudotime construction¶
To calculate diffusion pseudotime (DPT), first, the corresponding diffusion maps need to be calculated.
sc.tl.diffmap(adata)
尽管骨髓的分化层次结构已经被人们很好地理解,我们却并不知道发育过程究竟在哪个确切的细胞上开始,尽管我们知道它起始于造血干细胞。为了识别一个可能的初始细胞,我们研究了单个的扩散成分。我们通过一个维度(在本例中是第 3 维)上极值最高的扩散成分来识别干细胞。
# Setting root cell as described above
root_ixs = adata.obsm["X_diffmap"][:, 3].argmin()
sc.pl.scatter(
adata,
basis="diffmap",
color=["clusters"],
components=[2, 3],
)
adata.uns["iroot"] = root_ixs
sc.tl.dpt(adata)
不同的伪时间构建方法可能会得到不同的结果。有时,一种伪时间方法能够更准确地捕获潜在的发育过程。这里,我们比较了刚计算出的扩散伪时间 (DPT) 和预先计算的 Palantir 伪时间(有关教程请参见 这里)。比较不同伪时间的一种方法是对数据的低维嵌入(例如 t-SNE)进行染色。
在这里,与所有其他细胞类型相比,DPT 在 CLPs 簇中表现得非常高。相比之下,Palantir 伪时间随着发育成熟度而连续增加。这表明这两种伪时间可能捕获了发育过程的不同方面。
sc.pl.scatter(
adata,
basis="tsne",
color=["dpt_pseudotime", "palantir_pseudotime"],
color_map="gnuplot2",
)
为了避免在低维嵌入数据上进行染色,我们可以研究分配给每个细胞类型簇的伪时间值的分布。这种方法再次表明,对于 DPT 而言,CLP 簇是一个 outlier。此外,诸如 HSC_1 和 HSC_2 等簇包含一些具有较高伪时间的细胞。这些膨胀的伪时间值与我们先前对这些簇位于发育过程起始阶段的生物学认知相矛盾。
sc.pl.violin(
adata,
keys=["dpt_pseudotime", "palantir_pseudotime"],
groupby="clusters",
rotation=45,
order=[
"HSC_1",
"HSC_2",
"Precursors",
"Ery_1",
"Ery_2",
"Mono_1",
"Mono_2",
"CLP",
"DCs",
"Mega",
],
)
综合这些观察以及先前对骨髓发育的认知,我们决定继续使用 Palantir 伪时间进行后续分析。
Palantir 伪时间与我们对骨髓分化层次结构的理解(即 HSCs 位于分化层次的最开始)更加一致。DPT 分配给早期细胞类型簇 (HSC_1 和 HSC_2) 较高伪时间值,这与我们的预期不符。此外,DPT 突出了类淋巴祖细胞 (CLP) 簇,这在 Palantir 伪时间中没有观察到。因此,Palantir 伪时间似乎更好地捕获了骨髓细胞分化的方向性。
Key takeaways¶
伪时间推断 (TI) 方法需要预先知道(大致上)生物过程的起始点。
生物过程的性质决定了可以使用哪些 TI 算法。dynguidelines 可以帮助选择合适的 TI 方法。