从零开始的RNA Velocity

发表于2022-12-20更新于2023-03-10

RNA Velocity目前已经成为比较高级的单细胞动态分析，目前已有的方法，主要是最初的Velocyto和改良的scVelo，引入了metabolite label seq的dynamo，基于多组学的MultiVelo和在空间转录组层面应用的SIRV。目前基于RNA Velocity的单细胞动力学分析是最为可靠的单细胞层面的细胞动态分析，基于此还有一些衍生工具，比如CellRank，可以使用RNA Velocity数据推断细胞发育路径（这个其实挺迷惑，毕竟Velocity的信息已经能够很充分地进行细胞发育路径推断了。）

此外，由于真实的Velocity分析的计算开销非常恐怖，因此出现了一些基于纯转录组数据进行推断的工具，这些工具主要有两个方面的作用，一是向量场的推断，也就是全局层面的RNA Velocity，二是动力学推断，也即单个细胞的Velocity Vector的计算和推断。这些工具有VECTOR、scTour和DeepVelo，前者是一个R的脚本工具，后者是一个基于深度学习的python工具，二者从结果上讲，scTour更为出色一些，因为除了向量场推断之外，还顺便做了latent space的学习，可以理解为顺手做了个拟时序。除此之外，诸如PHATE降维，SPRING降维实际上已经有一些细胞向量方向上推断的意思了，只是当时还没有这个概念。

然而RNA Velocity的分析似乎长期处于一个教程不足的状态，原因很多，但最主要的是因为其涉及到多环境切换，如果要完完整整从FASTQ文件做一份RNA Velocity分析，需要涉及linux、python和R，目前来讲可以说是最为全面的分析范式，涉及到单细胞分析的全部方面是，所有数据，各种数据类型和全部分析步骤，如果只是为了得到一个简单的Velocity的结果，似乎有点不划算，因此通常在一开始就要决定好是否做这个分析，否则后期补的话那可以说是非常痛苦。

从各个工具的教程来看，典型的start都是给一个loom文件，然后用各种工具读入，然后就开始进行velocity的计算了，这个loom文件中，稍微知道一点的人就能发现，预处理、降维聚类、细胞注释乃至去批次，全都做好了，直接进行velocity的计算就行，但实际上呢？自己的数据，甚至从公司拿到的数据都不可能直接给一个loom文件，或者给一个loom文件但缺少一些关键信息，总之就是缺这缺那，网上各种教程都喜欢直接用官方给的做好的文件，这种意义不大。

那么，这个loom文件是怎么来的？哪些信息最为重要？这是最关键的两个问题：

  • loom文件从上游的bam文件获得，可以使用3个pipeline：Velocyto，loompy和kallisto+bustool
  • 最关键的信息：就是spliced/unspliced的metric，其余所有信息都是次要的，在其他工具里使用/获取到的

因此，这两个方面的信息怎么来？

首先，以10X提供的最为经典的pbmc3k的数据为例子，在10X的官网上可以下载到全部的数据，包括FASRQ，上游的bam文件，下游的10X标准输入文件等，这里从bam开始，因此最麻烦的就是得把这里的bam文件下载下来，如果是完全自己跑上游，那么在10X的outs文件夹里就有这个bam，直接使用即可。此外，由于需要对每一个细胞评估splice的情况，因此还需要barcodes.tsv这个文件，而且需要是filtered_gene_bc_matrices里面的barcodes。

需要特别注意的是，使用velocyto的时候，由于路径都是写死的，因此文件夹组织形式和文件名都是有特殊要求的，对于单个project，应当按照下面的文件夹组织，注意如果不会改动源码的话，连文件的名字都要是一样的：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
./pbmc3k
├── outs
│   ├── cellsorted_possorted_genome_bam.bam
│   ├── filtered_gene_bc_matrices
│   │   └── hg19
│   │       ├── barcodes.tsv
│   │       ├── genes.tsv
│   │       └── matrix.mtx
│   ├── possorted_genome_bam.bam
│   ├── possorted_genome_bam_index.bam.bai
│   └── raw_gene_bc_matrices
│       └── hg19
│           ├── barcodes.tsv
│           ├── genes.tsv
│           └── matrix.mtx
├── pbmc3k_filtered_gene_bc_matrices.tar.gz
├── pbmc3k_raw_gene_bc_matrices.tar.gz
└── velocyto
    └── pbmc3k.loom

