velocyto | scRNA速率计算方法

初学者还不太了解其中原理与步骤，谨记录本人的操作流程与学习心得，抛砖引玉供赏之。

官网教程：

Tutorial — velocyto 0.17.16 documentationhttp://velocyto.org/velocyto.py/tutorial/

大致流程：

1、将bam文件转化为loom文件

2、将loom文件转化为seruat对象

3、对数据进行标准化降维聚类分析

4、计算scRNA速率

5、展示分化方向

安装操作：

1、配置环境：（官网上提供几种下载方法，我选的是Install using PyPI）

conda install numpy scipy cython numba matplotlib scikit-learn h5py click

2、install velocyto：

pip install velocyto

#测试安装是否成功
velocyto --help

操作：

1、下载FASTA文件、GTF文件：从GENCODE(GENCODE - Home page)或者Ensembl(FTP Download)中下载基因注释文件。

GTF文件格式：（第三列-特征类型-外显子；第九列：转录名或基因ID或基因名或外显子编号）

 2、run之前输入的bam文件需要按照位置进行排序。（如果是cellranger的话，则已经执行这一步）

samtools sort -t CB -O BAM -o cellsorted_possorted_genome_bam.bam possorted_genome_bam.bam

3、根据不同产生数据的方式选择不同的run方式。run10x、run_smartseq2、run_dropest三种命令都隶属于run的命令亚组，因此也可直接用通用模式run。

#三种不同方式（10X、smartseq2、dropest）产生的数据进行处理
velocyto run10x -m repeat_msk.gtf mypath/sample01 somepath/refdata-cellranger-mm10-1.2.0/genes/genes.gtf

velocyto run_smartseq2 -o OUTPUT -m repeat_msk.gtf -e MyTissue plateX/*/*.bam mm10_annotation.gtf

velocyto run_dropest -o ~/mydata/SRR5945695_results -m rep_mask.gtf ~/mydata/SRR5945695_1/correct_SRR5945695_1Aligned.sortedByCoord.out.tagged.bam mm10_annotation.gtf

#通用run
velocyto run -b filtered_barcodes.tsv -o output_path -m repeat_msk_srt.gtf possorted_genome_bam.bam mm10_annotation.gtf
         

4、生成 .loom文件 并整合为一个.loom文件

files = ["Spleen_D0.loom","Spleen_D1.loom"]

loompy.combine(files,"merged.loom",key="Accession")

5、结合.loom文件，对原始数据进行提取

B.integrated <- readRDS("./B.cell.integrated.anno.rds")

#提取名字，并改成与.loom文件中的细胞名一致，并保存
cells.renamed <- Cells(B.integrated)
cells.renamed <- gsub("([A-Z]+)-1_1","Spleen_D0:\\1x",cells.renamed)
cells.renamed <- gsub("([A-Z]+)-1_2","Spleen_D1:\\1x",cells.renamed)
write.csv(cells.renamed, file = "./cellID_obs.csv", row.names = FALSE)

#提取umap及聚簇数据，改名，并保存
embeddings.renamed <- Embeddings(HSPC.integrated, reduction = "umap")
rownames(embeddings.renamed) <- cells.renamed
write.csv(embeddings.renamed, file = "./cell_embeddings.csv")
cell_clusters <- as.matrix(HSPC.integrated@meta.data$MajorType)
rownames(cell_clusters) <- cells.renamed
write.csv(cell_clusters, file = "./clusters.csv")

6、通过scVelo可视化RNA速率

RNA速率计算的三种模型：

1. 稳态/确定性模型（The steady-state / deterministic model）：在转录阶段（诱导和抑制）足够长时间来达到稳定平衡状态的假设下，把稳态下的比率与观测状态下的比率偏差作为RNA速率的量化值。mRNA的平衡水平可以大致看作是假定稳定状态的上下分数的线性回归，这就产生了两种基本的假设，即基因间的共同剪接率与稳定状态下的mRNA水平。       缺点：当一个种群包含多种异质亚种群状态时，速率估计与细胞状态可能会出错。
2. 随机模型（The stochastic model）：将转录、剪接、降解作为概率事件，用力矩方程（moment equations）来估计产生的马尔科夫过程。在二阶矩阵中，利用了剪接与未剪接的mRNA水平和它们的方差来构建随机模型，能够更好地捕获稳定状态。
3. 动态模型（The dynamical model）：（计算量最大、最强）不依赖于普通剪接率或稳定状态的限制来取样，使RNA速率能够适应广泛变化，对每个基因剪接动力学的完整动态进行捕获。基于最大似然法，迭代估计反应速率和潜在细胞特异性变量的参数。

scVelo包中的基本命令：数据预处理的调用（scv.pp.*），分析工具（scv.tl.*），绘图工具（scv.pl.*）。官网上有关于其基本用法。

#基于Python3.6

pip install -U scvelo
# pip install -U scvelo --user

import scvelo as scv
import pandas as pd
import anndata
import numpy as np
from matplotlib import pyplot as plt

adata = scv.read(merged.loom, cache=True)
sample_obs = pd.read_csv("./cellID_obs.csv")  

# 提取与sample_obs中ID相同的adata
adata_a = adata[np.isin(adata.obs.index,sample_obs["x"])]

umap & Majortype 由前面提取出来的embedding & clusters 筛选匹配

# 在基因选择与数据标准化后，计算RNA速率的一阶矩和二阶矩（均值与非中心方差）
scv.pp.filter_and_normalize(adata_a)
scv.pp.moments(adata_a, n_neighbors=30)
scv.tl.velocity(adata_a, mode = "steady_state")

# scv.tl.velocity(adata_a, mode = "stochastic")   #随机模型
## scv.tl.recover_dynamics(adata_a, max_iter=100, return_model=True)   #动力学模型
## scv.tl.velocity(adata_a, mode='dynamical')

scv.tl.velocity_graph(adata_a)
scv.pl.velocity_embedding_stream(adata_a, basis='umap',color=['MajorType'], save = './B_scvelo_steady_state.png')

scv.tl.paga(adata_a, groups='MajorType')
scv.pl.paga(adata_a, basis='umap', size=50, alpha=.1, min_edge_width=2, node_size_scale=1.5, save = './B_scVelo_dynamical_paga.pdf')

top_genes = adata_a.var.sort_values('fit_likelihood', ascending=False).index[:10]
scv.pl.scatter(adata_a, basis=top_genes,color=['MajorType'],ncols=5, save = './B_scVelo_dynamical_graph.pdf')