RNA velocity分析练习（四）聚类以及velocity可视化

系列回顾：
1.RNA velocity分析练习（一）文件下载以及预处理
2.RNA velocity分析练习（二）软件安装
3.RNA velocity分析练习（三）生成loom文件

这一篇文章就正式开始进行RNA velocity的练习。本篇文章代码参考：
velocyto简单使用

（一）加载数据并过滤

> input_loom <- "MyTissue.loom"
> adata <- read.loom.matrices(input_loom)
# Use the spliced data as the input data
> spliced_adata <- adata$spliced
# seperate into E12.5 and E13.5
# 因为文章里只用了E12.5天的数据来进行RNA velocity的分析，所以这里我也把它们区分一下
> E12.5_data <- spliced_adata[,1:384]
> E13.5_data <- spliced_adata[,385:768]

看一下spliced mRNA count的分布（下面的分析都用E12.5数据）：

# spliced expression magnitude distribution across genes:
> hist(log10(rowSums(E12.5_data)+1),col='wheat',xlab='log10[ number of reads + 1]',main='number of reads per gene')

这里可以看到很多的基因的count值是0,后面我们要把这些基因去掉

过滤前看一下多少基因：

> dim(E12.5_data)
[1] 26610   384 #一共是26610个基因

把在所有细胞里都不表达的基因删掉：

> E12.5_data <- E12.5_data[as.numeric(rowSums(E12.5_data)) != 0,]
> dim(E12.5_data)
[1] 23803   384  #还剩23803个基因

保留至少在10个以上的细胞里表达的基因：

> fdata_E12.5 = E12.5_data[apply(E12.5_data,1,function(x) sum(x>1) > 10),]
> dim(fdata_E12.5)
[1] 13767   384 #这里还剩13767个基因

现在再看一下reads的分布：

> hist(log10(rowSums(fdata_E12.5)+1),col='wheat',xlab='log10[ number of reads + 1]',main='number of reads per gene')

（二）标准化

上面我们取的是spliced的数据（我们平时做表达分析的时候，实际上用的都是spliced），现在用它来进行数据的标准化，以便后面的细胞聚类。

> library(pagoda2)
#这里要先过滤一下基因名，如果不过滤的话，后面标准化步骤会报错，告诉你基因名有重复，我也不知道怎么会有重复，很奇怪
> fdata_E12.5 <- fdata_E12.5[duplicated(rownames(fdata_E12.5))=="FALSE",]
> dim(fdata_E12.5)
[1] 13765   384　#这里只少了２个基因，说明有２个基因重复了
> fdata_E12.5_nor <- Pagoda2$new(fdata_E12.5,modelType='plain',trim=10,log.scale=T)
384 cells, 13765 genes; normalizing ... using plain model winsorizing ... log scale ... done.
#可视化
> fdata_E12.5_nor$adjustVariance(plot=T,do.par=T,gam.k=10)
calculating variance fit ... using gam 2396 overdispersed genes ... 2396persisting ... done.

（三）细胞聚类

# clustering cells,embedding analysis
> fdata_E12.5_nor$calculatePcaReduction(nPcs=100,n.odgenes=3e3,maxit=300)
running PCA using 3000 OD genes .... done
> fdata_E12.5_nor$makeKnnGraph(k=30,type='PCA',center=T,distance='cosine')
creating space of type angular done
adding data ... done
building index ... done
querying ... done
> fdata_E12.5_nor$getKnnClusters(method=multilevel.community,type='PCA',name='multilevel')
> fdata_E12.5_nor$getEmbedding(type='PCA',embeddingType='tSNE',perplexity=50,verbose=T)
#make figure(clustering)
> par(mfrow=c(1,2))
> fdata_E12.5_nor$plotEmbedding(type='PCA',embeddingType='tSNE',show.legend=F,mark.clusters=T,min.group.size=10,shuffle.colors=F,mark.cluster.cex=1,alpha=0.7,main='cell clusters')
> fdata_E12.5_nor$plotEmbedding(type='PCA',embeddingType='tSNE',colors=fdata_E12.5_nor$depth,main='depth')

（四）RNA velocity分析

现在我们需要把spliced和unspliced的数据导出来，和我们上面得到的标准化后的数据做个交集。然后把RNA velocity嵌合进t-SNE的图里。

#prepare spliced data and unspliced data
> emat <- adata$spliced
> nmat <- adata$unspliced

# filtered cells according to fdata_E12.5_nor$counts
> emat <- emat[,rownames(fdata_E12.5_nor$counts)]
> nmat <- nmat[,rownames(fdata_E12.5_nor$counts)]

#lable the clustered data对分类数据进行标记
> cluster.label <- fdata_E12.5_nor$clusters$PCA$multilevel # take the cluster factor that was calculated by fdata_E12.5_nor
> cell.colors <- pagoda2:::fac2col(cluster.label)

# take embedding from fdata_E12.5_nor
> emb <- fdata_E12.5_nor$embeddings$PCA$tSNE

# Calculate the distance between cells
> cell.dist <- as.dist(1-armaCor(t(fdata_E12.5_nor$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)))
[1] 13016

# Calculate RNA velocity
> fit.quantile <- 0.02
> rvel.cd <- gene.relative.velocity.estimates(emat,nmat,deltaT=1,kCells=25,cell.dist=cell.dist,fit.quantile=fit.quantile)
calculating cell knn ... done
calculating convolved matrices ... done
fitting gamma coefficients ... done. succesfful fit for 13016 genes
filtered out 2115 out of 13016 genes due to low nmat-emat correlation
filtered out 1543 out of 10901 genes due to low nmat-emat slope
calculating RNA velocity shift ... done
calculating extrapolated cell state ... done

可视化：

> 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)
delta projections ... sqrt knn ... transition probs ... done
calculating arrows ... done
grid estimates ... grid.sd= 0.3592013  min.arrow.size= 0.007184026  max.grid.arrow.length= 0.02704895  done

可视化指定的基因：

> gene <- "Serpine2"
> 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,old.fit=rvel.cd,do.par=T)

> gene <- "Chga"
> 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,old.fit=rvel.cd,do.par=T)

（五）与文献原图进行比较

到这里就基本完成RNA velocity的分析了。然后我想把自己的图跟文献里的对比一下，于是就放在了PPT里做了一张组合图：

上面一半是文献原图，下面一半是我自己做的图。一比较真的好丑，难怪人家能发Nature。。。

可以看到，原文里的细胞分化轨迹是“蓝>红（黄）>绿”的一个顺序，我的是“绿色> 粉（红）>蓝>黄”，颜色不要紧，主要是要有一个线性的方向。所以这里文献里的SCP细胞在我的图里就是绿色群，Differentiation bridge是红色/粉色群，chromaffin就是我的黄色的群（或者蓝色？）

原文里展示的两个基因，Serpine2主要表达在SCP群里(b左)，unspliced mRNA主要出现在SCP群里(b中，右)，也就是说在SCP里表达量很高，因为有很多新转录出来的unspliced mRNA；Chga基因主要表达在chromaffin群里，也就是分化末端里(c左)，这个基因的unspliced mRNA主要在bridge群和一部分的chromaffin群里(c中，右)，说明这个基因在分化过渡阶段和分化终端表达量高。大体上我的结果和文章里的还是差不多的。说明分析过程没有太大的问题。