前言
Gene在轉(zhuǎn)錄為mRNA的過(guò)程中會(huì)經(jīng)歷splicing，RNA剛轉(zhuǎn)錄出來(lái)（此時(shí)稱之為前體RNA）是沒(méi)有經(jīng)歷過(guò)splicing的移斩，而剪切過(guò)的RNA（此時(shí)稱之為mRNA）從生成時(shí)間上要晚一些竹观。即首先RNA轉(zhuǎn)錄出來(lái)為前體RNA藻糖，經(jīng)過(guò)剪切后形成成熟的mRNA江咳，因此在這個(gè)過(guò)程中存在時(shí)間差峭火。由于每個(gè)細(xì)胞的RNA速率不同爱致，因此可以從這個(gè)角度推測(cè)細(xì)胞的分化軌跡
其中烤送，α代表轉(zhuǎn)錄速率，β 代表剪切速率糠悯，u 代表unspliced mRNA帮坚，γ 代表成熟mRNA的降解速率，s 代表spliced mRNA（成熟mRNA）互艾。因此滿足于下式：
上述式子分為兩部分试和，du/dt 代表的是unspliced mRNA所能積累的速率，即轉(zhuǎn)錄出來(lái)的全部前體mRNA（α）等于剪切的部分 βu 加上積累的部分 du/dt 纫普；而 ds/dt 代表的是spliced mRNA所積累的速率阅悍，剪切的部分 βu 等于 spliced mRNA積累的部分加上降解的部分 γs ；u 和 s 分別表示unspliced mRNA和spliced mRNA所積累的量。
當(dāng) t ≤ ts 時(shí)节视，轉(zhuǎn)錄開(kāi)始晦墙；當(dāng) t ＞ ts 時(shí)，轉(zhuǎn)錄停止：
αon表示開(kāi)啟轉(zhuǎn)錄時(shí)的轉(zhuǎn)錄速率肴茄，αoff表示關(guān)閉轉(zhuǎn)錄時(shí)的轉(zhuǎn)錄速率=0晌畅。
模型解法
因此有兩種模型參數(shù)的估計(jì)方法：
(1).Constant velocity assumption:
令
則當(dāng)v < 0時(shí)，表現(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)實(shí)中可以利用 velocyto.py 來(lái)注釋 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)))
# 估計(jì)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代碼詳解
項(xiàng)目地址：傳送門
其中函數(shù)t.get.projected.delta()估計(jì)一段時(shí)間內(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);
## 對(duì)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模型計(jì)算一段時(shí)間內(nèi)s(t) 即spliced mRNA的變化
em*egt + (1-egt)*y/gamma - em
}
其中：
- em代表每個(gè)細(xì)胞拦坠，每個(gè)特定基因spliced mRNA的reads的初值 s0连躏，即從測(cè)序測(cè)到的reads中得到，測(cè)序測(cè)到的reads就定義為初值贞滨；em可以理解為由 ldat$spliced 而得
- nm代表每個(gè)細(xì)胞入热，每個(gè)特定基因unspliced mRNA的reads的初值 u0；這里用 y = nm-offset 做了個(gè)矯正晓铆；nm可以理解為由 ldat$unspliced 而得
- egt代表
- gamma代表 γ勺良，即RNA降解速率，估計(jì)gamma的方法在：傳送門的第[4]節(jié)骄噪，估計(jì)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聚類所篩選出來(lái)的每一個(gè)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]) # 對(duì)每一個(gè)基因的emat [u(t)] 和 nmat [(s(t))] 建立線性關(guān)系尚困，估計(jì)比例系數(shù) # 這里的 n 與 s 是某一個(gè)基因在不同細(xì)胞中 unspliced 和 spliced 的表達(dá)量（列為細(xì)胞steady.state.cells，行為某一個(gè)基因） 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描述的是在不同細(xì)胞下某個(gè)基因unspliced和spliced的表達(dá)量链蕊，因此假設(shè)周期時(shí)間很短的話事甜，每個(gè)基因在不同細(xì)胞中unspliced和spliced counts構(gòu)成的坐標(biāo)（如上圖的綠點(diǎn)表示的是不同的細(xì)胞）大致在一條過(guò)原點(diǎn)的直線上；這里的conv.nmat.norm滔韵，conv.emat.norm來(lái)源于nmat和emat逻谦，代碼如下：(值得注意的的是這里的基因和細(xì)胞都不是選擇全部的基因和細(xì)胞，而是利用knn篩選了部分的基因和細(xì)胞）# 值得注意的的是這里的基因和細(xì)胞都不是選擇全部的基因和細(xì)胞 # 而是利用knn篩選了部分的基因和細(xì)胞 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; }也就是說(shuō)當(dāng)估計(jì)gamma（γ）的值時(shí)陪蜻，作者將不同細(xì)胞中對(duì)于某個(gè)基因?qū)嶋H測(cè)到的unspliced和spliced counts畫在坐標(biāo)系中（上圖綠點(diǎn)代表不同細(xì)胞邦马，橫坐標(biāo)代表 spliced counts，縱坐標(biāo)代表 unspliced counts）囱皿，可以看到這些細(xì)胞（綠點(diǎn)）大致在一條直線上勇婴，并用過(guò)原點(diǎn)的直線擬合忱嘹，直線斜率為gamma（γ）
我們可以看到作者定義對(duì)單位時(shí)間（t=1）內(nèi)s(t)的變化嘱腥，即spliced mRNA的變化量定義為（本例中時(shí)間步長(zhǎng) t = 1）：
代碼為：
# 單位時(shí)間內(nèi) s(t) 的變化
deltaE <- em*egt + (1-egt)*y/gamma - em
因此，作者定義未來(lái)對(duì)一段時(shí)間內(nèi)s(t)的變化拘悦，即spliced mRNA的終量為：
# 定義基本的時(shí)間步長(zhǎng)為 delta=1
delta=1
# 對(duì) deltaE 進(jìn)行矩陣化
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))
## 對(duì)一段時(shí)間內(nèi)s(t)的變化量齿兔，即spliced mRNA的變化量：deltaE進(jìn)行篩選和矯正
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 代表時(shí)間步長(zhǎng), 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.定義單位時(shí)間spliced mRNA reads的變化量 deltaE 為
deltaE <- em*egt + (1-egt)*y/gamma - em
2.本例中時(shí)間步長(zhǎng)為 delta=1
- rz 為 deltaE 篩選和矯正后的矩陣，表示為一段時(shí)間內(nèi)s(t)的變化量，即spliced mRNA的變化量
- 一個(gè)單位時(shí)間后分苇，spliced mRNA reads初始量和spliced mRNA reads終值之間的關(guān)系：
emn <- emm + rz*delta添诉，emn 為spliced mRNA reads終量；emm 為spliced mRNA reads初始值医寿， rz 代表篩選后的 deltaE
因此在降維畫圖的時(shí)候栏赴，箭頭是如何表示出來(lái)的呢？以PCA降維為例：
從代碼上來(lái)看：
# 其中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)
# 對(duì)spliced mRNA reads的初始值進(jìn)行PCA降維構(gòu)建低維空間的cell坐標(biāo)靖秩，取PC1和PC2分別作為橫縱坐標(biāo)值须眷，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ù)進(jìn)行變化，epc@scores 代表經(jīng)過(guò)變化后低維空間的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ù)進(jìn)行變化沟突，x1.scores代表經(jīng)過(guò)變化后低維空間的spliced mRNA reads的終值
## 低維空間 spliced mRNA reads的終值 = 中心化后的高維空間spliced mRNA reads的終值 × spliced mRNA reads的初始值PCA降維后的載荷值
x1.scores <- t(x1.log - cent) %*% epc@loadings
# 計(jì)算spliced mRNA reads的終值和spliced mRNA reads的初始值之差花颗，定義為spliced mRNA reads的變化量
delta.pcs <- as.matrix(x1.scores-epc@scores)
× spliced mRNA reads的初始值PCA降維后的載荷值的目的是為了降維，因?yàn)?code>epc@loadings是一個(gè)2列矩陣惠拭，因此可以將維數(shù)降到二維
其中vel$current代表spliced mRNA reads的初始值 emm扩劝； 其中vel$projected代表spliced mRNA rea
那么箭頭指向的定義如下：
## pos代表每個(gè)cell spliced mRNA reads的初始值的二維坐標(biāo)，i 為每一個(gè)細(xì)胞
pos <- epc@scores[,c((i-1)+1,(i-1)+2)]
## ppos代表每個(gè)cell spliced mRNA reads的終值的二維坐標(biāo)职辅，i 為每一個(gè)細(xì)胞
ppos <- pos+delta.pcs[,c((i-1)+1,(i-1)+2)]
## 將每個(gè)細(xì)胞的二維坐標(biāo)畫在圖上
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)
## 對(duì)箭頭指向進(jìn)行矯正
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)
## 對(duì)箭頭指向進(jì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的初始值的橫坐標(biāo)棒呛；
## x1：代表spliced mRNA reads的終值的橫坐標(biāo)；
## y0：代表spliced mRNA reads的初始值的縱坐標(biāo)域携；
## y1：代表spliced mRNA reads的終值的縱坐標(biāo)条霜；
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 為每一個(gè)細(xì)胞
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)))
其中對(duì)于garrows里面的列元素來(lái)說(shuō)：
x0：代表spliced mRNA reads的初始值經(jīng)過(guò)PCA降維后的橫坐標(biāo)涵亏；
x1：代表spliced mRNA reads的終值經(jīng)過(guò)PCA降維后的橫坐標(biāo)宰睡；
y0：代表spliced mRNA reads的初始值經(jīng)過(guò)PCA降維后的縱坐標(biāo)；
y1：代表spliced mRNA reads的終值經(jīng)過(guò)PCA降維后的縱坐標(biāo)气筋；
那么對(duì)于每一個(gè)細(xì)胞來(lái)說(shuō)拆内，在確定箭頭指向的時(shí)候，并非是選取所有基因來(lái)做特征選取宠默，而是利用knn對(duì)基因進(jìn)行一定的篩選麸恍，利用表達(dá)量特征比較相似的基因來(lái)確定箭頭的指向
