scVelo于2020年发表在 Nature Biotechnology ，对最初的 RNA velocity研究 及其附带的软件 velocyto 进行了多项改进。
本教程是介绍如何使用python包scVelo对单细胞RNA-seq数据(scRNA-seq)进行RNA速率分析。

参考资料：

scvelo github网站： https://github.com/theislab/scvelo
scvelo官方文档： https://scvelo.readthedocs.io/index.html
Seurat to RNA-Velocity教程： https://github.com/basilkhuder/Seurat-to-RNA-Velocity#multiple-sample-integration
因为scvelo为python包，seurat数据在使用scvelo进行分析时，需要获取seurat对象中的meta-data数据，上面的教程可以指导如何生成scvelo输入文件。
scvelo实战教程： https://smorabit.github.io/tutorials/8_velocyto/

怎么理解RNA速率?

La Manno等人2018年在 Nature 杂志上发表《RNA velocity of single cells》，提出了RNA速率分析方法。

RNA丰度是分析单个细胞状态的有力指标。单细胞RNA测序通过定量的方式以高准确度、高灵敏度和高通量来揭示RNA丰度。然而，这种方法仅捕获细胞群某个时间点的静态快照，这对分析细胞群在时间维度的现象研究（例如胚胎发生或组织再生）构成了挑战。

作者提出了RNA速率，即基因表达状态的时间导数，在常见的单细胞RNA测序方案中，可以通过新生（未剪接）和成熟（剪接）mRNA的相对丰度来估计基因剪接和降解的速率，区分未剪接和剪接的mRNA，来直接估计细胞基因表达的动态分化。刚转录出的mRNA包含外显子和内含子，经过splicing切除内含子后，得到用于编码蛋白的spliced mRNA。spliced mRNA的丰度由未成熟mRNA的splicing速度和降解速率共同决定。通过计算未剪接转录本和剪接转录本之间的比率来推断细胞命运的状态（过渡与稳定）和方向性（轨迹）。

简而言之，RNA速率分析使我们能够使用转录动力学的数学模型推断在scRNA-seq实验中未直接观察到的转录动力学。我们可以使用RNA速率来确定在给定的目标细胞群中是否诱导或抑制了目标基因。此外，我们可以通过伪时间轨迹推断这些信息来预测细胞命运决定。

问题1：通过splicing和spliced mRNA信息如何推断分化方向？

在时间轴上，mRNA是未剪接到剪接的过程。scRNA-seq实验只是捕获细胞群某个时间点的静态快照，也就是细胞群一瞬间的转录组信息，时间点S1。假设分析A，B，C三种细胞类型，我们可以计算到在S1时间点，每个细胞类型的mRNA剪接率，ABC剪接率由低到高，怎么能得到分化方向是A->B->C这个结论?只是一个时间点，而不是一个时间轴，这块理解还是很断线...

分析步骤：

本教材将演示如何将处理过/标准化的Seurat对象与RNA速率分析结合使用。Seurat是基于R语言的，但所有RNA Velocity软件/包都是Python，因此我们将seurat单细胞数据迁移至python环境。

安装程序：

一、velocyto生成loom文件

首先，我们需要将Seurat分析中使用的每个单细胞样本生成loom文件（一种为基因组数据集设计的文件格式，例如单细胞）。loom文件与Seurat 中使用的文件格式不同。必须根据样品的原始FASTQ或BAM文件生成loom文件。生成Loom文件有两种方法是：Velocyto命令行和Kallisto Bustools(KB.)
Kallisto Bustools没有尝试过，我们用Velocyto实操。

velocyto github网站： https://github.com/velocyto-team/velocyto.py
velocyto文档： https://velocyto.org/velocyto.py

1.1 安装velocyto包的依赖

conda install numpy scipy cython numba matplotlib scikit-learn h5py click
1.2 安装velocyto
pip install velocyto
velocyto --help
1.3 示例脚本
velocyto run10x是针对10X的样本测序数据的命令，运行前要进行cellranger分析，cellranger分析完后生成bam文件。

velocyto将样本的bam文件生成loom文件，此步骤耗时，建议在HPC计算节点执行。

