RNA速率 | bam转loom+根据已有的Seurat对象可视化（避坑篇）

（一）velocyto生成loom文件（linux）

需要的是标准的基因组gtf注释文件和bam文件。

1.1 使用conda 构建一个velocyto软件分析环境

#创建velocyto环境
conda create -n velocyto 
#进入velocyto环境
conda activate velocyto 
#安装需要的包
conda install samtools
conda install numpy scipy cython numba matplotlib scikit-learn h5py click
pip install pysam
pip install velocyto

1.2 下载基因组gtf注释文件

#终端输入
curl -O https://cf.10xgenomics.com/supp/cell-exp/refdata-gex-GRCh38-2020-A.tar.gz
#解压得到指定的gtf文件
tar -zxvf /绝对路径/refdata-gex-GRCh38-2020-A.tar.gz refdata-gex-GRCh38-2020-A/genes/genes.gtf

1.3 bam文件转loom（耗时很长，耐心等待）

velocyto run /绝对路径/possorted_genome_bam.bam /绝对路径/refdata-gex-GRCh38-2020-A/genes/genes.gtf

参考：https://www.jianshu.com/p/17a9977ee2ca

（二）根据已有的Seurat对象可视化

2.1 loom和Seurat文件准备和预处理

注：Seurat对象是已经聚类分群过的。

进入服务器的R环境

#终端输入R，进入R环境
R
#载入需要的程辑包
library(velocyto.R)
library(Seurat)

注意！！！
velocyto.R包只能在linux环境使用，安装方法如下：

#在终端安装需要的包
conda install -c conda-forge openmp
conda install -c anaconda boost

#进入R环境
R
#需要dectools包安装github上的包
install.packages("devtools")
library(devtools)
#安装velocyto.R包，或者去官网下载通过devtools_install本地安装
install_github("velocyto-team/velocyto.R")

参考：https://blog.csdn.net/weixin_39568781/article/details/111249461
所有包就绪后

#读取loom文件
adata<-read.loom.matrices("/velocyto/possorted_genome_bam_SDVMC.loom")
#读取已经聚类分群过的seurat对象
load("/UMAP.Rdata")

查看一下Seurat对象

查看loom对象

统一loom对象和Seurat的细胞名与基因名

#注意：这里第二个参数是Seurat对象的行名字的尾巴（上面图片）
colnames(adata$spliced)<-paste0(colnames(adata$spliced),"_1")
#注意：这里第一个参数是loom的列名字的前面字符串（上面图片）
colnames(adata$spliced)<-gsub("possorted_genome_bam_SDVMC:","",colnames(adata$spliced))
colnames(adata$unspliced)<-colnames(adata$spliced)
colnames(adata$ambiguous)<-colnames(adata$spliced)

再次查看loom文件

2.2 计算RNA velocity

提取spliced与unspliced文件，并提取原有的Seurat的UAMP图

# 由于Seurat的对象筛选了数据，所以两个文件细胞并不相同，以Seurat对象为准
adata$spliced<-adata$spliced[,rownames(pbmc@meta.data)]
adata$unspliced<-adata$unspliced[,rownames(pbmc@meta.data)]
adata$ambiguous<-adata$ambiguous[,rownames(pbmc@meta.data)]
sp<-adata$spliced
unsp<-adata$unspliced
umap<-pbmc@reductions$umap@cell.embeddings
# 估计细胞的距离
cell.dist <- as.dist(1-armaCor(t(pbmc@reductions$umap@cell.embeddings)))
fit.quantile <- 0.02
rvel.cd <- gene.relative.velocity.estimates(sp,unsp,deltaT=2,kCells=10, cell.dist=cell.dist,fit.quantile=fit.quantile,n.cores=24)
#calculating cell knn ... done
#calculating convolved matrices ... done
#fitting gamma coefficients ... done. succesfful fit for 20405 genes
#filtered out 14436 out of 20405 genes due to low nmat-emat correlation
#filtered out 774 out of 5969 genes due to low nmat-emat slope
#calculating RNA velocity shift ... done
#calculating extrapolated cell state ... done

2.3 在UMAP聚类图上绘制RNA velocity

library(ggplot2)
pdf("/cell_velocity.pdf",height=6,width=8)
gg <- UMAPPlot(pbmc)
colors <- as.list(ggplot_build(gg)$data[[1]]$colour)
names(colors) <- rownames(umap)
p1 <- show.velocity.on.embedding.cor(umap,rvel.cd,n=30,scale='sqrt',cell.colors=ac(colors,alpha=0.5),cex=0.8,arrow.scale=2,show.grid.flow=T,min.grid.cell.mass=1.0,grid.n=50,arrow.lwd=1,do.par=F,cell.border.alpha =0.1,USE_OPENMP=1,n.cores=24,main="Cell Velocity")
dev.off()

在文件夹下查看pdf就是结果啦！
show.velocity.on.embedding.cor函数说明文档：https://rdrr.io/github/velocyto-team/velocyto.R/man/show.velocity.on.embedding.cor.html

参考：https://www.jianshu.com/p/ccb3403d49f8