scRNA——StaVIA细胞轨迹

最新推荐文章于 2025-05-20 18:44:14 发布

以待成追忆

最新推荐文章于 2025-05-20 18:44:14 发布

阅读量726

点赞数 13

分类专栏： scRNA 文章标签： python 数据分析

本文链接：https://blog.youkuaiyun.com/qq_24680545/article/details/145380140

版权

scRNA 专栏收录该内容

8 篇文章

订阅专栏

StaVIA

VIA是一种单细胞轨迹推理方法，可提供拓扑结构，伪时间，自动终末状态预测和沿谱系的时间基因动力学的自动绘制。

数据加载和预处理

import matplotlib.pyplot as plt
import omicverse as ov
import scvelo as scv
from omicverse.externel import VIA

这里使用测试用的训练集，scv.datasets.dentategyrus()这一训练集是齿状回神经发生(dentate gyrus neurogenesis)。

def prep():
    adata = scv.datasets.dentategyrus()
    adata = ov.pp.preprocess(adata, mode='shiftlog|pearson', n_HVGs=2000, )
    adata.raw = adata
    adata = adata[:, adata.var.highly_variable_features]
    ov.pp.scale(adata)
    ov.pp.pca(adata, layer='scaled', n_pcs=50)
    ov.pp.neighbors(adata, use_rep='scaled|original|X_pca', n_neighbors=15, n_pcs=30)
    ov.pp.umap(adata, min_dist=1)
    traj_embed = ov.pl.embedding(adata, basis='X_umap',
                    color=['clusters'],
                    frameon='small', cmap='Reds',show=False, return_fig=True)
    traj_embed.figure.savefig('dentategyrus_clusters.png', dpi=300, bbox_inches='tight')
    adata.write_h5ad('../data/dentategyrus.h5ad')

在linux和windows系统的不同，在windwos系统上需要将其放在if __name__ == '__main__''中才能在执行多进程任务的时候不会出现资源冲突，所以需要将其放在一个函数下，并将函数放在if __name__ == '__main__'里面执行

if __name__ == '__main__':
    prep()

在这里插入图片描述

构建模型

需要指定用于VIA推断的细胞特征向量use_rep，可以是X_pca, X_scVI或者X_glue，这次使用的是X_pca，也就是以主成分分析的结果作为特征向量。同时还要指定一些参数：

components，这个参数不能超过特征向量的长度，例如，如果X_pca有50个主成分，用户不能指定使用100个，过多可能导致过拟合，过少可能丢失信息。
clusters，这个参数用于VIA计算，要求在分析前，数据已经被注释细胞类型了。
root_user，这个参数是指定根细胞，而我们使用的是神经细胞，因此指定的是nIPC。

def via():
    adata = ov.read('../data/dentategyrus.h5ad')
    ncomps = 30
    knn = 15
    v0_random_seed = 4
    root_user = ['nIPC']
    memory = 10
    dataset = ''
    use_rep = 'scaled|original|X_pca'
    clusters = 'clusters'
	v0 = VIA.core.VIA(data=adata.obsm[use_rep][:, 0:ncomps],
                      true_label=adata.obs[clusters],
                      edgepruning_clustering_resolution=0.15, cluster_graph_pruning=0.15,
                      knn=knn, root_user=root_user, resolution_parameter=1.5,
                      dataset=dataset, random_seed=v0_random_seed, memory=memory)

    v0.run_VIA()
    # 查看
    fig, ax, ax1 = VIA.core.plot_piechart_viagraph_ov(adata, clusters='clusters', dpi=80,
                                                      via_object=v0, ax_text=False, show_legend=False)
    fig.set_size_inches(8, 4)
    fig.savefig('dentategyrus_via.png', dpi=300, bbox_inches='tight')

