使用scVelo进行RNA velocity分析-01

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文件，每个样本单独生成；这步前提要cellranger分析；
• 第二步：获取seurat分析中的meta-data数据，一般只对关注的某些细胞类型做RNA速率分析，获取umap坐标信息，cluster和celltype信息；这步前提要seurat聚类分析；
• 第三步：运行python程序，调用scvelo，整合loom文件和meta-data数据，运行RNA速率，绘图；
  如果有多个样本，需要整合多个样本的loom文件，生成combined.loom；

安装程序：

• scVelo（用于RNA Velocity）
• Velocyto 或 Kallisto Bustools（用于生成我们的初始RNA Velocity对象）
• Anndata（用于操作我们的RNA Velocity对象）
• Seurat
• Samtools -- 可选（Velocyto运行Samtools排序未分类.bam)

---

一、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_...

image.png

三、整合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
)

image.png

（python运行）seurat数据对象转adata数据结构数据

import scanpy as sc
import anndata
from scipy import io
from scipy.sparse import coo_matrix, csr_matrix
import numpy as np
import os
import pandas as pd

# load sparse matrix:
X = io.mmread("counts.mtx")

# create anndata object
adata = anndata.AnnData(
    X=X.transpose().tocsr()
)

# load cell metadata:
cell_meta = pd.read_csv("metadata.csv")

# 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')

image.png

Step 4: 使用scVelo计算RNA速率

我们使用scVelo来计算RNA速率。scVelo允许用户使用2018年原始出版物中的稳态模型（steady-state model）以及 2020年出版物中更新的动态模型（ dynamical model）。

首先，我们检查每个细胞类型 中spliced和unspliced的占比。接下来我们执行一些预处理，然后我们使用稳态模型（随机选项）计算RNA速率。

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

image.png

# pre-process 
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.png

scv.pl.velocity_embedding_grid(adata, basis='umap', color='celltype', save='embedding_grid.pdf', title='', scale=0.25)

image

scv.pl.velocity_embedding_stream(adata, basis='umap', color=['celltype', 'condition'], save='embedding_stream.pdf', title='')

image

我们还可以可视化我们感兴趣的基因的动态变化。在这里，我们展示了Ptgds基因的未剪接转录本与剪接转录本的比率，以及在feature plots图中显示基因的速率和基因表达值。

# plot velocity of a selected gene
scv.pl.velocity(adata, var_names=['Ptgds'], color='celltype')

image.png

---

Morabito教程是把seurat转成python环境的Anndata对象，再和loom数据进行merge；
而我们常见的做法是将seurat的meta-data和loom数据进行merge；
后面的下游分析和对RNA速率算法理解留到下一篇文章。未完待续...