4.基于python的scRNA-seq轨迹推断-伪时间计算

参考：
[1] https://github.com/Starlitnightly/single_cell_tutorial
[2] https://github.com/theislab/single-cell-best-practices

拟时序计算分类

拟时序（伪时间）根据细胞在发展过程中的阶段对细胞进行排序。不成熟的细胞被分配小值，成熟的细胞被分配大值。例如，研究骨髓样本时，造血干细胞被分配低拟时序值，而红细胞则被分配高拟时序值。并且，在后续分析拟时序时，通常需要指定一个初始细胞，或者叫根细胞，然后从起始点推测细胞的动态过程。

拟时序计算遵循以下步骤，首先将高维细胞数据投影到低维表示，然后根据生物过程本身的性质选择构造拟时序的方法，生物过程本身的性质需要我们提前了解，比如是线性的，循环的，分叉的等。

拟时序方法分为下面几类：

• 聚类方法：基于树的方法----首先对细胞表示进行聚类，然后识别簇之间的连接。这些簇可以被排序，从而构建伪时间。经典的聚类算法包括K-Means，Leiden，或者层次聚类。然后用最小生成树MST来连接这些簇，得到排序
• 基于流形学习的方法：用主曲线或图来估计潜在的轨迹，代表方法是Slingshot
• 基于状态转移概率的方法：每个转移概率量化了参考细胞是其他细胞祖先的可能性。这些概率定义了用于定义伪时间的随机过程。例如，DPT定义为随机游走的连续状态之间的差异。这种方法需要指定一个根细胞，然后伪时间本身是相对于这个细胞计算的

得到伪时间后，可以用PAGA和跟细胞信息对带有时序的拓扑结构进行抽象，得到轨迹。或者直接用PAGA进行聚类得到图摘要，此时是一个没有方向的无向图。

加载数据

我们使用小鼠齿状回神经元的单细胞RNA-seq测序数据，进而研究nIPC的分化过程。我们使用omicverse预处理数据，选择前3000个基因作为高可变基因。另外，我们应该考虑多少个主成分来计算细胞之间邻域关系的信息。根据经验，通常粗略估计主成分的数量就可以满足需求。一般而言前50个主成分即涵盖大部分的数据变异性。：

import omicverse as ov
import scanpy as sc
import scvelo as scv

adata = sc.read("./data/DentateGyrus/10X43_1.h5ad")
print(adata)
"""
AnnData object with n_obs × n_vars = 2930 × 13913
    obs: 'clusters', 'age(days)', 'clusters_enlarged'
    uns: 'clusters_colors'
    obsm: 'X_umap'
    layers: 'ambiguous', 'spliced', 'unspliced'
"""

adata=ov.pp.preprocess(adata,mode='shiftlog|pearson',n_HVGs=3000,)
adata.raw = adata
adata = adata[:, adata.var.highly_variable_features]
ov.pp.scale(adata)
ov.pp.pca(adata,layer='scaled',n_pcs=50)

DPT扩散+PAGA

PAGA可以估计图摘要，PAGA可以保留数据的全局拓扑，并以不同分辨率去分析，PAGA需要先定义发育的跟细胞，然后用扩散方法计算伪时间：

Traj=ov.single.TrajInfer(adata,basis='X_umap',use_rep='scaled|original|X_pca',n_comps=50)
Traj.set_origin_cells('nIPC')
Traj.inference(method='diffusion_map')

我们可以可视化每个细胞的拟时序值，进而观察细胞的发育过程：

import numpy as np
adata.obs['dpt_pseudotime']=adata.obs['dpt_pseudotime'].fillna(0)
adata.obs['dpt_pseudotime'].replace([np.inf], 1, inplace=True)
adata.obs['dpt_pseudotime'].replace([-np.inf], 0, inplace=True)

ov.pl.embedding(adata,basis='X_umap',
                   color=['clusters','dpt_pseudotime'],
                   frameon='small',cmap='Reds', wspace=0.45)

我们还可以研究分配给每个细胞簇的伪时间值的分布，而不是给数据的低维表示着色：

sc.pl.violin(
    adata,
    keys=["dpt_pseudotime"],
    groupby="clusters",
    rotation=90,
    order=[
        "nIPC",
        "Neuroblast",
        "Granule immature",
        "Granule mature",
        "Mossy",
        "Radial Glia-like",
        "Astrocytes",
        "Microglia",
        "OPC",
        "OL",
        'Cajal Retzius', 'Cck-Tox',"GABA","Endothelial"
    ],
)

我们还可以计算每一类细胞类型转换的状态转移概率矩阵，PAGA 图的摘要被标准化为轨迹推理的最佳方法。它提供了一个数据拓扑的图形化映射，其加权边对应于两个簇之间的连通性：

ov.utils.cal_paga(adata,use_time_prior='dpt_pseudotime',vkey='paga',
                 groups='clusters')
ov.utils.plot_paga(adata,basis='umap', size=50, alpha=.1,title='PAGA LTNN-graph',
            min_edge_width=2, node_size_scale=1.5,show=False,legend_loc=False)

从上面结果看出，分化轨迹似乎是向两个方向扩展。

其余方法的使用标准和上面的类似：比如上面的是DPT+PAGA，可以换成流形方法+PAGA，总之要计算出伪时间，然后用PAGA获得图摘要用于模拟轨迹。

补充：仅使用PAGA

仅使用PAGA获得图摘要：

adata = sc.read("./data/DentateGyrus/10X43_1.h5ad")
print(adata)

adata=ov.pp.preprocess(adata,mode='shiftlog|pearson',n_HVGs=3000,)
adata.raw = adata
adata = adata[:, adata.var.highly_variable_features]
ov.pp.scale(adata)
ov.pp.pca(adata,layer='scaled',n_pcs=50)

sc.pp.neighbors(adata, n_neighbors=4, n_pcs=30)
sc.tl.paga(adata, groups='clusters')
sc.pl.paga(adata, color=['clusters'], frameon=False, fontsize=6, edge_width_scale=0.25, threshold=0.1, node_size_scale=3)

我们发现如果没有伪时间，PAGA的图摘要能反映一定的拓扑连接，但是不够清楚。