velocyto.py的命令行参数：参考
velocyto run10x -m repeat_msk.gtf mypath/sample01 somepath/refdata-cellranger-mm10-1.2.0/genes/genes.gtf
-s, --metadatatable FILE 样本的注释文件，csv格式，行是样本，列是注释信息。(跟seurat，monocle等R包中的metadata格式是一样的)
-m, --mask FILE 去掉一些重复表达的基因， .gtf文件，可以从UCSE gene browser上下载
-l, --logic TEXT 用于筛选的逻辑值（这个参数还没有研究过）
-M, --multimap 考虑非唯一映射（作者不建议这样）
-@, --samtools-threads INTEGER 用samtools对bam文件进行sort时用到的进程数
--samtools-memory INTEGER 每个进程分配的内存（单位Mb）
-t, --dtype TEXT loom文件层的dtype。如果每个细胞的每个基因超过6000个molecules/reads，这个参数建议设置为32。（默认是16）
mypath/sample01：sample01的cellranger结果目录，该目录包含outs、outs/analys和outs/filtered_feature_bc_matrix子文件夹，包含bam文件。

genes.gtf：cellranger分析流程的参考基因组的注释文件

path/cellranger/reference/refdata-gex-mm10-2020-A/genes/genes.gtf

path/cellranger/reference/refdata-gex-GRCh38-2020-A/genes/genes.gtf

repeat_msk.gtf：重复表达的基因.gtf文件，可以从UCSE gene browser上下载

human_allTracks.gtf：小鼠的链接

mm10_allTracks.gtf：小鼠的链接
百度云盘下载地址：

path/to/RNA_velocity/human_allTracks.gtf

path/to/RNA_velocity/mm10_allTracks.gtf

