介绍：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
所以输入关键三个要素：
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')
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)
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)
欢迎关注~
参考：
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
