参考链接:
https://github.com/velocyto-team/velocyto.R
http://velocyto.org/velocyto.py/index.html
http://pklab.med.harvard.edu/velocyto/notebooks/R/chromaffin2.nb.html
https://github.com/velocyto-team/velocyto.R/issues/16
https://www.cnblogs.com/raisok/p/12425258.html
目录
RNA速率：使用Seurat的结果做RNA velocity
导入Seurat以及loom对象
>library(Seurat)
## remotes::install_github('satijalab/seurat-wrappers')
>library(SeuratWrappers)
## 导入Seurat对象,之前分析的结果
>load("wang.rds")
>load('wang-loom.rds')
统一loom对象和Seurat的细胞名与基因名
> wt$spliced[1:3,1:3]
3 x 3 sparse Matrix of class "dgCMatrix"
WANG:AAAGTAGAGATGTTAGx WANG:AAACCTGTCAGCATGTx WANG:AAAGCAACATTTGCTTx
AT1G01020 . . .
AT1G01030 . . .
AT1G03993 . . .
> wang@meta.data[1:3,1:3]
orig.ident nCount_RNA nFeature_RNA
AAACCTGAGAATTCCC-1 zxz 3756 2158
AAACCTGAGGGCACTA-1 zxz 2774 1669
AAACCTGAGTAATCCC-1 zxz 2463 1290
> colnames(wt$spliced)<-gsub("x","-1",colnames(wt$spliced))
> colnames(wt$spliced)<-gsub("WANG:","",colnames(wt$spliced))
> colnames(wt$unspliced)<-colnames(wt$spliced)
> colnames(wt$ambiguous)<-colnames(wt$spliced)
计算velocity
提取spliced与unspliced文件，并提取原有的Seurat的UAMP图
## 由于Seurat的对象筛选了数据，所以两个文件细胞并不相同，以Seurat对象为准
> wt$spliced<-wt$spliced[,rownames(wang@meta.data)]
> wt$unspliced<-wt$unspliced[,rownames(wang@meta.data)]
> wt$ambiguous<-wt$ambiguous[,rownames(wang@meta.data)]
> sp <- wt$spliced
> unsp <- wt$unspliced
> WTumap <- wang@reductions$umap@cell.embeddings
## 估计细胞和细胞的距离
> cell.dist <- as.dist(1-armaCor(t(wang@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)
在UMAP聚类图上绘制RNA velocity
library(ggplot2)
pdf("cell_velocity.pdf",height=6,width=8)
gg <- UMAPPlot(wang)
ggplot_build(gg)$data
colors <- as.list(ggplot_build(gg)$data[[1]]$colour)
names(colors) <- rownames(WTumap)
p1 <- show.velocity.on.embedding.cor(WTumap,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()
这结果根本看不出啥啊，我选RC试试
id<-c('4','5','16','19')
Cell.sub <- subset(wang@meta.data,seurat_clusters==id)
scRNAsub <- subset(wang, cells=row.names(Cell.sub))
##后续的操作相同
倒是有个轨迹，但是并不明显
换成cluster12.14.19试试
这一次就有一个明显的轨迹了
总体来讲，感觉植物做RNA velocity的结果并没有动物的那么好，是自己分析的问题还是其他原因有待商榷。
转载请注明>>>周小钊的博客
