scVelo是现今使用最为广泛的、基于RNA velocity进行细胞轨迹（trajectory）的工具之一，它基于python。而seurat则是如今使用最广泛的对scRNA-seq数据进行预处理以及分群的工具，它基于R。因此在实践过程中，我们经常会遇见的一个问题，即如何在seurat上进行预处理和分群后，在scVelo上进行细胞轨迹（trajectory）分析。
在开始分析之前，我们需要准备以下数据/软件：
1. 用cellranger对10x测序数据分析后获得的结果文件夹cellranger_output
2. 已包含分群（clustering）结果的seruat object
3. 安装velocyto（http://velocyto.org/velocyto.py/install/index.html）及scVelo （https://scvelo.readthedocs.io/installation/）
一、获得loom文件scVelo的分析基于细胞中spliced mRNA和unspliced mRNA的reads数量。因此我们要先得到包含这一信息的loom文件。这一步可以通过velocyto完成。现在最主流的scRNA测序技术是10x，因此我们以10x获得的结果举例。关于其它类型的scRNA-seq数据，可以参考：http://velocyto.org/velocyto.py/tutorial/cli.html#running-velocyto
velocyto run10x -m
1. run10x表示我们使用数据为cellranger分析后的结果文件夹，在这，我们把这个文件夹命名为cellranger_output/
2. hg38_rmsk.gtf为注释文件，包含一些重复表达序列的信息，有助于获得更准确的结果。可以从UCSC（https://genome.ucsc.edu/cgi-bin/hgTables?hgsid=611454127_NtvlaW6xBSIRYJEBI0iRDEWisITa&clade=mammal&org=Mouse&db=mm10&hgta_group=allTracks&hgta_track=rmsk&hgta_table=0&hgta_regionType=genome&position=chr12%3A56694976-56714605&hgta_outputType=primaryTable&hgta_outputType=gff&hgta_outFileName=mm10_rmsk.gtf）中获得。在挑选好我们想要的genome（红框）之后，就可以点击get output（红箭头）得到。
3. refdata-gex-GRCh38-2020-A/genes/genes.gtf为参考基因组的注释文件，我们用和cellranger分析时一样的即可。也可以从GENCODE（https://www.gencodegenes.org/）获得。如果是人的基因组，先选择human，然后下载第一个gtf文件
4. <>中为需要替换的部分，下同
最终我们会获得位于cellranger_output/velocyto/这个文件夹里的loom文件。
二、从Seurat到scVelo接下来就是将seurat上的分析结果迁移到python上用于后续scVelo的分析。scVelo用于分析的数据类型是AnnData我们可以选择一些包将seurat object变成AnnData。例如：https://mojaveazure.github.io/seurat-disk/articles/convert-anndata.html。但在实际情况中，小L发现这些工具要么只适用特定版本的seurat object和AnnData，要么无法选择想要保留在AnnData中的数据导致最后无法与loom文件中的数据匹配。因而更为稳定且方便的方法是直接将seurat object中我们需要用到的数据保存成txt或者csv文件，再添加到读取loom文件后获得的AnnData中。
我们需要从seurat中获取的信息有三样，我们可以通过以下代码在R中将它们提取出来并保存，其中seuratObj为我们的seurat object：
1. 我们经过preproceess之后保留的细胞
2. 我们计算出的UMAP1、UMAP2的值，用于画图
3. 对于细胞的分类
## cell ID write.csv(Cells(), file = "cellID_obs.csv", row.names = FALSE) ## UMAP write.csv(Embeddings( , reduction = "umap"), file = "cell_embeddings.csv") ## seurat clustering write.csv(Idents( ), "clusters.csv")
接着，我们需要在python中读取loom文件
import anndata import scvelo as scv import pandas as pd import numpy as np ldata = scv.read(<"10x.loom">, cache=True) ldata.var_names_make_unique()
修改ldata中barcode以及在上一步中通过R保存的文件里的barcode，使它们相匹配，并将UMAP的信息和seurat中导出的分群信息添加到AnnData中。如果我们用velocyto run10x来获得loom文件的话，loom文件中的细胞数量会比我们通过seurat读filtered_feature_bc_matrix得到的细胞数多，因此我们要从AnnData中根据barcode挑选出我们在seurat中保留的细胞。
## 把AnnData里barcode中的'v3:|x'去掉
obs_name_rep = ldata.obs_names.str.replace('v3:|x','')
ldata.obs_names = obs_name_rep
## 读取barcode、UMAP以及seurat的分群结果
cellID_obs = pd.read_csv("cellID_obs.csv")
umap_cord = pd.read_csv("cell_embeddings.csv")
cell_clusters = pd.read_csv("clusters.csv")
## 把seurat导出文件里barcode中的'-1'去掉
cellID_obs['x'] = cellID_obs['x'].str.replace('-1', '')
## 根据barcode选出我们要的细胞
filtered_ldata = ldata[cellID_obs['x']].copy()
## 添加UMAP信息
ldata_index = pd.DataFrame(filtered_ldata.obs.index)
umap_cord = umap_cord.rename(columns = {'Unnamed: 0':'CellID'})
umap_cord['CellID'] = umap_cord['CellID'].str.replace('-1', '')
umap_ordered = ldata_index.merge(umap_cord, on = "CellID")
umap_ordered = umap_ordered.iloc[:,1:]
filtered_ldata.obsm['X_umap'] = umap_ordered.values
## 添加分群信息
cell_clusters = cell_clusters.rename(columns = {'Unnamed: 0':'CellID'})
cell_clusters['CellID'] = cell_clusters['CellID'].str.replace('-1', '')
cell_clusters = ldata_index.merge(cell_clusters, on = "CellID")
cell_clusters = cell_clusters.iloc[:,1:]
filtered_ldata.obs['seurat_clusters'] = cell_clusters.values
接着，我们就可以开始用scVelo进行细胞的轨迹（trajectory）分析。这部分可以参考scVelo的tutorial：RNA Velocity Basic （https://scvelo.readthedocs.io/VelocityBasics/），从Preprocess the Data开始。因为有官方更详细的tutorial，在这小L就只展示一个简单的范例
scv.pp.filter_and_normalize(filtered_ldata, min_shared_counts=20, n_top_genes=2000) scv.pp.moments(filtered_ldata, n_pcs=30, n_neighbors=30) scv.tl.velocity(filtered_ldata) scv.tl.velocity_graph(filtered_ldata) scv.pl.velocity_embedding_stream(filtered_ldata, basis="umap", color="seurat_clusters")
祝大家吃好喝好睡好，科研快乐~
欢迎关注微信公众号 “小L的读博日常”，第一时间获得更多和生物信息学相关的小tips。
Ref：
http://velocyto.org/velocyto.py/tutorial/cli.html#run10x-run-on-10x-chromium-samples