链接：https://pan.baidu.com/s/1-dk1CQUa-39UHIxJ0wWQpg 提取码：1234
1.4 运行示例项目
# mouse
velocyto run10x -m path/to/RNA_velocity/mm10_allTracks.gtf /path/to/scrna_project/Cellranger_result/sample1/ path/to/cellranger/reference/refdata-gex-mm10-2020-A/genes/genes.gtf -@ 10 --samtools-memory 24000
# human
velocyto run10x -m path/to/RNA_velocity/human_allTracks.gtf /path/to/scrna_project/Cellranger_result/sample1/ path/to/cellranger/reference/refdata-gex-GRCh38-2020-A/genes/genes.gtf -@ 10 --samtools-memory 24000
运行上面的程序老是出现报错：
no such file or directory: 'samtools'
出现的问题说明：没有安装samtools和libcrypto.so.1.0.0
安装samtools
conda install samtools=1.5
samtools出现的报错
samtools: error while loading shared libraries: libcrypto.so.1.0.0
尝试过的方法：
conda install -c bioconda samtools openssl=1.0
conda install samtools openssl=1.0
最终的解决方法：
conda config --add channels https://mirrors.bfsu.edu.cn/anaconda/pkgs/main
conda config --add channels https://mirrors.bfsu.edu.cn/anaconda/pkgs/free
conda install samtools openssl
samtools
重新运行，成功，写成python脚本去集群投递任务。
二、从seurat对象中提取meta-data
（R语言运行）已经在R中使用Seurat执行了主要的数据处理（过滤、归一化、聚类、批次处理、降维，细胞类型注释），接下来，提取meta-data的主要信息：
Filtered Cell Ids
UMAP or TSNE coordinates
Clusters/Celltype (Optional)
Clusters/Celltype Colors (Optional)
makedirs <- function(dir_list){ 
  for (i_dir in dir_list){ 
    if (!dir.exists(i_dir) ){ 
      dir.create(i_dir, recursive = TRUE) 
#--------------------------------------------------
makedirs(file.path(project_work_dir, "rna_velocity")) 
setwd(file.path(project_work_dir, "rna_velocity"))
table(seurat_obj$celltype)
# 获得每个细胞的UMAP或TSNE坐标，使用 Embeddings函数
write.csv(Embeddings(seurat_obj, reduction = "umap"), file = "cell_embeddings.csv")
# 获取每个细胞的barcode
write.csv(Cells(seurat_obj), file = "cellID_obs.csv", row.names = FALSE)
# 提取每个细胞的cluster信息
write.csv(seurat_obj@meta.data[, 'cluster', drop = FALSE], file = "cell_clusters.csv")
# 提取每个细胞的celltype信息
write.csv(seurat_obj@meta.data[, 'celltype', drop = FALSE], file = "cell_celltype.csv")
# 获取celltype的颜色信息
hue_pal()(length(levels(seurat_obj$celltype)))
# 获取cluster的颜色信息
hue_pal()(length(levels(seurat_obj$cluster)))
# 绘制umap图，与RNA速率图对比看
Idents(seurat_obj) <- "celltype"
alpha.use <- 0.8
p1 <- DimPlot(object = seurat_obj, reduction = "umap", label = TRUE, label.size = 5, pt.size=0.4, raster=FALSE)
p1$layers[[1]]$mapping$alpha <- alpha.use
save_plot("UMAP_clustering.png", p1, base_height = 8, base_aspect_ratio = 1.3, base_width = NULL, dpi=600)
save_plot("UMAP_clustering.pdf", p1, base_height = 8, base_aspect_ratio = 1.3, base_width = NULL)
seurat对象中的barcode是有样本名称前缀的，sample1_, sample2_...
三、整合Loom文件和meta-data
（python运行）整合loom文件和meta-data数据

代码中有对loom文件barcode的处理，与seurat的barcode一致。
# 整合多个loom文件
import loompy
files=['path/to/sample_one','path/to/sample_two']
output_filename='/save/path/to/combined.loom'
loompy.combine(files, output_filename, key="Accession")
# 整合loom文件和meta-data数据
import scvelo as scv
import pandas as pd
import numpy as np
import os
loom_data = scv.read('/save/path/to/combined.loom', cache=False)
loom_data.obs
# barcode名字去重后缀x，与seurat导出的barcode名称一致
loom_data.obs = loom_data.obs.rename(index = lambda x: x.replace(':', '-').replace('x', ''))
loom_data.obs.head()
# 读取seurat中的meta信息
meta_path = "path/to/rna_velocity"
sample_obs = pd.read_csv(os.path.join(meta_path, "cellID_obs.csv"))
cell_umap= pd.read_csv(os.path.join(meta_path, "cell_embeddings.csv"), header=0, names=["Cell ID", "UMAP_1", "UMAP_2"])
cell_clusters = pd.read_csv(os.path.join(meta_path, "cell_clusters.csv"), header=0, names=["Cell ID", "cluster"])
cell_celltype = pd.read_csv(os.path.join(meta_path, "cell_celltype.csv"), header=0, names=["Cell ID", "celltype"])
# 对细胞文件和RNA剪切速率文件取交集，保留关注的细胞类型
sample_one = loom_data[np.isin(loom_data.obs.index, sample_obs)]
sample_one.obs.head()
# 构建umap坐标, cluster, celltype信息数据框
sample_one_index = pd.DataFrame(sample_one.obs.index)
sample_one_index = sample_one_index.rename(columns = {0:'Cell ID'})
umap_ordered = sample_one_index.merge(cell_umap, on = "Cell ID")
umap_ordered.head()
celltype_ordered = sample_one_index.merge(cell_celltype, on = "Cell ID")
celltype_ordered.head()
clusters_ordered = sample_one_index.merge(cell_clusters, on = "Cell ID")
clusters_ordered.head()
# 将umap坐标与cluster信息加入sample_one
umap_ordered = umap_ordered.iloc[:,1:]
clusters_ordered = clusters_ordered.iloc[:,1:]
celltype_ordered = celltype_ordered.iloc[:,1:]
sample_one.obsm['X_umap'] = umap_ordered.values
sample_one.uns['clusters'] = clusters_ordered.values
sample_one.obs['celltype'] = celltype_ordered.values
adata = sample_one
# some gene labels are duplicated (Ensembl IDs are still unique!!)
adata.var_names_make_unique()
# save model to file
adata.write('Allcelltype_dynamicModel.h5ad', compression = 'gzip')
# read h5ad file
adata= scv.read('Allcelltype_dynamicModel.h5ad')
四、运行RNA速率
（python运行）可视化绘图
# Running RNA Velocity
scv.pp.filter_and_normalize(adata,min_shared_counts=30, n_top_genes=2000)
scv.pp.moments(adata, n_pcs=30, n_neighbors=30)
scv.tl.velocity(adata, mode = "stochastic")
scv.tl.velocity_graph(adata)
scv.pl.velocity_embedding(adata, basis='X_umap', arrow_size=5)
ident_colours = ["#F8766D", "#A3A500", "#00BF7D", "#00B0F6", "#E76BF3"]
scv.pl.velocity_embedding_stream(adata, basis='X_umap',color = "celltype", palette = ident_colours)
scv.pl.velocity_embedding_stream(adata, basis='X_umap',color = "celltype", palette = ident_colours, size = 20,alpha =0.8)
scv.pl.velocity_embedding_stream(adata, basis='X_umap',color = "celltype", palette = ident_colours, size = 20,alpha =0.8, save="UMAP_stream.png", figsize=(7,5), legend_fontsize = 9, show=False, title='')
scv.pl.velocity_embedding_stream(adata, basis='X_umap',color = "celltype", palette = ident_colours, size = 20,alpha =0.8, save="UMAP_stream.pdf", figsize=(7,5), legend_fontsize = 9, show=False, title='')
以上是我常用的分析流程，下面学习下smorabit的教材，看哪些内容可以借鉴。
Morabito教程： RNA velocity analysis with scVelo
Step-1：数据从Seurat 转换为Python/anndata
将Seurat对象转换为可在scVelo中可使用的数据文件，共导出4个文件，分别是：1）seurat对象的metadata文件；2）细胞的基因表达矩阵counts.mtx，文件很大；3）pca降维的pca坐标文件pca.csv，有多列，选择了多少PCA，就有多少列。4）基因名称列表gene_names.csv；
makedirs(file.path(project_work_dir, "rna_velocity")) 
setwd(file.path(project_work_dir, "rna_velocity"))
table(seurat_obj$celltype)
# save metadata table:
seurat_obj$barcode <- colnames(seurat_obj)
seurat_obj$UMAP_1 <- seurat_obj@reductions$umap@cell.embeddings[,1]
seurat_obj$UMAP_2 <- seurat_obj@reductions$umap@cell.embeddings[,2]
write.csv(seurat_obj@meta.data, file='metadata.csv', quote=F, row.names=F)
# write expression counts matrix
library(Matrix)
counts_matrix <- GetAssayData(seurat_obj, assay='RNA', slot='counts')
writeMM(counts_matrix, file='counts.mtx')
# write dimesnionality reduction matrix, in this example case pca matrix
write.csv(seurat_obj@reductions$pca@cell.embeddings, file='pca.csv', quote=F, row.names=F)
# write gene names
write.table(
  data.frame('gene'=rownames(counts_matrix)),file='gene_names.csv',
  quote=F,row.names=F,col.names=F
# load gene names:
with open("gene_names.csv", 'r') as f:
    gene_names = f.read().splitlines()
# set anndata observations and index obs by barcodes, var by gene names
adata.obs = cell_meta
adata.obs.index = adata.obs['barcode']
adata.var.index = gene_names
# load dimensional reduction:
pca = pd.read_csv("pca.csv")
pca.index = adata.obs.index
# set pca and umap
adata.obsm['X_pca'] = pca.to_numpy()
adata.obsm['X_umap'] = np.vstack((adata.obs['UMAP_1'].to_numpy(), adata.obs['UMAP_2'].to_numpy())).T
# plot a UMAP colored by sampleID to test:
sc.pl.umap(adata, color=['celltype'], frameon=False, save=True)
# save dataset as anndata format
adata.write('my_data.h5ad')
# reload dataset
adata = sc.read_h5ad('my_data.h5ad')
Step 2: 构建spliced和unspliced计数矩阵
与我们在Seurat中使用的相同的基于UMI的基因计数矩阵不同，我们需要一个splicing和spliced转录本的矩阵。 我们可以使用velocy命令行工具或使用Kallisto-Bustools构建这个矩阵。在这里我使用了  velocyto command line tool,，仅仅是因为我个人从来没有尝试过Kallisto-Bustools。
velocy命令行工具直接从cellranger输出目录开始运行，需要提供.bam文件，它也可以在任何单细胞平台上使用。我们还必须为你的物种提供参考.gtf注释（此处使用mm10），并且你可以选择提供.gtf以掩盖重复区域（velocy推荐这种方式）。
repeats="/path/to/repeats/mm10_rmsk.gtf"
transcriptome="/path/to/annoation/file/gencode.vM25.annotation.gtf"
cellranger_output="/path/to/cellranger/output/"
velocyto run10x -m $repeats \
                   $cellranger_output \
                   $transcriptome
Step 3: 加载数据
前面我们已经生成了scvelo的输入数据，可以将其加载到python中。Velocy为每个样本创建了一个单独的拼接和未拼接矩阵，因此我们首先必须将不同的样本合并为一个对象。此外，需要处理loom文件的细胞barcode名称，以使anndata对象与基因表达矩阵数据相匹配。
import scvelo as scv 
import scanpy as sc 
import cellrank as cr 
import numpy as np 
import pandas as pd 
import anndata as ad 
scv.settings.verbosity = 3 
scv.settings.set_figure_params('scvelo', facecolor='white', dpi=100, frameon=False) 
cr.settings.verbosity = 2 
adata = sc.read_h5ad('my_data.h5ad') 
# load loom files for spliced/unspliced matrices for each sample: 
ldata1 = scv.read('Sample1.loom', cache=True) 
ldata2 = scv.read('Sample2.loom', cache=True) 
ldata3 = scv.read('Sample3.loom', cache=True) 
Variable names are not unique. To make them unique, call `.var_names_make_unique`. 
Variable names are not unique. To make them unique, call `.var_names_make_unique`. 
Variable names are not unique. To make them unique, call `.var_names_make_unique`. 
# rename barcodes in order to merge: 
barcodes = [bc.split(':')[1] for bc in ldata1.obs.index.tolist()] 
# 查看seurat的barcode命名方式，也可以是'Sample1_'
barcodes = [bc[0:len(bc)-1] + '_10' for bc in barcodes] 
ldata1.obs.index = barcodes 
barcodes = [bc.split(':')[1] for bc in ldata2.obs.index.tolist()] 
barcodes = [bc[0:len(bc)-1] + '_11' for bc in barcodes] 
ldata2.obs.index = barcodes 
barcodes = [bc.split(':')[1] for bc in ldata3.obs.index.tolist()] 
barcodes = [bc[0:len(bc)-1] + '_9' for bc in barcodes] 
ldata3.obs.index = barcodes 
# make variable names unique 
ldata1.var_names_make_unique() 
ldata2.var_names_make_unique() 
ldata3.var_names_make_unique() 
# concatenate the three loom 
ldata = ldata1.concatenate([ldata2, ldata3]) 
# merge matrices into the original adata object 
adata = scv.utils.merge(adata, ldata) 
# plot umap to check 
sc.pl.umap(adata, color='celltype', frameon=False, legend_loc='on data', title='', save='_celltypes.pdf')
Step 4: 使用scVelo计算RNA速率
我们使用scVelo来计算RNA速率。scVelo允许用户使用2018年原始出版物中的稳态模型（steady-state model）以及 2020年出版物中更新的动态模型（ dynamical model）。
首先，我们检查每个细胞类型 中spliced和unspliced的占比。接下来我们执行一些预处理，然后我们使用稳态模型（随机选项）计算RNA速率。
scv.pl.proportions(adata, groupby='celltype_full')
scv.pp.filter_and_normalize(adata) 
scv.pp.moments(adata) 
WARNING: Did not normalize X as it looks processed already. To enforce normalization, set `enforce=True`. 
Normalized count data: spliced, unspliced. 
WARNING: Did not modify X as it looks preprocessed already. 
computing neighbors 
    finished (0:00:14) --> added 
    'distances' and 'connectivities', weighted adjacency matrices (adata.obsp) 
computing moments based on connectivities 
    finished (0:00:17) --> added 
    'Ms' and 'Mu', moments of un/spliced abundances (adata.layers)
# compute velocity 
scv.tl.velocity(adata, mode='stochastic') 
scv.tl.velocity_graph(adata) 
computing velocities 
    finished (0:01:15) --> added 
    'velocity', velocity vectors for each individual cell (adata.layers) 
computing velocity graph 
    finished (0:15:00) --> added 
    'velocity_graph', sparse matrix with cosine correlations (adata.uns)
Step 5:  可视化RNA速率场
我们可以通过可视化二维降维顶部的矢量场，对所有基因和所有细胞的RNA速率进行广泛的可视化。
scv.pl.velocity_embedding(adata, basis='umap', frameon=False, save='embedding.pdf')
computing velocity embedding
    finished (0:00:06) --> added
    'velocity_umap', embedded velocity vectors (adata.obsm)
saving figure to file ./figures/scvelo_embedding.pdf
image
我们还可以可视化我们感兴趣的基因的动态变化。在这里，我们展示了Ptgds基因的未剪接转录本与剪接转录本的比率，以及在feature plots图中显示基因的速率和基因表达值。
# plot velocity of a selected gene
scv.pl.velocity(adata, var_names=['Ptgds'], color='celltype')