前言
Gene在轉(zhuǎn)錄為mRNA的過程中會經(jīng)歷splicing，RNA剛轉(zhuǎn)錄出來（此時稱之為前體RNA）是沒有經(jīng)歷過splicing的，而剪切過的RNA（此時稱之為mRNA）從生成時間上要晚一些。即首先RNA轉(zhuǎn)錄出來為前體RNA，經(jīng)過剪切后形成成熟的mRNA，因此在這個過程中存在時間差。由于每個細胞的RNA速率不同，因此可以從這個角度推測細胞的分化軌跡
其中，α代表轉(zhuǎn)錄速率，β 代表剪切速率，u 代表unspliced mRNA，γ 代表成熟mRNA的降解速率，s 代表spliced mRNA（成熟mRNA）。因此滿足于下式：
上述式子分為兩部分，du/dt 代表的是unspliced mRNA所能積累的速率，即轉(zhuǎn)錄出來的全部前體mRNA（α）等于剪切的部分 βu 加上積累的部分 du/dt ；而 ds/dt 代表的是spliced mRNA所積累的速率，剪切的部分 βu 等于 spliced mRNA積累的部分加上降解的部分 γs ；u 和 s 分別表示unspliced mRNA和spliced mRNA所積累的量。
當 t ≤ ts 時，轉(zhuǎn)錄開始；當 t ＞ ts 時，轉(zhuǎn)錄停止：
αon表示開啟轉(zhuǎn)錄時的轉(zhuǎn)錄速率，αoff表示關(guān)閉轉(zhuǎn)錄時的轉(zhuǎn)錄速率=0。
模型解法
因此有兩種模型參數(shù)的估計方法：
(1).Constant velocity assumption:
令
則當v < 0時，表現(xiàn)為該基因被下調(diào)；其中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
現(xiàn)實中可以利用 velocyto.py 來注釋 spliced 和 unspliced reads
library(velocyto.R)
# 加載數(shù)據(jù)
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代碼詳解
項目地址：傳送門
其中函數(shù)t.get.projected.delta()估計一段時間內(nèi)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模型計算一段時間內(nèi)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]節(jié)，估計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))] 建立線性關(guān)系，估計比例系數(shù) # 這里的 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構(gòu)成的坐標（如上圖的綠點表示的是不同的細胞）大致在一條過原點的直線上；這里的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（γ）的值時，作者將不同細胞中對于某個基因?qū)嶋H測到的unspliced和spliced counts畫在坐標系中（上圖綠點代表不同細胞，橫坐標代表 spliced counts，縱坐標代表 unspliced counts），可以看到這些細胞（綠點）大致在一條直線上，并用過原點的直線擬合，直線斜率為gamma（γ）
我們可以看到作者定義對單位時間（t=1）內(nèi)s(t)的變化，即spliced mRNA的變化量定義為（本例中時間步長 t = 1）：
代碼為：
# 單位時間內(nèi) s(t) 的變化
deltaE <- em*egt + (1-egt)*y/gamma - em
因此，作者定義未來對一段時間內(nèi)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 的數(shù)據(jù)
gn <- intersect(rownames(deltae),rownames(rz))
## 對一段時間內(nèi)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 篩選和矯正后的矩陣，表示為一段時間內(nèi)s(t)的變化量，即spliced mRNA的變化量
- 一個單位時間后，spliced mRNA reads初始量和spliced mRNA reads終值之間的關(guān)系：
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降維構(gòu)建低維空間的cell坐標，取PC1和PC2分別作為橫縱坐標值，epc
## x0.log-cent 代表中心化數(shù)據(jù)
epc <- pcaMethods::pca(t(x0.log-cent),center=F,nPcs=ifelse(is.na(norm.nPcs),nPcs,norm.nPcs))
# 繼續(xù)進行變化，epc@scores 代表經(jīng)過變化后低維空間的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;
# 繼續(xù)進行變化，x1.scores代表經(jīng)過變化后低維空間的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列矩陣，因此可以將維數(shù)降到二維
其中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的初始值經(jīng)過PCA降維后的橫坐標；
x1：代表spliced mRNA reads的終值經(jīng)過PCA降維后的橫坐標；
y0：代表spliced mRNA reads的初始值經(jīng)過PCA降維后的縱坐標；
y1：代表spliced mRNA reads的終值經(jīng)過PCA降維后的縱坐標；
那么對于每一個細胞來說，在確定箭頭指向的時候，并非是選取所有基因來做特征選取，而是利用knn對基因進行一定的篩選，利用表達量特征比較相似的基因來確定箭頭的指向
