官网的代码演示差不多已经能够完成我的部分需求，但是有部分函数他介绍的不太完善，我进行了一些修改，并对脚本进行了封装。
step1 数据读取
1.读取loom文件，将所有的loom文件都放在一个名为loom_file的文件夹里，并进行读取。
loom_data = {}
index_sample = []
# 遍历 loom 文件列表
for file in loom_file:
# 检查文件是否为 loom 文件
if '.loom' in file:
# 打印文件路径
print(f"{loom_dir}/{file}")
# 读取 loom 文件数据
loom_temp_data = sc.read(f"{loom_dir}/{file}", cache=True)
# 确保基因名称唯一
loom_temp_data.var_names_make_unique()
# 替换观测名称中的特殊字符,这段不固定，需要根据你的loom里面的barcode和单细胞数据里面的barcode进行统一
loom_temp_data.obs_names = loom_temp_data.obs_names.str.replace("x", "")
# 将样本索引添加到列表中
index_sample += list(loom_temp_data.obs_names)
# 将数据存储在 loom 数据字典中
loom_data[file] = loom_temp_data
2.读取rds或者h5ad文件，rds文件需要提前转换成h5ad文件
sc_data = scv.read('test.h5ad')
loom_data_merge = loom_data[list(loom_data.keys())[0]]
# 遍历 loom_data 中的其他 loom 数据
for i in range(1, len(loom_data)):
# 获取下一个 loom 数据
next_data = loom_data[list(loom_data.keys())[i]]
# 将当前 loom 数据与下一个 loom 数据连接
loom_data_merge = loom_data_merge.concatenate(next_data)
loom_data_merge.obs_names=index_sample
adata = scv.utils.merge(sc_data, loom_data_merge)
3.读取一下每个cluster的剪切和未剪切的比例,也可以观察其他的信息，比如sample，group，celltype等
groupby='clusters'
output='scvelocity_outdir'
scv.pl.proportions(adata, groupby=groupby,
save=f'{output}/spliced_proportions_groupby_{groupby}.pdf')
step2 数据处理
1.数据预处理，有些人可能之前已经用scanpy做过一些数据的预处理，会疑问这一步是否还有必要进行，作者给出的答案是就算之前已经做过了预处理，但是也没有对spliced和unspliced进行处理，还是要运行这一步，该函数会自动判断数据是否进行过标准化，如今已经进行了，则不会再log处理。
scv.pp.filter_and_normalize(adata，min_shared_counts=20, n_top_genes=2000)
scv.pp.moments(adata, n_pcs=30, n_neighbors=30)
step3 速度计算
scvelo支持动态模型和稳态模型，使用动态模型比稳态模型能够获得更多的信息。
通过代表基因推算细胞的未来分化方向，并投影到低纬度上。
scv.tl.recover_dynamics(adata)
scv.tl.velocity(adata,vkey='dynamical_velocity', mode='dynamical',n_jobs=5)
scv.tl.velocity_graph(adata,vkey='dynamical_velocity',n_jobs=5)
step4 结果可视化
作者提供三种速度投影到umap的展示方式
细胞水平 scv.pl.velocity_embedding，
网格线 scv.pl.velocity_embedding_grid
流线型 scv.pl.velocity_embedding_stream
#第一张
scv.pl.velocity_embedding_stream(adata, basis='umap',key='dynamical_velocity',
title=f'velocity_embedding_stream_groupby_{groupby}',color=groupby,
legend_fontsize=18, legend_loc='right margin',colorbar=TRUE,dpi=300,
save=f'{output}/velocity_embedding_stream_groupby{groupby}.svg')
image.png
#第二张
scv.pl.velocity_embedding(adata, basis='umap',vkey='dynamical_velocity',
title=f'velocity_embedding_groupby_{groupby}',
arrow_length=2,color=groupby,arrow_size=1, legend_loc='right margin',
dpi=300,save=f'{output}/velocity_embedding_groupby_{groupby}.svg')
image.png
#第三张
scv.pl.velocity_embedding_grid(adata, basis='umap',vkey='dynamical_velocity',
color=groupby,dpi=300,title=f'velocity_embedding_grid_groupby_{groupby}',
legend_loc='right margin',save=f'{output}/velocity_embedding_grid_groupby_{groupby}.pdf')
image.png
step5 latent time
#计算潜伏时间
scv.tl.recover_latent_time(adata, vkey='dynamical_velocity')
scv.pl.scatter(adata,basis='umap',
color='latent_time',
color_map='gnuplot',
colorbar=True,
perc=[2, 98], #指定连续着色的百分比。
rescale_color=[0,1], #颜色重新缩放的边界，设置颜色条的最小/最大值。
dpi = 300,
save=f'{output}/latent_time_based_dynamical.pdf')
image.png
#取出Top300的likelihood genes信息，并保存，绘制热图
TopGenes = adata.var['fit_likelihood'].sort_values(ascending=False).index[:300]
df = pd.DataFrame(TopGenes)
df=df.transpose()
df.to_csv(f'{output}/Top300_likelihood_genes_data.csv')
#绘制热图
scv.pl.heatmap(adata, var_names=TopGenes,
sortby='latent_time',
col_color=groupby,#选择作为注释条的列名
n_convolve=100,#平滑曲线的核数
dpi=400, #这个图片很大，不要保存太大的dpi
save = f'{output}/latent_time_heatmap_groupby_{groupby}.pdf')
image.png
step6 Top-likelihood genes
动态模型中的高似然基因可以认为是驱动基因。这些基因展现出明显的动态特征，并且可以明显的地被检测到
scv.pl.scatter(adata2,
basis=top_genes[:15],
ncols=5,
use_raw=False,
frameon=False,
legend_loc='right margin',
save=f'{output}/Top-likelihood_genes_groupby_{groupby}.pdf')