如果有多个项目，则每一个都要按照这个组织形式，存放在不同的projectname的文件夹里，然后才能使用，如果是smart-seq或者drop-seq的数据，则同样需要bam和类似的组织形式，只是10X用的最多罢了。

这里能看到有个velocyto文件夹，这个里面的pbmc3k.loom就是最后输出的结果了。最开始是没有这个文件夹的。

除了单细胞的数据，由于这里需要切分RNA的splice情况，因此需要基因组的注释文件，其一是在做上游的时候，由10X官方提供的reference里的gene/gene.gtf，这个文件也可以从ensembl网站自己下载，只是但凡跑了上游则必然存在这个文件，容易获得。

第二个注释文件则是masked sequence的gtf注释文件，这个文件可以从GENECODE、ensembl网站获取，但最方便的则是从UCSC Browser中获取，从Table Browser中进行获取下载，需要注意选择好物种、基因组版本、mask sequence选项，输入文件名并使用gzip下载，然后解压就能获得这个文件了。

然后使用velocyto就能获得这个上游的loom文件了。velocyto的安装也有点坑爹，mac和windows上安装非常麻烦，千万别用pip安装，涉及到编译，由于目前的0.17.7版本是19年上传的，如果编译系统太新了还没法编译，非常烦。直接使用github上的源码安装就行了。

1
2
3
4
5
#依赖包安装，建议用conda或者mamba安装
conda install numpy scipy cython numba matplotlib scikit-learn h5py click
git clone https://github.com/velocyto-team/velocyto.py.git
cd velocyto.py
pip install -e .  # note the trailing dot

这里只要没有缺依赖库，基本都能装上，连windows都装得上去，然后激活对应的conda环境就能直接使用velocyto了。

实际上，在clone下来velocyto.py后，进入/velocyto.py/velocyto/commands/ 里面就能看到使用的的命令了，如果想改动上游读取数据的方式，在对应的命令里就能进行改动了，比如希望使用全部的raw barcodes而不是filtered之后的barcodes，在run10x.py里改动对应的文件夹名字就好了，或者笨一点的办法是把文件夹的名字改了，只要自己记得就行了。

然后就能使用velocyto命令，生成这里的loom文件了：

1
2
3
4
5
rmsk_gtf=$YOUR_REF_DIR/hg38_repeat_rmsk.gtf
cellranger_outDir=pbmc3k
cellranger_gtf=$YOUR_REF_DIR/refdata-gex-GRCh38-2020-A/genes/genes.gtf

nohup velocyto run10x -m $rmsk_gtf  $cellranger_outDir $cellranger_gtf &

然后等着就行了…..个人的渣电脑跑了大概4小时，bam数据15.6G，最后输出一个20+MB的loom…….这一步极其消耗时间，但是一个多线程的步骤，因此建议使用集群运行，大大节省时间空间。

有时候上游的bam没排序或者忘了，则需要手动先排序：

1
2
cd $cellranger_outDir/out 
nohup samtools sort -@ 10  -t CB -O BAM -o cellsorted_possorted_genome_bam.bam possorted_genome_bam.bam &

sort之后再进行运行即可，但通常只要samtools安装好了是不需要手动做sort这一步的。

多个样品的话，需要对每一个样品进行同样的操作，可以写bash的循环先手动sort：

1
2
3
4
ls  */outs/possorted_genome_bam.bam|while read id;do  new=${id/possorted_genome_bam.bam/cellsorted_possorted_genome_bam.bam}
echo $new 
nohup samtools sort -@ 4  -t CB -O BAM -o $new   $id  & 
 done

再批量跑velocyto即可,然后每个project都会有一个velocyto文件夹，里面放着对应的loom文件

1
2
3
ls -d *-*|while read cellranger_outDir;do 
nohup velocyto run10x -m $rmsk_gtf  $cellranger_outDir $cellranger_gtf & 
done

如果使用kallisto+bustools，则相对命令上简单一些，这个pipeline支持直接从上游生成loom文件：

1
2
3
4
5
kb ref -i index.idx -g t2g.txt -f1 cdna.fa -f2 intron.fa -c1 cdna_t2c.txt -c2 intron_t2c.txt --workflow lamanno -n 4 \
fasta.fa \
gtf.gtf

