单细胞转录组测序 (scRNA-seq)中拟时序分析是很常见的，是用来推断不同的细胞类型间的发育关系，有一类通过RNA速率来推断轨迹的方法，在此介绍一下。
RNA velocity
- RNA速率 (RNA velocity)，是基因表达状态的时间倒数，可以通过剪切和未剪切的mRNA数目来推断；
- RNA速率是一个高维向量，可以用来推断一定时间后单个细胞的未来转态；
- RNA速率可以帮助分析发育谱系和细胞动力学。
RNA velocity求解
1 对转录动力学进行数学建模
未剪切速率和剪切数随时间的变化分别定义如下：
其中，α, β和γ分别指转录、剪切和降解的反应率。
2 两个模型
- 假设剪切数随时间的变化是个常数，则：
- 假设未剪切的分子数是个常数，则：
3 求解数学模型
-
为了简化模型，我们假设：
- 细胞最终会达到一个稳定的状态，此时ds(t)/dt = 0
- RNA速率定义为偏离稳态的程度
- 所有基因的剪切反应率一致，即β=1
- 因此对于稳态细胞，可以通过γ=u(t)/s(t)来计算降解反应率γ，并由此计算RNA速率
计算γ参数
对于每个基因，通过最小二乘拟合以下线性方程来得到：
u和s矫正过后的未剪切和剪切数目。-
计算RNA速率
- 在model I中，RNA速率是个常数，且
- 在modle II中，RNA速率为
此处，u为未剪切数目，s0为初始剪切数目。
- 在model I中，RNA速率是个常数，且
4 改善模型表现
为了改善模型表现，作者引入了如下操作：
- 对剪切及未剪切的数目进行矫正
- 引入基因特异性的补偿
- 根据k最相似细胞对转录本数目进行pooling
- 根据相关性基因对转录本数目进行pooling
velocyto实践
文章[1]作者开发了可用于实际使用的软件velocyto.py和velocyto.R，本人使用过velocyto.R (v 0.6)，在此汇总一下。
步骤
- velocyto run根据bam比对文件生成unspliced与spliced矩阵如下：
gene S1_AAAACTAGACCG S1_AAACAATTCAGT S1_AAACCGCAGGGT
Pcmt1 0 2 0
Zc2hc1b 0 1 0
Aig1 0 0 1
- filter.genes.by.cluster.expression()函数根据基因表达均值过滤基因
filter.genes.by.cluster.expression <- function(emat,clusters,min.max.cluster.average=0.1) {
if(!any(colnames(emat) %in% names(clusters))) stop("provided clusters do not cover any of the emat cells!")
vc <- intersect(colnames(emat),names(clusters))
cl.emax <- apply(do.call(cbind,tapply(vc,as.factor(clusters[vc]),function(ii) Matrix::rowMeans(emat[,ii]))),1,max)
vi <- cl.emax>min.max.cluster.average;
emat[vi,]
}
- gene.relative.velocity.estimates()函数为主函数，主要进行以下步骤
- 矫正unspliced和spliced矩阵，通过knn近邻细胞和基因的表达进行矫正，并去除文库大小影响等
- 计算每个基因的降解速率γ，通过lm()函数
- get the transcriptional change matrix，即velocities
- 预测经过一个单位时间后的spliced count matrix
结果概览
- 原始的unspliced count matrix, 行为基因，列为细胞，数值为reads数
gene S1_AAAACTAGACCG S1_AAACAATTCAGT S1_AAACCGCAGGGT
Pcmt1 0 2 0
Zc2hc1b 0 1 0
Aig1 0 0 1
- 原始的spliced count matrix
gene S1_AAAACTAGACCG S1_AAACAATTCAGT S1_AAACCGCAGGGT
Pcmt1 2 0 0
Zc2hc1b 0 0 0
Aig1 0 0 0
- cellKNN matrix，k近邻则为1，否则为0
S1_AAAACTAGACCG S1_AAACAATTCAGT S1_AAACCGCAGGGT S1_AAACTTGAGCAA
S1_AAAACTAGACCG 1 . . .
S1_AAACAATTCAGT . 1 . .
S1_AAACCGCAGGGT . . 1 .
S1_AAACTTGAGCAA . . . 1
S1_AAATACGGAGAT . . . .
S1_AAATAACATGAT . . . .
S1_AAATACGGAGAT S1_AAATAACATGAT
S1_AAAACTAGACCG . .
S1_AAACAATTCAGT . .
S1_AAACCGCAGGGT . .
S1_AAACTTGAGCAA . .
S1_AAATACGGAGAT 1 .
S1_AAATAACATGAT . 1
- normalized spliced count matrix
S1_AAAACTAGACCG S1_AAACAATTCAGT S1_AAACCGCAGGGT S1_AAACTTGAGCAA
Pcmt1 1.44731535 1.20195086 0.7124199 0.36209745
Zc2hc1b 0.94905925 1.29440862 0.0000000 0.90524362
Aig1 0.00000000 0.00000000 0.0000000 0.00000000
1700020N01Rik 0.23726481 0.06934332 2.5647115 0.02263109
Rsph3 0.04745296 0.09245776 0.1424840 0.04526218
Arhgap35 0.02372648 0.04622888 0.0000000 0.02263109
S1_AAATACGGAGAT S1_AAATAACATGAT
Pcmt1 0.39392234 0.82674335
Zc2hc1b 0.84411930 0.61107117
Aig1 0.05627462 0.03594536
1700020N01Rik 0.28137310 1.07836089
Rsph3 0.16882386 0.25161754
Arhgap35 0.02813731 0.03594536
- 降解反应率γ
Pcmt1 Zc2hc1b Aig1 1700020N01Rik Rsph3
9.513799 2.984177 15.547821 1.663814 18.872405
Arhgap35
9.660245
- transcriptional change matrix: log2 observed unspliced count/ expected unspliced count ratio
S1_AAAACTAGACCG S1_AAACAATTCAGT S1_AAACCGCAGGGT S1_AAACTTGAGCAA
Pcmt1 -0.380501503 -0.35166118 -0.6124347 0.560501410
Zc2hc1b -0.638238515 -0.42208600 0.0000000 -0.562667628
Aig1 0.004363889 0.00000000 0.0000000 0.000000000
1700020N01Rik 0.273165122 -0.05620864 -1.7473054 0.217611106
Rsph3 0.076477001 -0.03211146 -0.1424840 0.001545896
Arhgap35 -0.023724968 -0.04622593 0.0235771 -0.022629647
S1_AAATACGGAGAT S1_AAATAACATGAT
Pcmt1 0.36471764 -0.23093740
Zc2hc1b -0.42503141 0.06064056
Aig1 -0.03492573 0.07786126
1700020N01Rik -0.22807674 -0.70481700
Rsph3 -0.08377039 -0.09708474
Arhgap35 -0.02813552 -0.03594307
- 预测的经过一个单位时间后的spliced count matrix
S1_AAAACTAGACCG S1_AAACAATTCAGT S1_AAACCGCAGGGT S1_AAACTTGAGCAA
Pcmt1 8.959464895 0.000000 0.3870690 0.561127391
Zc2hc1b 0.000000000 5.611424 0.0000000 0.000000000
Aig1 0.004361053 0.000000 0.0000000 0.000000000
1700020N01Rik 0.272987547 0.000000 5.1815394 0.217854139
Rsph3 0.076427287 0.000000 0.8444961 0.001547622
Arhgap35 0.000000000 0.000000 0.0229488 0.000000000
S1_AAATACGGAGAT S1_AAATAACATGAT
Pcmt1 0.3587747 0.00000000
Zc2hc1b 0.0000000 0.05853433
Aig1 0.0000000 0.07515691
1700020N01Rik 0.0000000 0.00000000
Rsph3 0.0000000 6.90097579
Arhgap35 0.0000000 0.00000000
-
速率可视化
image.png
后续模型[1-3]
velocyto.R运用的是steady-state model，面对复杂的数据，还有其他进行优化的模型如下：
- estimation based on gene structure, considering
- the number of exons
- total intronic length
- the number of predicted internal priming sites
- stochastic model
- to better capture the steady states
- treat transcription, splicing and degradation as probabilistic events
- dynamical model
- not simplify the model and
- solve the model in a likelihood-based expectation-maximization framework
- adapts RNA velocity to widely varying specifications such as non-stationary populations
参考文献
[1] RNA velocity of single cells
[2] http://velocyto.org/
[3] Generalizing RNA velocity to transient cell states through dynamical modeling.
