单细胞测序最好的教程（十七）：使用StaVIA进行轨迹推断
作者按
本章节详细讲解了基于随机二阶游走的拟时序模型StaVIA在OmicVerse中的扩展应用。该算法可较为精确，快速地计算拟时序值，同时还有很多富有表现力的图像。本教程首发于单细胞最好的中文教程，未经授权许可，禁止转载。
全文字数|预计阅读时间: ～4000|10min
——Starlitnightly（星夜）
VIA是一种单细胞轨迹推断方法，提供拓扑结构构建、伪时间、自动终端状态预测以及沿谱系的时间基因动态的自动绘图功能。在这里，我们改进了原作者的着色逻辑和用户习惯，使用户可以直接使用anndata对象进行分析。
我们使用原VIA作者提供的分析完成了本教程。
论文：使用VIA在单细胞组学数据中进行广义和可扩展的轨迹推断
代码：https://github.com/ShobiStassen/VIA
Colab可重现性：https://colab.research.google.com/drive/1A2X23z_RLJaYLbXaiCbZa-fjNbuGACrD?usp=sharing
import scanpy as sc
import omicverse as ov
from omicverse.externel import VIA
import matplotlib.pyplot as plt
ov.plot_set()
____ _ _ __
/ __ \____ ___ (_)___| | / /__ _____________
/ / / / __ `__ \/ / ___/ | / / _ \/ ___/ ___/ _ \
/ /_/ / / / / / / / /__ | |/ / __/ / (__ ) __/
\____/_/ /_/ /_/_/\___/ |___/\___/_/ /____/\___/
Version: 1.6.9, Tutorials: https://omicverse.readthedocs.io/
数据预处理
作为示例，我们将差异动力学分析应用于海马齿状回神经发生，该过程包含多个异质性亚群。
import scvelo as scv
adata=scv.datasets.dentategyrus()
adata
在这里，我们加载了一个包含2930个观测值和13913个变量的AnnData对象。接下来，我们对数据进行预处理。
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)
这段代码对数据进行了以下操作：
- 识别并选择高变异基因（HVGs）。
- 对数据进行大小标准化和对数转换。
- 进行主成分分析（PCA）以降维。
接下来，我们构建邻居图并计算UMAP嵌入。
ov.pp.neighbors(adata,use_rep='scaled|original|X_pca',n_neighbors=15,n_pcs=30)
ov.pp.umap(adata,min_dist=1)
这段代码通过计算邻居来构建数据的图结构，并将其投影到UMAP空间中进行可视化。
ov.pl.embedding(adata,basis='X_umap',
color=['clusters'],
frameon='small',cmap='Reds')
这一步通过UMAP可视化了不同簇的分布。
模型构建与运行
我们需要指定用于VIA推断的细胞特征向量use_rep，可以使用X_pca, X_scVI或X_glue，取决于分析目的。在这里，我们直接使用X_pca。还需要指定使用多少个组件（不应超过向量的长度）。
此外，我们需要指定用于着色和计算的clusters。如果root_user为None，则会自动计算根细胞。
我们还需要设置存储在adata.obsm中的basis参数。一个示例是设置为tsne，因为它存储在obsm: 'tsne', 'MAGIC_imputed_data', 'palantir_branch_probs', 'X_pca'中。
其他参数和属性的解释可以在https://pyvia.readthedocs.io/en/latest/notebooks/ViaJupyter_scRNA_Hematopoiesis.html找到。
StaVia用于时间序列
StaVia用于空间-时间分析
ncomps=30
knn=15
v0_random_seed=4
root_user = ['nIPC'] #属于nIPC细胞类型的一个细胞的索引
memory = 10
dataset = ''
use_rep = 'scaled|original|X_pca'
clusters = 'clusters'
basis='X_umap'
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()
可视化和分析
在后续分析之前，我们需要指定每个簇的颜色。在这里，我们使用sc.pl.embedding自动为每个簇着色，如果需要指定自己的颜色，可以指定调色板参数。
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)
adata.obs['pt_via']=v0.single_cell_pt_markov
ov.pl.embedding(adata,basis='X_umap',
color=['pt_via'],
frameon='small',cmap='Reds')
轨迹投影
以不同方式可视化投影到二维嵌入（如UMAP, Phate, TSNE等）上的VIA轨迹。
- 在嵌入上绘制高级簇图抽象；
- 绘制高边分辨率有向图；
- 绘制细粒度细胞方向性的向量场/流图。
关键参数：
- scatter_size
- scatter_alpha
- linewidth
- draw_all_curves（如果过于拥挤，设置为False）
clusters='clusters'
color_true_list=adata.uns['{}_colors'.format(clusters)]
fig, ax, ax1 = VIA.core.plot_trajectory_curves_ov(adata,clusters='clusters',dpi=80,
via_object=v0,embedding=adata.obsm['X_umap'],
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=80,
extra_title_text='Atlas View colored by pseudotime')
fig.set_size_inches(4,4)
# edge plots can be made with different edge resolutions. Run hammerbundle_milestone_dict() to recompute the edges for plotting. Then provide the new hammerbundle as a parameter to plot_edge_bundle()
# it is better to compute the edges and save them to the via_object. this gives more control to the merging of edges
decay = 0.6 #increasing decay increasing merging of edges
i_bw = 0.02 #increasing bw increases merging of edges
global_visual_pruning = 0.5 #higher number retains more edges
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)
# via_streamplot() requires you to either provide an ndarray as embedding as an input parameter OR for via to have an embedding attribute
fig, ax = VIA.core.via_streamplot_ov(adata,'clusters',
v0, embedding=adata.obsm['X_umap'], dpi=80,
density_grid=0.8, scatter_size=30,
scatter_alpha=0.3, linewidth=0.5)
fig.set_size_inches(5,5)
#Colored by pseudotime
fig, ax = VIA.core.via_streamplot_ov(adata,'clusters',
v0,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=80,)
fig.set_size_inches(5,5)
概率路径和记忆
可视化从根到终端状态的概率路径，这由谱系可能性指示。谱系可能性越高，特定细胞向终端状态分化的潜力就越大。改变记忆参数将改变谱系路径的特异性。这可以在单细胞水平可视化，也可以与可视化细胞-细胞连接性和路径的Atlas View结合。
关键参数：
- marker_lineages（列表）终端簇
fig, axs= VIA.core.plot_sc_lineage_probability(via_object=v0, dpi=80,
embedding=adata.obsm['X_umap'])
fig.set_size_inches(12,12)
fig, axs= VIA.core.plot_atlas_view(via_object=v0, dpi=80,
lineage_pathway=[5,25,4],
fontsize_title = 12,
fontsize_labels = 12,
)
fig.set_size_inches(12,4)
可视化基因/特征图
查看基因表达在VIA图上的变化。我们使用VIA中计算的HNSW小世界图来加速基因插补计算（使用类似于MAGIC的方法）：
import pandas as pd
gene_list_magic =['Tmsb10', 'Hn1', ]
df = adata.to_df()
df_magic = v0.do_impute(df, magic_steps=3, gene_list=gene_list_magic) #optional
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)
基因动态
VIA自动推断所有检测到的谱系的基因动态。这些可以解释为沿任何给定谱系的基因表达变化。
关键参数
- n_splines
- spline_order
- gene_exp（数据框）选择基因的单细胞水平基因表达（基因插补是可选的预步骤）
fig, axs=VIA.core.get_gene_expression(via_object=v0,
gene_exp=df_magic[gene_list_magic])
fig.set_size_inches(14,4)
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)
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='result/animation_test.gif')
from IPython.display import Image
with open('result/animation_test.gif','rb') as file:
display(Image(file.read(),width=400,height=400))
VIA.core.animate_atlas(via_object=v0,
extra_title_text='test animation',
n_milestones=None,
saveto='result/edgebundle_test.gif')
from IPython.display import Image
with open('result/edgebundle_test.gif','rb') as file:
display(Image(file.read(),width=500,height=500))
adata.obs['pt_via']=v0.single_cell_pt_markov
adata.write('data/stavia_res.h5ad',compression='gzip')
以上是关于使用StaVIA进行轨迹推断的完整教程，涵盖了数据预处理、模型构建、轨迹可视化、概率路径分析以及基因动态分析等多个方面。
本文来自博客园，作者：Starlitnightly，转载请注明原文链接：https://www.cnblogs.com/starlitnightly/articles/18614638
浙公网安备 33010602011771号