前言
前文跨物种RNA-seq标准化及差异表达介绍了跨物种标准化的方法，这次我们研究下其对应的R包SCBN的源代码。SCBN的用户手册地址为SCBN Tutorial
使用
假设我们对两个物种进行标准化，分别为 1 物种和 2 物种
sim_data
上图的 geneLength1 代表物种 1 每个基因的长度，count1 代表物种 1 对应基因的count数；geneLength2 代表物种 2 每个基因的长度，count2 代表物种 2 对应基因的count数
data(sim_data)
factor <- SCBN(orth_gene=sim_data, hkind=1:1000, a=0.05)
## hkind代表两个物种保守的基因在sim_data中的位置, 该例子中为第一个到第1000个是两个物种的保守基因
## α 代表p_val的阈值
运行完以后，factor 为：
我们可以简单的看一下SCBN的结构
SCBN <- function(orth_gene, hkind, a=0.05)
{
if (all(!is.na(orth_gene)) == FALSE) {
stop("The dataset of orthologous genes has NA values.", call.=TRUE)
} else if (all(hkind %in% (seq_len(nrow(orth_gene)))) == FALSE){
stop("The conserved genes are not included in dataset of
orthologous genes.", call.=TRUE)
}
scale <- MediancalcNorm(orth_gene, hkind)
scbn_val <- Iter_optimal(scale=scale, orth_gene=orth_gene, hkind=hkind, a=a)
list(median_val=scale, scbn_val=scbn_val)
}
不难发现它是由两个主要的函数构成，一个是 MediancalcNorm() 主要计算 median_val；另一个是 Iter_optimal() 主要计算 scbn_val
MediancalcNorm()
首先我们看一下 MediancalcNorm() 的功能：
MediancalcNorm <- function(orth_gene, hkind)
{
##对于物种 1 , 用物种 1 的count除以对应基因的长度, 再乘10^9 / (物种 1 所有基因count数之和
RPKMh <- orth_gene[, 2]/orth_gene[, 1]*(10^9/sum(orth_gene[, 2]))
##对于物种 2 , 用物种 2 的count除以对应基因的长度, 再乘10^9 / (物种 2 所有基因count数之和
RPKMm <- orth_gene[, 4]/orth_gene[, 3]*(10^9/sum(orth_gene[, 4]))
##挑选物种 1 和 2 的保守基因
cRPKMh <- RPKMh[hkind]
cRPKMm <- RPKMm[hkind]
## 标准化因子scale为两个物种保守基因矫正后count值的商
scale <- median(cRPKMh)/median(cRPKMm)
return(scale)
}
利用p_val判断两个物种的差异基因
计算p_val的方法页很简单，按照说明书上的方法：
data(sim_data)
orth_gene <- sim_data
## hkind 代表保守基因在sim_data的位置
hkind <- 1:1000
##scale 代表上面所计算的scbn_val
scale <- factor$scbn_val
x <- orth_gene[, 2]
y <- orth_gene[, 4]
lengthx <- orth_gene[, 1]
lengthy <- orth_gene[, 3]
n1 <- sum(x)
n2 <- sum(y)
p_value <- sageTestNew(x, y, lengthx, lengthy, n1, n2, scale)
通过上述命令就可以计算出每个基因对应的p_val来了
我们对应的看一下核心函数 sageTestNew() 的功能：
sageTestNew<-function (x, y, lengthx, lengthy, n1, n2, scale)
{
if (any(is.na(x)) || any(is.na(y)))
stop("missing values not allowed")
if (any(x < 0) || any(y < 0))
stop("x and y must be non-negative")
if (length(x) != length(y))
stop("x and y must have same length")
## x 代表物种 1 的count数
x <- round(x)
## y 代表物种 2 的count数
y <- round(y)
## n1(nn) 代表物种 1 总的count数
nn <- round(n1)
## n2(mm) 代表物种 2 总的count数
mm <- round(n2)
## size 代表物种 1 和物种 2 对应的每个基因(保守基因)count数之和(比方说保守的 gene A 在物种 1 的count数加上在物种 2 中的count数)
size <- x + y
pValue <- rep(1, length(size))
## 定义rate (scale为之前计算的scbn_val), lengthx为物种 1 每个保守基因的长度, lengthy为物种 2 对应每个保守基因的长度,
## 每个基因都有一个rate
rate <- scale*lengthx*nn/lengthy/mm
## 定义概率, 每一个基因都有个prob, 这个prob
prob <- rate/(1+rate)
#接下来需要做判断, 当size大于10000时, 采用卡方检验计算p_val
if (any(big <- size > 10000)) {
ibig <- (seq_along(x))[big]
## 每个 i 代表两个物种中每一个保守的基因
for (i in ibig)
pValue[i] <- chisq.test(matrix(c(x[i], y[i],
nn-x[i], mm-y[i]), 2, 2))$p.value
}
## 筛选出两个物种对应基因count之和大于0小于10000的基因
size0 <- size[size > 0 & !big]
## 筛选出两个物种对应基因count之和大于0小于10000的基因count
x1 <- x[size > 0 & !big]
## 筛选出两个物种对应基因count之和大于0小于10000的基因概率
prob1 <- prob[size > 0 & !big]
pValue1 <- rep(1, length(size0))
if (length(size0))
## 依次遍历每一个满足条件的基因
for(i in seq_len(length(size0))){
## 对每一个基因进行二项分布检验, 计算每次实验结果为物种 1 count的
## size0[i] (size0[i]是一个数值) 为第 i 个基因两个物种count数总和, 可以理解为二项分布试验总的量
## prob1[i] 为第 i 个基因二项分布的概率
## 0:size0[i] 可以理解为二项分布实验的次数为 0~size0[i] 次, 物种 1 的count数为 k
p <- dbinom(0:size0[i], prob=prob1[i], size=size0[i])
## 按l排序
o <- order(p)
## 做累积求和, 求极端情况的概率(p_val)
cumsump <- cumsum(p[o])[order(o)]
## x1[i] 为物种 1 第 i 个基因的count值
pValue1[i] <- cumsump[x1[i] + 1]
}
pValue[size > 0 & !big] <- pValue1
return(pValue)
}
对于size > 10000 的基因来说
这里作者采用的是卡方检验来完成p_val的计算，卡方检验的原理很简单，事实上就是列联表检验，其目的是检验列联表横纵遍历是否相关，如下图：
其中括号内是期望值，也就是说对于基因 1 来说物种 1 和 2 的count数应该是对半分的，也就是基因 1 在两个物种之间是没差异的，那么根据拟合优度公式计算p_val ：
如果p_val不显著，代表基因 1 与物种（物种1，2）这个变量之间是没有关系的，翻译为生物学现象即为基因 1 在两个物种之间不是差异表达的；
如果p_val显著，代表基因 1 与物种（物种1，2）这个变量之间是有关系的，翻译为生物学现象即为基因 1 在两个物种之间是差异表达的；
对于size < 10000 的基因来说
这一段比较有意思的是二项分布检验：
二项分布
我们知道离散型二项分布检验相当于是在做摸球游戏，比方说 gene A，在物种 1 的count数为10(相当于10个红球)，在物种 2 中的count数为20(相当于20个白球)，假设抽取到物种 1 的count的概率为 1/3，那么 0:size0[i](size0[i] =30) 相当于一共做了30次实验；P{ X = k }相当于做30次实验，摸到物种 1 的count(摸到红球)次数为 k(k取值0-10次) 的概率
那么对应的P即为选取物种 1 count为 k 时的概率，那么对于某物种，当k（次数）取值越小，而对应的概率越大时，代表该物种对应基因的count数越少；那么相应的另外一个物种的count数就会很多，从而产生差异表达，那么p_val定义为当选取物种 1 的count数(k)为10的时候，累积的概率P的和 ( ∑ P{ X < k })。当p_val显著，则代表该gene在物种 1 和 2 中差异表达
举个例子
##物种 1, count数为 1
dbinom(0:10, prob=0.1, size=10)
[1] 0.3486784401 0.3874204890 0.1937102445 0.0573956280 0.0111602610 0.0014880348 0.0001377810 0.0000087480 0.0000003645
[10] 0.0000000090 0.0000000001
##物种 2, count数为 9
dbinom(0:10, prob=0.9, size=10)
[1] 0.0000000001 0.0000000090 0.0000003645 0.0000087480 0.0001377810 0.0014880348 0.0111602610 0.0573956280 0.1937102445
[10] 0.3874204890 0.3486784401
>
我们可以看到两种检验结果成相反的结果，物种1的count较少，所以 k 取值小的时候，对应的概率大；而物种2的count较多，所以 k 取值大的时候，对应的概率大
不过这里有两个疑问：
- 为什么作者在做卡方检验的时候不对总的count做一个文库归一化处理
- 在二项分布定义 prob(rate) 的时候为什么考虑的是基因长度而不是表达量
Iter_optimal()
接着我们再看下Iter_optimal()的功能，这种标准化的模式是作者自己提出来的：
Iter_optimal <- function(scale, orth_gene, hkind, a) {
if (scale > 0.5) {
fW1 <- seq(scale-0.5, scale + 0.5, 0.1)
} else {
fW1 <- seq(0.02, scale +0.5, 0.1)
}
n1 <- length(fW1)
fdr1 <- rep(0, n1)
for(j in seq_len(n1)){
sMm1 <- sageTestNew(orth_gene[, 2], orth_gene[, 4], orth_gene[, 1],
orth_gene[, 3], sum(orth_gene[, 2]),
sum(orth_gene[, 4]), fW1[j])
fdr1[j] <- sum(sMm1[hkind] < a)/length(hkind)
}
fdr11 <- abs(fdr1-a)
counts1 <- length(which(fdr11 == min(fdr11)))
if (counts1 %% 2 == 1) {
fw2 <- median(fW1[which(fdr11 == min(fdr11))])
} else {
fw2 <- fW1[which(fdr11 == min(fdr11))][counts1/2+1]
}
if (fw2 > 0.5) {
fW3 <- seq(fw2-0.5, fw2+0.5, 0.05)
} else {
fW3 <- seq(0.02, fw2+0.5, 0.05)
}
n2 <- length(fW3)
fdr2 <- rep(0,n2)
for(j in seq_len(n2)){
sMm2 <- sageTestNew(orth_gene[, 2], orth_gene[, 4], orth_gene[,1],
orth_gene[, 3], sum(orth_gene[, 2]),
sum(orth_gene[, 4]), fW3[j])
fdr2[j] <- sum(sMm2[hkind] < a)/length(hkind)
}
fdr22 <- abs(fdr2-a)
counts2 <- length(which(fdr22 == min(fdr22)))
if (counts2%%2 == 1) {
fw4 <- median(fW3[which(fdr22 == min(fdr22))])
} else {
fw4 <- fW3[which(fdr22 == min(fdr22))][counts2/2+1]
}
if (fw4 > 0.25) {
fW5 <- seq(fw4-0.25,fw4+0.25,0.005)
} else {
fW5 <- seq(0.02, fw4+0.25, 0.005)
}
n3 <- length(fW5)
fdr3 <- rep(0,n3)
for (j in seq_len(n3)) {
sMm3 <- sageTestNew(orth_gene[, 2], orth_gene[, 4], orth_gene[, 1],
orth_gene[, 3], sum(orth_gene[, 2]),
sum(orth_gene[, 4]), fW5[j])
fdr3[j] <- sum(sMm3[hkind] < a)/length(hkind)
}
fdr33 <- abs(fdr3-a)
counts3 <- length(which(fdr33 == min(fdr33)))
if (counts3%%2 == 1) {
fw6 <- median(fW5[which(fdr33 == min(fdr33))])
} else {
fw6 <- fW5[which(fdr33 == min(fdr33))][counts3/2+1]
}
factor <- fw6
return(factor)
}
其核心思想就是根据显著性阈值（α）构造一列数组
fW1 <- seq(scale-0.5, scale + 0.5, 0.1)
## 或
fW1 <- seq(0.02, scale +0.5, 0.1)
##计算p_val
sMm1 <- sageTestNew(orth_gene[, 2], orth_gene[, 4], orth_gene[, 1],
orth_gene[, 3], sum(orth_gene[, 2]),
sum(orth_gene[, 4]), fW1[j])
## fW1[j] 为不同的显著性阈值
然后根据一系列方法去刷选，最后得到标准化因子 scbn_val