在后续分析之前，我们需要指定每个簇的颜色。在这里，我们使用sc.pl.embedding自动为每个簇着色，如果需要指定自己的颜色，可以指定调色板参数。

	fig, ax, ax1 = VIA.core.plot_piechart_viagraph_ov(adata, clusters='clusters', dpi=300,
                                                      via_object=v0, ax_text=False, show_legend=False)
    fig.set_size_inches(8, 4)
    fig.savefig('dentategyrus_via.png', dpi=300, bbox_inches='tight')
    adata.obs['pt_via'] = v0.single_cell_pt_markov
    pt_via = ov.pl.embedding(adata, basis='X_umap',
                    color=['pt_via'],
                    frameon='small', cmap='Reds', show=False, return_fig=True)
    pt_via.savefig('dentategyrus_via_embedding.png', dpi=300, bbox_inches='tight')

dentategyrus_via.png
在这里插入图片描述
 dentategyrus_via_embedding.png

轨迹投影

以不同的方式可视化投影到2D嵌入（UMAP，PHATE，TSNE等）上的整体轨迹。

	fig, ax, ax1 = VIA.core.plot_trajectory_curves_ov(adata, clusters='clusters', dpi=300,
                                                      via_object=v0, embedding=adata.obsm['X_umap'],
                                                      draw_all_curves=False)
    fig.set_size_inches(16, 8)
    fig.savefig('dentategyrus_via_curves.png', dpi=300, bbox_inches='tight')

在这里插入图片描述
当上面的文字过于拥挤时，可以将draw_all_curves设置为False。

阿特拉斯分析图

	v0.embedding = adata.obsm['X_umap']
    fig, ax = VIA.core.plot_atlas_view(via_object=v0,
                                       n_milestones=150,
                                       sc_labels=adata.obs['clusters'],
                                       fontsize_title=12,
                                       fontsize_labels=12, dpi=300,
                                       extra_title_text='Atlas View colored by pseudotime')
    fig.set_size_inches(6, 6)
    fig.savefig('dentategyrus_via_atlas.png', dpi=300, bbox_inches='tight')

在这里插入图片描述
直观地看到细胞随着伪时间的推移而发生变化。

绘制边缘图

    decay = 0.6  # 增加衰减的边缘融合
    i_bw = 0.02  # 增加BW增加了边缘的合并
    global_visual_pruning = 0.5  # 较高的数字保留更多的边缘
    n_milestones = 200

    v0.hammerbundle_milestone_dict = VIA.core.make_edgebundle_milestone(
    												via_object=v0,
                                                    n_milestones=n_milestones,
                                                    decay=decay, 
                                                    initial_bandwidth=i_bw,
													global_visual_pruning=global_visual_pruning)
	fig, ax = VIA.core.plot_atlas_view(via_object=v0,
                                       add_sc_embedding=True,
                                       sc_labels_expression=adata.obs['clusters'],
                                       cmap='jet', sc_labels=adata.obs['clusters'],
                                       text_labels=True,
                                       extra_title_text='Atlas View by Cell type',
                                       fontsize_labels=3, fontsize_title=3, dpi=300
                                       )
    fig.set_size_inches(6, 4)
    fig.savefig('dentategyrus_via_edgebundle.png', dpi=300, bbox_inches='tight')

在这里插入图片描述
可以看到，非常的的杂乱，所以我们可以改一下代码。

    fig, ax = VIA.core.via_streamplot_ov(adata, 'clusters',
                                         v0, embedding=adata.obsm['X_umap'], dpi=300,
                                         density_grid=0.8, scatter_size=30,
                                         scatter_alpha=0.3, linewidth=0.5)
    fig.set_size_inches(5, 5)
    fig.savefig('dentategyrus_via_streamplot2.png', dpi=300, bbox_inches='tight')

在这里插入图片描述
然后我们可以绘制一下根据伪时间轨迹的边缘图。

    fig, ax = VIA.core.via_streamplot_ov(adata, 'clusters',
                                         v0,embedding=adata.obsm['X_umap'], 
                                         density_grid=0.8, 
                                         scatter_size=30, 
                                         color_scheme='time', 
                                         linewidth=0.5,
                                         min_mass=1, 
                                         cutoff_perc=5, 
                                         scatter_alpha=0.3, 
                                         marker_edgewidth=0.1,
                                         density_stream=2, 
                                         smooth_transition=1, 
                                         smooth_grid=0.5, 
                                         dpi=300, )
    fig.set_size_inches(5, 5)
    fig.savefig('dentategyrus_via_streamplot1.png', dpi=300, bbox_inches='tight')