kb count -i transcriptome.idx -g t2g.txt -x 10xv2 --workflow lamanno --loom -c1 cdna_t2c.txt -c2 intron_t2c.txt read_1.fastq.gz read_2.fastq.gz

运行完kb count就会直接出来loom，只要加了—loom参数即可。

这时候，已经距离官方的示例文件SCG71.loom很近了，下游可以跑了。由于velocyto是最原始的方法，结果上不是很稳定，因此现在通常只用它来获取输入，下游采用scVelo或者dynamo进行分析。

scVelo直接读这里生成的loom即可，如果有多个loom则挨个读进去，然后使用scv.utils.merge(adata, ldata……)merge为一个anndata即可，这里只有一个loom，因此暂不需要merge，merge之后的流程和单个文件流程完全一致。

1
2
3
4
5
6
7
8
9
10
11
12
import scvelo as scv
import scanpy as sc

scv.logging.print_version()
scv.settings.verbosity = 3  # show errors(0), warnings(1), info(2), hints(3)
scv.settings.set_figure_params('scvelo')  # for beautified visualization

loomf = './pbmc3k/velocyto/pbmc3k.loom'
adata = scv.read(loomf, cache=False)

adata
adata.var_names_make_unique

这里如果熟悉scanpy的会发现，adata中只有两个layer，什么X_umap，PCA降维所有的信息一概没有，因为我们读入的，的确只是一个原始的loom文件，里面的确不包含任何降维聚类注释等信息，因此，这里有两种办法，其一是在进行velocity之前，使用scanpy的流程，把降维聚类先做了，至少要获得PCA、umap/tsne和cell cluster信息，细胞是否进行注释则并不需要，但通常都习惯于直接注释，再继续用scvelo读入即可无缝继续分析，这时候这个loom文件就跟官方提供的SCG71.loom一致，包含了所有的信息，直接开跑velocity，不会有任何障碍，个人推荐采用这个操作，因为这样无论是筛细胞还是跑降维聚类，都是一致的。

第二种则是使用外部工具，比如seurat处理之后，提取关键信息后传入，我们需要提取的信息主要有细胞的barcode，因为Seurat可能对细胞进行了筛选，而且本身二者细胞数就可能是不一致的；第二则是umap的坐标，第三则是cell的annotation，或者其他任何的meta信息。但这样有代价，非常的麻烦，由于python生态对于barcode的处理不太一样，因此需要手动转换掉全部的barcode，并且进行匹配，然后再将meta信息和embedding信息读入后传入anndata，其实还挺麻烦的，当然也可以直接将loom读入R中，针对splice矩阵进行正常的降维聚类分析，然后使用sceasy再将seuratobject转为anndata…..这就有点迷惑操作了，直接用scanpy做不就好了。

1
2
3
write.csv(Cells(seurat_object), file = "cellID_obs.csv", row.names = FALSE)
write.csv(Embeddings(seurat_object, reduction = "umap"), file = "cell_embeddings.csv")
write.csv(seurat_object@meta.data$seurat_clusters, file = "clusters.csv")

然后就得在py里操作了，如果这里涉及到多个样本，通常会在barcode前面加诸如sample1_, sample2_，这样的前缀用来切分所有的细胞文件，anndata里需要单独映射

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
import anndata
import scvelo as scv
import pandas as pd
import numpy as np
import matplotlib as plt
%load_ext rpy2.ipython

sample_one = anndata.read_loom("sample_one.loom")
.... 
sample_n = anndata.read_loom("sample_n.loom")

sample_obs = pd.read_csv("cellID_obs.csv")
umap_cord = pd.read_csv("cell_embeddings.csv")
cell_clusters = pd.read_csv("clusters_obs.csv")

cellID_obs_sample_one = cellID_obs[cellID_obs_sample_one[0].str.contains("sample1_")]
cellID_obs_sample_two = cellID_obs[cellID_obs_sample_two[0].str.contains("sample2_")]
...
cellID_obs_sample_n = cellID_obs[cellID_obs_sample_two[0].str.contains("samplen_")]

sample_one = sample_one[np.isin(sample_one.obs.index, cellID_obs_sample_one)]
sample_two = sample_two[np.isin(sample_two.obs.index, cellID_obs_sample_two)]
...
sample_n = sample_n[np.isin(sample_n.obs.index, cellID_obs_sample_n)]

