用Seurat做RNA Velocity

介绍：https://www.jianshu.com/p/63071b368be5
安装：注意：velocyto 需要 python 版本>=3.6.0

git clone https://github.com/velocyto-team/velocyto.py.git
cd velocyto.py
pip install -e .  # note the trailing dot

其他安装方式https://velocyto.org/velocyto.py/install/index.html#install

使用方法以10X为例

1.Running the CLI

Usage: velocyto run10x [OPTIONS] SAMPLEFOLDER GTFFILE

          Runs the velocity analysis for a Chromium 10X Sample

          10XSAMPLEFOLDER specifies the cellranger sample folder

          GTFFILE genome annotation file

        Options:
          -s, --metadatatable FILE        Table containing metadata of the various samples (csv fortmated rows are samples and cols are entries)
          -m, --mask FILE                 .gtf file containing intervals to mask
          -l, --logic TEXT                The logic to use for the filtering (default: Default)
          -M, --multimap                  Consider not unique mappings (not reccomended)
          -@, --samtools-threads INTEGER  The number of threads to use to sort the bam by cellID file using samtools
          --samtools-memory INTEGER       The number of MB used for every thread by samtools to sort the bam file
          -t, --dtype TEXT                The dtype of the loom file layers - if more than 6000 molecules/reads per gene per cell are expected set uint32 to avoid truncation (default run_10x: uint16)
          -d, --dump TEXT                 For debugging purposes only: it will dump a molecular mapping report to hdf5. --dump N, saves a cell every N cells. If p is prepended a more complete (but huge) pickle report is printed (default: 0)
          -v, --verbose                   Set the vebosity level: -v (only warinings) -vv (warinings and info) -vvv (warinings, info and debug)
          --help                          Show this message and exit.

