11.1 前言
单细胞测序分析提供生物组织的高分辨率测量。因此,此类技术可以帮助破译和理解细胞异质性和生物过程的动态。相应的研究包括量化细胞命运以及识别驱动该过程的基因。然而,由于在经典单细胞 RNA 测序(scRNA-seq)方案中进行测序时细胞会被破坏,因此无法跟踪它们的发育情况,例如随时间变化的基因表达谱。因此,需要根据测量的快照snapshot数据来估计底层动态过程。
尽管传统上样品是从单个实验时间点采集的,但仍可以观察到多种细胞类型。这种多样性源于生物过程的异步性质。因此,可以观察到一系列的发育过程。重建发展过程是该领域的轨迹推断(trajectory inference,TI)的目标。这项任务是通过根据发育过程对观察到的细胞状态进行排序来实现的。通过将离散注释映射到连续域(即所谓的伪时间),状态沿着发展方向对齐。
伪时间根据细胞在发育过程中各自的阶段对细胞进行相对排列。不太成熟的细胞分配较小的值,成熟的细胞分配较大的值。例如,在研究骨髓样本时,造血干细胞被指定为低伪时间,而红细胞被指定为高伪时间。对于单细胞RNA测序数据,分配是基于细胞的转录组图谱。此外,构建通常需要指定整个过程开始的初始或等效的根细胞。
11.2 拟时序的建立
伪时间构建通常遵循一个常见的工作流程:第一步,将超高维单细胞数据投影到较低维的表示上。这一过程通过观察动态过程在低维流形上进展而得到证实。在实践中,伪时间方法可能依赖于主成分或扩散成分(例如扩散伪时间 (DPT) 。接下来,伪时间是根据以下原则之一构建的。
首先对观察结果进行聚类,然后确定这些聚类之间的连接。可以对簇进行排序,从而构建伪时间。今后,我们将这种方法称为集群方法。经典的聚类算法包括k-means,Leiden或层次聚类。簇可以基于相似性或通过构建最小生成树(MST)来连接。
图方法graph approach首先找到观测值的低维表示之间的联系。此过程定义了一个图表,基于该图表定义了集群,从而定义了排序。例如,PAGA将图划分为Leiden簇并估计它们之间的连接。直观上,这种方法保留了数据的全局拓扑,同时以较低的分辨率进行分析。因此,提高了计算效率。
基于流形学习的方法Manifold-learning based approaches与集群方法类似。然而,簇之间的连接是通过使用主曲线或图形来估计基础轨迹来定义的。主曲线找到一条连接高维空间中的细胞观察结果的一维曲线。这种方法的一个著名代表是Slingshot。
概率框架将转移概率分配给有序的细胞对。每个转移概率都量化参考细胞是另一个细胞的祖先的可能性。这些概率定义了用于定义伪时间的随机过程。例如,DPT被定义为随机游走的连续状态之间的差异。相比之下,Palantir将轨迹本身建模为马尔可夫链。虽然这两种方法都依赖于概率框架,但它们需要指定根细胞。伪时间本身是相对于该细胞计算的。
TI是一个经过深入研究的领域,提供了丰富的方法。要应用适当的方法来分析单细胞数据集,首先需要了解生物过程本身。这种理解尤其包括该过程的性质,即,例如,它是否是线性的、环状的或支化的。类似地,同一数据集中的正交过程限制了适用的TI方法。为了帮助识别合适的工具,dynguidelines[Deconinck et al., 2021]提供了算法及其特征的详尽概述。
11.3 下游任务及展望
尽管TI和伪时间已经可以提供有价值的见解,但它们通常充当更细粒度分析的垫脚石。例如,识别终态是一个可以研究的经典生物学问题。同样,可以根据TI和伪时间来识别谱系分歧和驱动命运决定的基因。可以回答哪个问题以及如何找到答案通常是特定于方法的。例如,Palantir将终端状态识别为其构建的马尔可夫链的吸收状态。
轨迹推断的成功有据可查,因此,人们提出了许多方法。然而,随着测序技术的进步,新的信息来源变得可用。例如,ATAC-seq、CITE-seq和 DOGMA-seq可测量转录组以外的其他模式。谱系追踪和代谢标记甚至提供给定细胞的(可能的)未来状态。因此,未来的TI工具将能够包含更多信息来更准确、更稳健地估计轨迹和伪时间,并允许回答新问题。例如,RNA velocity是一种使用未剪接和剪接mRNA来推断超越经典静态信息的定向动态信息的技术。
11.4 推断成人骨髓的伪时间
为了展示如何构建伪时间并比较不同的伪时间,我们研究了成人骨髓的数据集[Setty et al., 2019]。
11.4.1 环境设置
from pathlib import Path
import scanpy as sc
11.4.2 常规设置
DATA_DIR = Path("../../data/")
DATA_DIR.mkdir(parents=True, exist_ok=True)
FILE_NAME = DATA_DIR / "bone_marrow.h5ad"
11.4.3 数据加载
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")
11.4.4 拟时序构建
为了计算扩散拟时间(diffusion pseudotime,DPT),首先需要计算相应的扩散图。
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)进行着色。在这里,与所有其他细胞类型相比,CLP簇中的DPT极高。相比之下,Palantir假时间随着发育成熟而不断增加。
sc.pl.scatter(
adata,
basis="tsne",
color=["dpt_pseudotime", "palantir_pseudotime"],
color_map="gnuplot2",
)
我们可以研究分配给每个细胞类型簇的伪时间值的分布,而不是对数据的低维表示进行着色。该表示再次表明CLP簇在DPT情况下形成异常值。此外,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拟时间。