sample_one = sample_one.concatenate(sample_two, sample_three, sample_four,...)
umap = pd.read_csv("umap.csv")

sample_one_index = pd.DataFrame(sample_one.obs.index)
sample_one_index = sample_one_index.rename(columns = {0:'Cell ID'})

umap = umap.rename(columns = {'Unnamed: 0':'Cell ID'})

umap_ordered = sample_one_index.merge(umap, on = "Cell ID")

umap_ordered = umap_ordered.iloc[:,1:]
sample_one.obsm['X_umap'] = umap_ordered.values

sample_one.uns['Cluster_colors']

#继续velocity，本质上核心就是一个个做好barcode的映射在连接成同一个anndata
scv.pp.filter_and_normalize(sample_one)
scv.pp.moments(sample_one)
scv.tl.velocity(sample_one, mode = "stochastic")
scv.tl.velocity_graph(sample_one)
scv.pl.velocity_embedding(sample_one, basis = 'umap')

color = sample_one.uns['Cluster_colors']

当然，现在还有取巧的办法，那就是使用SCP或者seurat-wrapper等工具，直接RunSCVELO完事…..纯原生R环境，无需担心做映射…..可谓分析里的豪杰…..就是需要有一定的在R中构建py环境的小经验。

1
2
3
4
5
6
pancreas_sub <- RunSCVELO(
  srt = pancreas_sub, group_by = "SubCellType",
  linear_reduction = "PCA", nonlinear_reduction = "UMAP", return_seurat = TRUE
)
VelocityPlot(srt = pancreas_sub, reduction = "UMAP", group_by = "SubCellType")
VelocityPlot(srt = pancreas_sub, reduction = "UMAP", plot_type = "stream")

我只能说真的np……

目前看来RNA Velocity的分析，基本上就是cellranger+Velocyto+scanpy/Seurat+scVelo/dynamo/Deepvelo流程，或者Kallisto+Bustools+scanpy/Seurat+scVelo/dynamo/Deepvelo流程，基本已经固定了，这里强烈推荐全python流程，因为不需要匹配来匹配去，环境更加统一，分析完这一步后，再迁移到R环境画图都行。

要学的还有很多，下一步得继续学习多组学的分析了，特别是表观组学的分析，又是一个漫长的旅程。

更新：Velocyt.R和python的安装

在linux或者mac上这俩安装非常简单，一个直接安装一下boost库就能搞定，一个直接从git下载velocyto.py.zip然后用pip本地安装就能装好，但在windows上就是另外一回事了，特别是velocyto.R的安装，由于windows默认没有boost库，因此难点就在于怎么搞到编译好的boost库文件，这样velocyto在R安装的时候才能进行调用编译然后才能安装成功。

首先，感谢知乎的这份教程，这里提供了基础的解决思路，那就是搞到一份编译好的boost库，然后改动编译调用的pkg lib的位置即可。

第一时间想到的，当然是直接从linux上复制一份boost文件过来看行不行，很可惜，不行，boost库本身也需要编译，不同的平台编译出来的并不通用，后来直呼自己憨批，如果能通用那这问题早就不存在了。

但是对于不会搞windows编译的，从VS从头编译boost库八成要挂，还有什么办法呢？有没有什么地方能直接搞到boost.exe这种安装呢？

很幸运，有，就在SourceForge上，在这里：boost。这里面就有各个版本的boost binary的文件，我选了个比较老的版本，然后进去就找到了熟悉的exe，那还说什么，直接下载安装一条龙，然后进入安装的文件夹，就能直接发现里面有boost这个子文件夹，然后按照知乎教程来就行，下载velocyto.R，在里面新建inst/include文件夹，然后改动src/makevars文件夹里的语句，删除”PKG_LIBS”行，添加以下语句：

1
2
PKG_CPPFLAGS = -I../inst/include
PKG_LIBS= -lstdc++ $(LAPACK_LIBS) $(BLAS_LIBS)  $(SHLIB_OPENMP_CFLAGS) $(FLIBS)

然后去R里面，用本地安装即可毫无任何问题的安装上….

1
install.package('path/to/velocyto.R', repos = NULL, type = 'source')

完全没任何难度了，也不需要自己去编译boost，更不需要浪费一大堆硬盘空间安装一个VS编译器。非常nice的解决方案~