前言

生信宝库之前的一篇推文：单细胞分析流程之RNA velocity（一）中，简单介绍了一下RNA velocity的安装和使用。有关为啥后面没有继续更新呢，是因为Immugent发现目前有很多RNA velocity相关的推文了。

但是Immugent平时在做RNA velocity分析时发现，它的结果并不是都完全符合自己的生物学预期，甚至有的时候和已知的结果相反。至于为什么会有这样不稳定的结果，Immugent推测在一定程度上和10x数据的dropout有关。因为使用smart-seq2的测序数据进行RNA velocity分析，结果就相对比较符合预期。

下面Immugent就通过代码实操展示，如何使用Supercell整合后的数据进行RNA velocity分析。

主要内容

为了方便大家对结果进行复现，Immugent还是使用SuperCell自带的pancreas数据作为输入。

library(SuperCell)
library(velocyto.R)

data("pancreas")
SC_pancreas <- SCimplify_for_velocity(
  pancreas$emat, 
  pancreas$nmat, 
  gamma = 10)

一行代码计算 RNA velocity。

Vel <- supercell_estimate_velocity(SC_pancreas$emat, SC_pancreas$nmat)
## Warning in supercell_estimate_velocity(SC_pancreas$emat, SC_pancreas$nmat): supercell_size was replaced with 1

## calculating cell knn ... done
## calculating convolved matrices ... done
## fitting gamma coefficients ... done. succesfful fit for 14807 genes
## filtered out 3459 out of 14807 genes due to low nmat-emat correlation
## filtered out 1069 out of 11348 genes due to low nmat-emat slope
## calculating RNA velocity shift ... done
## calculating extrapolated cell state ... done

下面就是出图了。

# Assign clusters to metacells based on the given single clustering
clusters <- supercell_assign(pancreas$meta$clusters, SC_pancreas$membership)

# Set up the color scheme
N.clusters <- length(unique(clusters))
pal <- setNames(colorRampPalette(RColorBrewer::brewer.pal(8, name= "Set1"))(N.clusters), 
                as.character(unique(clusters))) 
color <- setNames(pal[as.character(clusters)], names(clusters))

map_cluster_to_cluster_name <- c('Ngn3 low EP', 'Alpha', 'Delta', 'Beta', 'Pre-endocrine', 'Ngn3 high EP', 'Ductal', 'Epsilon')

size <- sqrt(1+table(SC_pancreas$membership))/sqrt(2)/2
set.seed(3)
tsne10 <- tSNE.velocity.plot(Vel, 
                             cell.colors = ac(color, alpha = 0.6), 
                             scale = "log", do.par = T,
                             delta.norm = FALSE, nPcs = 15, norm.nPcs = 15 * 10, perplexity = 50,
                             show.grid.flow = FALSE, grid.n = 20, grid.sd = NULL, min.grid.cell.mass = 1,
                             pcount = 1, verbose = TRUE, min.arrow.median.ratio = 1/10,
                             max.arrow.quantile = 0.9, arrow.scale = 1, arrow.lwd = 1,
                             xlab = "tSNE-x", ylab = "tSNE-y", size.norm = FALSE, cex = size
)
legend("topleft", legend=map_cluster_to_cluster_name,
       fill=pal,  cex=0.8)