在这里插入图片描述
从这个时间图上可以看到，趋势与伪时间(pseudotime)相同。

细胞轨迹的概率建模（Probabilistic pathways and Memory）

可视化从根状态到终末状态的概率路径，路径的谱系可能性越高，表示该细胞向目标终末状态分化的潜力越大。调整记忆参数会改变谱系路径的特异性。此分析可以在单细胞层面进行可视化，也可以与 Atlas View 结合，展示细胞间的连接性和路径。

	fig, axs = VIA.core.plot_sc_lineage_probability(via_object=v0, dpi=300,
                                                    embedding=adata.obsm['X_umap'])
    fig.set_size_inches(12, 12)
    fig.savefig('dentategyrus_via_lineage_probability.png', dpi=300, bbox_inches='tight')

在这里插入图片描述
当然也可以指定clusters的编号，也要在clusters里面.

    fig, axs = VIA.core.plot_atlas_view(via_object=v0, dpi=300,
                                        lineage_pathway=[5, 25, 4],# 这几个在clusters里面
                                        fontsize_title=12,
                                        fontsize_labels=12,
                                        )
    fig.set_size_inches(12, 4)
    fig.savefig('dentategyrus_via_lineage_pathway.png', dpi=300, bbox_inches='tight')

在这里插入图片描述
可以看到，各个cluster都有比较明显的概率分布位置。

可视化基因图

假设想要看Tmsb10和Hn1的基因表达情况，可以查看沿着 VIA 图谱查看基因表达。我们使用 VIA 中计算的 HNSW 小世界图来加速基因补充计算（采用类似 MAGIC 的方法）。

    gene_list_magic = ['Tmsb10', 'Hn1', ]
    df = adata.to_df()
    df_magic = v0.do_impute(df, magic_steps=3, gene_list=gene_list_magic)
    # 查看基因表达在VIA图上的变化。我们使用VIA中计算的HNSW小世界图来加速基因插补计算
    fig, axs = VIA.core.plot_viagraph(via_object=v0,
                                      type_data='gene',
                                      df_genes=df_magic,
                                      gene_list=gene_list_magic[0:3], arrow_head=0.1)
    fig.set_size_inches(12, 4)
    fig.savefig('dentategyrus_via_gene.png', dpi=300, bbox_inches='tight')

在这里插入图片描述

基因动力图/热力图

VIA 会自动推断所有检测到的谱系沿伪时间的基因动态。这些变化可以解释为沿任意谱系的基因表达变化。

    # 基因动力图
    fig, axs = VIA.core.get_gene_expression(via_object=v0,  # cmap_dict=color_dict,
                                            gene_exp=df_magic[gene_list_magic])
    fig.set_size_inches(14, 4)
    fig.savefig('dentategyrus_via_gene_expression.png', dpi=300, bbox_inches='tight')

    # 可视化基因表达趋势的热图
    fig, axs = VIA.core.plot_gene_trend_heatmaps(via_object=v0,
                                                 df_gene_exp=df_magic, cmap='plasma',
                                                 marker_lineages=[34, 5])
    fig.set_size_inches(5, 5)
    fig.savefig('dentategyrus_via_gene_trend_heatmaps.png', dpi=300, bbox_inches='tight')

在这里插入图片描述

这个图的图例怎么放外面还没找到，先将就看吧。

在这里插入图片描述

轨迹图动态演示

via提供懂了一个函数animate_streamplot_ov可以查看轨迹图的动态过程，说白了就是gif文件。

VIA.core.animate_streamplot_ov(adata, 'clusters', v0, embedding=adata.obsm['X_umap'],
                                   cmap_stream='Blues',
                                   scatter_size=200, scatter_alpha=0.2, marker_edgewidth=0.15,
                                   density_grid=0.7, linewidth=0.1,
                                   segment_length=1.5,
                                   saveto=r'temp\animation_test.gif')