需要准备的文件：
1.下载genome annotation file
Download a genome annotation (.gtf file) for example from GENCODE or Ensembl. If you use the cellranger pipeline, you should download the gtf that comes prepackaged with it here.
2.下载 expressed repeats annotation(可选）
You might want to mask expressed repetitive elements, since those count could constitute a confounding factor in the downstream analysis. To do so you would need to download an appropriate expressed repeat annotation (for example from UCSC genome browser and make sure to select GTF as output format).
我的数据是大鼠的，所以选择rn6

image.png

所以输入关键三个要素：
1.路径：outs的上级目录 (e.g.包含 outs, outs/analys and outs/filtered_gene_bc_matrices的总目录).
2.genes.gtf is the genome annotation file provided with the cellranger pipeline.
3.repeat_msk.gtf is the repeat masker file described in the Preparation section above.

例子：

#对bam进行sort
samtools sort -m 3G -t CB -O BAM -o yourpath_to_outs/outs/cellsorted_possorted_genome_bam.bam yourpath_to_outs/outs/possorted_genome_bam.bam

velocyto run10x -m yourpath/repeat_msk.gtf  yourpath_to_outs/NRVC yourpath/cellranger_rn6/genes/genes.gtf

如果出现samtools的报错：可以参考以下解决方案
https://github.com/velocyto-team/velocyto.py/issues/192#issuecomment-545349435

2.Estimating RNA velocity

This guide covers the analysis and assumes that you have produced a .loom file using the velocyto CLI (follow the guide above).

关于velocyto run10x输入的loom文件的介绍

A valid .loom file is simply an HDF5 file that contains specific groups representing the main matrix as well as row and column attributes. Because of this, .loom files can be created and read by any language that supports HDF5.
.loom files can be easily handled using the loompy package.

Get started with the analysis

At this point you are ready to start analyzing your .loom file. To get started read our analysis tutorial and have a look at the notebooks examples.

library(devtools)
install_github("velocyto-team/velocyto.R")
library(velocyto.R)
library(pagoda2)
ldat <- read.loom.matrices("velocyto/NRVC.loom")
#没有数据的话可以download 10X43_1.loom from the following URL: http://pklab.med.harvard.edu/velocyto/DG1/10X43_1.loom
#ldat <- read.loom.matrices("10X43_1.loom")

#Using spliced expression matrix as input to pagoda2.
emat <- ldat$spliced
hist(log10(colSums(emat)),col='wheat',xlab='cell size')

# filter 表达量低的样本
emat <- emat[,colSums(emat)>=1e3]
#将重复的行名去除
emat<-emat[!duplicated(rownames(emat)),]

#需要用到无重复的基因的基因表达矩阵
r <- Pagoda2$new(emat,modelType='plain',trim=10,log.scale=T)
r$adjustVariance(plot=T,do.par=T,gam.k=10)
#generate cell embedding and clustering, visualize
r$calculatePcaReduction(nPcs=100,n.odgenes=3e3,maxit=300)
r$makeKnnGraph(k=30,type='PCA',center=T,distance='cosine');
r$getKnnClusters(method=multilevel.community,type='PCA',name='multilevel')
r$getEmbedding(type='PCA',embeddingType='tSNE',perplexity=50,verbose=T)

#Plot embedding, labeling clusters and gene expression 
par(mfrow=c(1,2))
r$plotEmbedding(type='PCA',embeddingType='tSNE',show.legend=F,mark.clusters=T,min.group.size=10,shuffle.colors=F,mark.cluster.cex=1,alpha=0.3,main='cell clusters')
r$plotEmbedding(type='PCA',embeddingType='tSNE',colors=r$counts[,"Cdh5"],main='Cdh5')  

image.png image.png

Velocity estimation

#Prepare matrices and clustering data
emat <- ldat$spliced
nmat <- ldat$unspliced
emat <- emat[,rownames(r$counts)]; 
nmat <- nmat[,rownames(r$counts)]; # restrict to cells that passed p2 filter
# take cluster labels
cluster.label <- r$clusters$PCA[[1]]
cell.colors <- pagoda2:::fac2col(cluster.label)
# take embedding
emb <- r$embeddings$PCA$tSNE

#将cell-cell相关性转化为距离
cell.dist <- as.dist(1-armaCor(t(r$reductions$PCA)))
#过滤表达量低的基因
emat <- filter.genes.by.cluster.expression(emat,cluster.label,min.max.cluster.average = 0.5)
nmat <- filter.genes.by.cluster.expression(nmat,cluster.label,min.max.cluster.average = 0.05)
length(intersect(rownames(emat),rownames(emat)))

#Estimate RNA velocity (using gene-relative model with k=20 cell kNN pooling and using top/bottom 2% quantiles for gamma fit)
fit.quantile <- 0.02
rvel.cd <- gene.relative.velocity.estimates(emat,nmat,deltaT=1,kCells=20,cell.dist=cell.dist,fit.quantile=fit.quantile)

#Visualize velocity on the t-SNE embedding, using velocity vector fields
show.velocity.on.embedding.cor(emb,rvel.cd,n=300,scale='sqrt',cell.colors=ac(cell.colors,alpha=0.5),cex=0.8,arrow.scale=5,show.grid.flow=TRUE,min.grid.cell.mass=0.5,grid.n=40,arrow.lwd=1,do.par=F,cell.border.alpha = 0.1)

image.png

Visualize a fit for a particular gene (we reuse rvel.cd to save on calcualtions here):

gene <- "Myh7"
gene.relative.velocity.estimates(emat,nmat,deltaT=1,kCells = 20,kGenes=1,fit.quantile=fit.quantile,cell.emb=emb,cell.colors=cell.colors,cell.dist=cell.dist,show.gene=gene,old.fit=rvel.cd,do.par=T)

image.png

欢迎关注~

参考：
https://velocyto.org/velocyto.py/tutorial/cli.html#running-velocyto
http://velocyto.org
https://doi.org/10.1038/s41586-018-0414-6
http://pklab.med.harvard.edu/velocyto/notebooks/R/SCG71.nb.html
http://pklab.med.harvard.edu/velocyto/notebooks/R/DG1.nb.html