从velocyto到scvelo，拟时序分析

 

个人背景：首先我使用的工具是scVelo， 我使用的编程背景是python, jupyter, scanpy。

在进行RNA velocity分析之前，我们还需要获得一个loom文件，因为不同于单细胞分析平台scanpy和seurat等需要的是基因表达矩阵。时序分析需要额外一个剪切信息，对于测序得到的reads计算其成熟度以此来估计细胞分化方向。

1、velocyto生成loom文件

安装就不在这里说了，官网很详细。运行需要准备文件：

1.1 bam文件在cellranger的文件夹下有一个outs的文件夹，其下有个名为possorted_genome_bam.bam的文件，在双组学的输出命名可能叫gex_possorted_bam.bam

1.2 cellranger使用过的基因组注释文件，genes.gtf，我使用的是refdata-gex-GRCh38-2024-A中的注释

1.3 下载表达重复注释，hg38_rmsk.gtf，从UCSC genome browser下载，Group为All Tracks，其他默认，注意你分析什么物种就下什么物种，我这个是人的。注意空格和换行，我的示例有换行方便好看而已。

velocyto run 
-b /DataPath/Cellranger/sample/outs/filtered_feature_bc_matrix/barcodes.tsv.gz 
-o /OutPath/Cellranger/sample/ 
-m hg38_rmsk.gtf /DataPath/Cellranger/sample/outs/yourfile.bam genes.gtf

这行代码就是在终端中运行的，如果你是10x单转录组学可以直接使用run10x命令更简单，我这个是更全面的命令。

velocyto run10x 
-m hg38_rmsk.gtf /DataPath/Cellranger/sample/outs/yourfile.bam genes.gtf

大约几十分钟到几小时就能输出完了，在你的输出目录下能看到一个sample.loom文件。

2、输入scanpy与anndata数据结合

各种安装包和依赖就不在这里说了，去官网下载import scvelo as scv 

可以事先准备好一个已经做完降维和注释的文件

adata = scv.read('/datapath/yourfile.h5ad', cache=True)

输入loom文件

ldata = scv.read(/datapath/.loom, cache=True)
ldata.obs.index = ldata.obs.index.str.修改为和adata中的一致
ldata.var_names_make_unique() 

如果是多个loom文件，使用

ldatas = [ldata1,ldata2,......]
ldata_combined = ldatas[0].concatenate(ldatas[1:])

一般来说，你所分析的adata都是筛选过的细胞，而loom文件下是从cellranger得到的全部细胞，所以需要根据adata的索引选取对应的细胞

cellselect = adata.obs.index.tolist()
cellselect

filtered_ldata = ldata_combined[ldata_combined.obs.index.isin(cellselect)]
filtered_ldata

这样细胞数就一致且对齐了，然后整合在一起

adata = scv.utils.merge(adata, filtered_ldata)

3、计算时序，可视化

可视化你的注释

sc.pl.umap(adata, color='cell_type', frameon=False, legend_loc='on data')

查看细胞类型中的剪切与未剪切情况

scv.pl.proportions(adata, groupby='cell_type')

处理数据

scv.pp.filter_and_normalize(adata, min_shared_counts=20, n_top_genes=2000)
scv.pp.moments(adata, n_pcs=50, n_neighbors=30)

 如果这一步报错write_knn_indices，就在第一步输入adata之后紧接着运行

sc.pp.pca(adata)
sc.pp.neighbors(adata, n_pcs=50, n_neighbors=30)

计算时序

# compute velocity 
scv.tl.velocity(adata, mode='stochastic') 
scv.tl.velocity_graph(adata) 

大致就能看到不错的结果了，还可以优化一下可视化

scv.pl.velocity_embedding_stream(adata, basis='umap',color="cell_type", min_mass=0, cutoff_perc= 0, palette = palette, legend_loc="right margin", size=3, alpha=1, arrow_size=0.5, linewidth=0.8, dpi=400)

至此在你的umap就有期刊中那种好看的时序方向了，还想要更多的可视化就去官网进一步学习吧！