RNA velocity原理此前已经介绍过，参考单细胞测序的轨迹推断。
1. loom文件准备
由于RNA velocity分析的前提是要我们从单细胞RNA-seq的数据中区分出未成熟的mRNA(unspliced)和成熟的mRNA(spliced)，所以需要从fastq文件开始，与基因组进行比对后得到sam文件，从sam文件转成bam文件，再从bam文件中提取spliced，unspliced和ambiguous信息。得到.loom为后缀的文件。
（loom是scanpy常用的保存单细胞数据的格式）
2. Velocyto.R的使用练习（基于PAGODA2）
本练习基于：教程
练习数据下载：SCG71.loom
The example below starts with a loom file produced by velocyto.py, uses pagoda2 to obtain cell clusters/embedding, and then estimate/visualize velocity.
1.加载数据
library(velocyto.R)
input_loom <- "SCG71.loom"
ldat <- read.loom.matrices(input_loom)
View(ldat)
2. 准备pagoda2的输入数据
#使用剪切位点的表达量作为pagoda2的输入
emat <- ldat$spliced
hist(log10(colSums(emat)),col='wheat',xlab='cell size') #做直方图查看数据分布
emat <- emat[,colSums(emat)>=1e3]
#对数据进行过滤，滤掉剪切事件在1000以下的细胞，也就是上图中横轴小于3的细胞被滤掉了
#如果过滤了可以不进行
dim(emat)
# [1] 24421 2600
# 原先的6667个细胞只剩下2600个了
3. 使用Pagoda2 processing（标准化和细胞聚类）
PAGODA（pathway and gene set overdispersion analysis）是一个分析单细胞测序的方法，主要特点是在已知的重要信号通路基础上对细胞进行分类，以提高统计效力并揭示可能的功能性解释。pagoda2可以用来进行细胞聚类，生成细胞-细胞距离矩阵等（其他软件如Seurat2也可以进行同样的操作）。
由于RNA速率分析velocyto.R是基于pagoda的cluster和tsne，因此有必要学一下pagoda。PAGODA的Nature Methods原文链接：characterizing transcriptional heterogeneity through pathway and gene set overdispersion analysis
3.1 读入数据，构建pagoda对象并进行标准化
library(pagoda2) # 导入pagoda2包
# 构建Pagoda2对象
r <- Pagoda2$new(emat,modelType='plain',trim=10,log.scale=T)
3.2 对表达量差异很大的基因对下游分析所占比重进行调整
r$adjustVariance(plot=T,do.par=T,gam.k=10)
3.3 对细胞进行降维聚类和细胞嵌合分析tsne
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)
聚类结果可视化
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$depth,main='depth')
#不能绘制标题，也不能设置mfrow是为什么？
在tsne图的基础上观察某些特异基因的表达情况
gene <-"Ccr2"
r$plotEmbedding(type='PCA',embeddingType='tSNE',colors=r$counts[,gene],shuffle.colors=F,mark.cluster.cex=1,alpha=0.8,main=gene)
计算每个cluster的差异基因并可视化，比如画出cluster2中高表达基因的热图
r$getDifferentialGenes(type='PCA',verbose=T)
de=r$diffgenes$PCA$multilevel$`2`
r$plotGeneHeatmap(genes=rownames(de)[1:15],groups=r$clusters$PCA[[1]])
4. 速率估计
准备矩阵和聚类数据
emat <- ldat$spliced
nmat <- ldat$unspliced #忽略跨剪切位点的数据（数目太少）
#通过p2对细胞进行过滤
emat <- emat[,rownames(r$counts)]
nmat <- nmat[,rownames(r$counts)]
#对分类数据进行标记
cluster.label <- r$clusters$PCA$multilevel # take the cluster factor that was calculated by p2
library(sccore)
cell.colors <- fac2col(cluster.label)
# take embedding form p2
emb <- r$embeddings$PCA$tSNE
计算细胞间的距离（除了聚类和tSNE嵌合，在p2(pagoda2)过程中也可以得到一个细胞细胞距离矩阵，而这个矩阵比velocyto.R中通常使用的全转录组相关距离矩阵要好。
cell.dist <- as.dist(1-armaCor(t(r$reductions$PCA)))
基于最小平均表达量筛选基因（至少在一个簇中），输出产生的有效基因数
emat <- filter.genes.by.cluster.expression(emat,cluster.label,min.max.cluster.average = 0.2)
nmat <- filter.genes.by.cluster.expression(nmat,cluster.label,min.max.cluster.average = 0.05)
length(intersect(rownames(emat),rownames(nmat)))
计算RNA速率(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=25,cell.dist=cell.dist,fit.quantile=fit.quantile)
在tsne上可视化RNA速率结果
show.velocity.on.embedding.cor(emb,rvel.cd,n=200,scale='sqrt',cell.colors=ac(cell.colors,alpha=0.5),cex=0.8,arrow.scale=3,show.grid.flow=TRUE,min.grid.cell.mass=0.5,grid.n=40,arrow.lwd=1,do.par=F,cell.border.alpha = 0.1)
可视化特定的基因
gene <- "Camp"
gene.relative.velocity.estimates(emat,nmat,deltaT=1,kCells = 25,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)
增加k的数目
gene <- "Camp"
gene.relative.velocity.estimates(emat,nmat,deltaT=1,kCells = 100,kGenes=1,fit.quantile=fit.quantile,cell.emb=emb,cell.colors=cell.colors,cell.dist=cell.dist,show.gene=gene,do.par=T)
参考：http://pklab.med.harvard.edu/velocyto/notebooks/R/SCG71.nb.html
