前言
Gene在转录为mRNA的过程中会经历splicing，RNA刚转录出来（此时称之为前体RNA）是没有经历过splicing的，而剪切过的RNA（此时称之为mRNA）从生成时间上要晚一些。即首先RNA转录出来为前体RNA，经过剪切后形成成熟的mRNA，因此在这个过程中存在时间差。由于每个细胞的RNA速率不同，因此可以从这个角度推测细胞的分化轨迹
其中，α代表转录速率，β 代表剪切速率，u 代表unspliced mRNA，γ 代表成熟mRNA的降解速率，s 代表spliced mRNA（成熟mRNA）。因此满足于下式：
上述式子分为两部分，du/dt 代表的是unspliced mRNA所能积累的速率，即转录出来的全部前体mRNA（α）等于剪切的部分 βu 加上积累的部分 du/dt ；而 ds/dt 代表的是spliced mRNA所积累的速率，剪切的部分 βu 等于 spliced mRNA积累的部分加上降解的部分 γs ；u 和 s 分别表示unspliced mRNA和spliced mRNA所积累的量。
当 t ≤ ts 时，转录开始；当 t ＞ ts 时，转录停止：
αon表示开启转录时的转录速率，αoff表示关闭转录时的转录速率=0。
模型解法
因此有两种模型参数的估计方法：
(1).Constant velocity assumption:
令
则当v < 0时，表现为该基因被下调；其中s0代表spliced mRNA的初值条件
(2).Constant unspliced molecules assumption:
其中s0代表spliced mRNA的初值条件，u0代表unspliced mRNA的初值条件，γ代表将解速率
velocyto.R pipeline
摘抄自：http://pklab.med.harvard.edu/velocyto/notebooks/R/chromaffin2.nb.html
现实中可以利用 velocyto.py 来注释 spliced 和 unspliced reads
library(velocyto.R)
# 加载数据
ldat <- readRDS(url("http://pklab.med.harvard.edu/velocyto/chromaffin/ldat.rds"))
cell.colors <- readRDS(url("http://pklab.med.harvard.edu/velocyto/chromaffin/cell.colors.rds"))
emb <- readRDS(url("http://pklab.med.harvard.edu/velocyto/chromaffin/embedding.rds"))
# 提取unspliced mRNA和spliced mRNA
# exonic read (spliced) expression matrix
emat <- ldat$spliced;
# intronic read (unspliced) expression matrix
nmat <- ldat$unspliced
# spanning read (intron+exon) expression matrix
smat <- ldat$spanning;
# filter expression matrices based on some minimum max-cluster averages
emat <- filter.genes.by.cluster.expression(emat,cell.colors,min.max.cluster.average = 5)
nmat <- filter.genes.by.cluster.expression(nmat,cell.colors,min.max.cluster.average = 1)
smat <- filter.genes.by.cluster.expression(smat,cell.colors,min.max.cluster.average = 0.5)
# look at the resulting gene set
length(intersect(rownames(emat),rownames(nmat)))
# 估计RNA速率
fit.quantile <- 0.05;
rvel.qf <- gene.relative.velocity.estimates(emat,nmat,deltaT=1,kCells = 5,fit.quantile = fit.quantile)
# 可视化
pca.velocity.plot(rvel.qf,nPcs=5,plot.cols=2,cell.colors=ac(cell.colors,alpha=0.7),cex=1.2,pcount=0.1,pc.multipliers=c(1,-1,-1,-1,-1))
R代码详解
项目地址：传送门
其中函数t.get.projected.delta()估计一段时间内s(t) 即spliced mRNA的变化采用的既是Constant unspliced molecules assumption模型
t.get.projected.delta <- function(em,nm,gamma,offset=rep(0,length(gamma)),delta=0.5) {
# adjust rownames
gn <- intersect(names(gamma),rownames(em));
if(is.null(names(offset))) { names(offset) <- names(gamma); }
em <- em[gn,]; nm <- nm[gn,]; gamma <- gamma[gn]; offset <- offset[gn];
# time effect constant
egt <- exp(-gamma*delta);
## 对nm（unspliced mRNA reads 做矫正）
y <- nm-offset; y[y<0] <- 0; # zero out entries with a negative n levels after offset adjustment
## Constant unspliced molecules assumption模型计算一段时间内s(t) 即spliced mRNA的变化
em*egt + (1-egt)*y/gamma - em
}
其中：
- em代表每个细胞，每个特定基因spliced mRNA的reads的初值 s0，即从测序测到的reads中得到，测序测到的reads就定义为初值；em可以理解为由 ldat$spliced 而得
- nm代表每个细胞，每个特定基因unspliced mRNA的reads的初值 u0；这里用 y = nm-offset 做了个矫正；nm可以理解为由 ldat$unspliced 而得
- egt代表
- gamma代表 γ，即RNA降解速率，估计gamma的方法在：传送门的第[4]节，估计gamma的代码如下：
# 作者定义steady.state.cells为 exonic read (spliced) expression cells steady.state.cells=colnames(emat) sfit <- data.frame(do.call(rbind,parallel::mclapply(sn(rownames(conv.emat.norm)),function(gn) { # gn 为gene knn聚类所筛选出来的每一个gene df <- data.frame(n=(conv.nmat.norm[gn,steady.state.cells]), e=(conv.emat.norm[gn,steady.state.cells]), s=conv.smat.norm[gn,steady.state.cells]) # 对每一个基因的emat [u(t)] 和 nmat [(s(t))] 建立线性关系，估计比例系数 # 这里的 n 与 s 是某一个基因在不同细胞中 unspliced 和 spliced 的表达量（列为细胞steady.state.cells，行为某一个基因） d <- lm(n~e,data=df,weights=pw) # note: re-estimating offset here return(c(o=as.numeric(coef(d)[1]), g=as.numeric(coef(d)[2]), r=cor(df$e,df$n,method='spearman'))) # omit genes for which sfit didn't work ko <- ko[rownames(ko) %in% rownames(sfit),]; vi <- sfit$r > min.nmat.smat.correlation ko <- ko[vi,] gamma <- ko$g
conv.nmat.norm和conv.emat.norm描述的是在不同细胞下某个基因unspliced和spliced的表达量，因此假设周期时间很短的话，每个基因在不同细胞中unspliced和spliced counts构成的坐标（如上图的绿点表示的是不同的细胞）大致在一条过原点的直线上；这里的conv.nmat.norm，conv.emat.norm来源于nmat和emat，代码如下：(值得注意的的是这里的基因和细胞都不是选择全部的基因和细胞，而是利用knn筛选了部分的基因和细胞）# 值得注意的的是这里的基因和细胞都不是选择全部的基因和细胞 # 而是利用knn筛选了部分的基因和细胞 if(!is.null(old.fit)) { cellKNN <- old.fit[['cellKNN']]} knn.maxl <- 1e2 if(kCells>1) { if(is.null(cellKNN)) { cat("calculating cell knn ... ") if(is.null(cell.dist)) { cellKNN <- balancedKNN(emat.log.norm,kCells,kCells*knn.maxl,n.threads=n.cores); } else { cellKNN <- balancedKNN(emat.log.norm,kCells,kCells*knn.maxl,n.threads=n.cores,dist=cell.dist); } diag(cellKNN) <- 1; resl$cellKNN <- cellKNN; cat("done\n") } rm(emat.log.norm); # smoothed matrices cat("calculating convolved matrices ... ") conv.emat <- emat %*% cellKNN[colnames(emat),colnames(emat)] conv.nmat <- nmat %*% cellKNN[colnames(nmat),colnames(nmat)] conv.emat.cs <- (emat.cs %*% cellKNN[colnames(emat),colnames(emat)])[1,] conv.nmat.cs <- (nmat.cs %*% cellKNN[colnames(nmat),colnames(nmat)])[1,] cat("done\n") } else { conv.emat <- emat; conv.nmat <- nmat; cellKNN <- NULL; conv.emat.cs <- emat.cs; conv.nmat.cs <- nmat.cs; }也就是说当估计gamma（γ）的值时，作者将不同细胞中对于某个基因实际测到的unspliced和spliced counts画在坐标系中（上图绿点代表不同细胞，横坐标代表 spliced counts，纵坐标代表 unspliced counts），可以看到这些细胞（绿点）大致在一条直线上，并用过原点的直线拟合，直线斜率为gamma（γ）
我们可以看到作者定义对单位时间（t=1）内s(t)的变化，即spliced mRNA的变化量定义为（本例中时间步长 t = 1）：
代码为：
# 单位时间内 s(t) 的变化
deltaE <- em*egt + (1-egt)*y/gamma - em
因此，作者定义未来对一段时间内s(t)的变化，即spliced mRNA的终量为：
# 定义基本的时间步长为 delta=1
delta=1
# 对 deltaE 进行矩阵化
deltae=as.matrix(deltaE)
target.mult=1e3
model.mult=1e3
# 初始化
rz <- matrix(0,nrow=nrow(em),ncol=ncol(em))
colnames(rz) <- colnames(em)
rownames(rz) <- rownames(em)
# 获得 deltae 的数据
gn <- intersect(rownames(deltae),rownames(rz))
## 对一段时间内s(t)的变化量，即spliced mRNA的变化量：deltaE进行筛选和矫正
rz[match(gn,rownames(rz)),colnames(deltae)] <- deltae[gn,];
rz <- t(t(rz)*cellSize)
## emm为spliced mRNA reads的初始值
emm <- t(t(em)*cellSize)
## spliced mRNA reads的终量为spliced mRNA reads的初始值加上 rz*delta
## delta=1 代表时间步长, rz 代表筛选后的 deltaE
emn <- emm + rz*delta
## emn为spliced mRNA reads的终值
emn[emn<0] <- 0;
newCellSize <- (cellSize+Matrix::colSums(emn-emm)/mult)
emn <- t(t(emn)/newCellSize)
1.定义单位时间spliced mRNA reads的变化量 deltaE 为
deltaE <- em*egt + (1-egt)*y/gamma - em
2.本例中时间步长为 delta=1
- rz 为 deltaE 筛选和矫正后的矩阵，表示为一段时间内s(t)的变化量，即spliced mRNA的变化量
- 一个单位时间后，spliced mRNA reads初始量和spliced mRNA reads终值之间的关系：
emn <- emm + rz*delta，emn 为spliced mRNA reads终量；emm 为spliced mRNA reads初始值， rz 代表筛选后的 deltaE
因此在降维画图的时候，箭头是如何表示出来的呢？以PCA降维为例：
从代码上来看：
# 其中vel$current代表spliced mRNA reads的初始值 emm
x0 <- vel$current
# 其中vel$projected代表spliced mRNA reads的终值 emn
x1 <- vel$projected
# transform
x0.log <- log2(x0+pcount)
x1.log <- log2(x1+pcount)
cent <- rowMeans(x0.log)
# 对spliced mRNA reads的初始值进行PCA降维构建低维空间的cell坐标，取PC1和PC2分别作为横纵坐标值，epc
## x0.log-cent 代表中心化数据
epc <- pcaMethods::pca(t(x0.log-cent),center=F,nPcs=ifelse(is.na(norm.nPcs),nPcs,norm.nPcs))
# 继续进行变化，epc@scores 代表经过变化后低维空间的spliced mRNA reads的初始值
## epc@loadings 代表 PCA 的载荷
epc@loadings <- t(t(epc@loadings)*pc.multipliers)
epc@scores <- scale(epc@completeObs,scale=F,center=T) %*% epc@loadings;
# 继续进行变化，x1.scores代表经过变化后低维空间的spliced mRNA reads的终值
## 低维空间 spliced mRNA reads的终值 = 中心化后的高维空间spliced mRNA reads的终值 × spliced mRNA reads的初始值PCA降维后的载荷值
x1.scores <- t(x1.log - cent) %*% epc@loadings
# 计算spliced mRNA reads的终值和spliced mRNA reads的初始值之差，定义为spliced mRNA reads的变化量
delta.pcs <- as.matrix(x1.scores-epc@scores)
× spliced mRNA reads的初始值PCA降维后的载荷值的目的是为了降维，因为
epc@loadings是一个2列矩阵，因此可以将维数降到二维
其中vel$current代表spliced mRNA reads的初始值 emm； 其中vel$projected代表spliced mRNA rea
那么箭头指向的定义如下：
## pos代表每个cell spliced mRNA reads的初始值的二维坐标，i 为每一个细胞
pos <- epc@scores[,c((i-1)+1,(i-1)+2)]
## ppos代表每个cell spliced mRNA reads的终值的二维坐标，i 为每一个细胞
ppos <- pos+delta.pcs[,c((i-1)+1,(i-1)+2)]
## 将每个细胞的二维坐标画在图上
plot(pos,bg=cell.colors[rownames(pos)],pch=21,
col=ac(1,alpha=0.3),lwd=0.5,xlab=paste("PC",(i-1)+1),
ylab=paste("PC",(i-1)+2),axes=T,
main=paste('PC',(i-1)+1,' vs. PC',(i-1)+2,sep=''));
box();
## 定义箭头指向
ars <- data.frame(pos[,1],pos[,2],ppos[,1],ppos[,2])
colnames(ars) <- c('x0','y0','x1','y1')
## 定义箭头的变化量
arsd <- data.frame(xd=ars$x1-ars$x0,yd=ars$y1-ars$y0)
rownames(ars) <- rownames(arsd) <- rownames(pos)
## 对箭头指向进行矫正
rx <- range(c(range(ars$x0),range(ars$x1)))
ry <- range(c(range(ars$y0),range(ars$y1)))
gx <- seq(rx[1],rx[2],length.out=grid.n)
gy <- seq(ry[1],ry[2],length.out=grid.n)
## 对箭头指向进行矫正
garrows <- do.call(rbind,lapply(gx,function(x) {
cd <- sqrt(outer(pos[,2],-gy,'+')^2 + (x-pos[,1])^2)
cw <- dnorm(cd,sd=grid.sd)
gw <- Matrix::colSums(cw)
cws <- pmax(1,Matrix::colSums(cw));
gxd <- Matrix::colSums(cw*arsd$xd)/cws
gyd <- Matrix::colSums(cw*arsd$yd)/cws
al <- sqrt(gxd^2+gyd^2);
vg <- gw>=min.grid.cell.mass & al>=min.arrow.size
cbind(rep(x,sum(vg)),gy[vg],x+gxd[vg],gy[vg]+gyd[vg])
}))
## x0：代表spliced mRNA reads的初始值的横坐标；
## x1：代表spliced mRNA reads的终值的横坐标；
## y0：代表spliced mRNA reads的初始值的纵坐标；
## y1：代表spliced mRNA reads的终值的纵坐标；
colnames(garrows) <- c('x0','y0','x1','y1')
## 箭头的画图为
# plot
alen <- pmin(max.grid.arrow.length,sqrt( ((garrows[,3]-garrows[,1]) * par('pin')[1] / diff(par('usr')[c(1,2)]) )^2 + ((garrows[,4]-garrows[,2])*par('pin')[2] / diff(par('usr')[c(3,4)]) )^2))
## arrows画箭头指向，i 为每一个细胞
suppressWarnings(lapply(1:nrow(garrows),function(i) arrows(garrows[i,1],garrows[i,2],garrows[i,3],garrows[i,4],length=alen[i],lwd=arrow.lwd)))
其中对于garrows里面的列元素来说：
x0：代表spliced mRNA reads的初始值经过PCA降维后的横坐标；
x1：代表spliced mRNA reads的终值经过PCA降维后的横坐标；
y0：代表spliced mRNA reads的初始值经过PCA降维后的纵坐标；
y1：代表spliced mRNA reads的终值经过PCA降维后的纵坐标；
那么对于每一个细胞来说，在确定箭头指向的时候，并非是选取所有基因来做特征选取，而是利用knn对基因进行一定的筛选，利用表达量特征比较相似的基因来确定箭头的指向
