前言

前文跨物種RNA-seq標(biāo)準(zhǔn)化及差異表達(dá)介紹了跨物種標(biāo)準(zhǔn)化的方法，這次我們研究下其對(duì)應(yīng)的R包SCBN的源代碼。SCBN的用戶(hù)手冊(cè)地址為SCBN Tutorial

使用

假設(shè)我們對(duì)兩個(gè)物種進(jìn)行標(biāo)準(zhǔn)化，分別為 1 物種和 2 物種

sim_data

data(sim_data)
factor <- SCBN(orth_gene=sim_data, hkind=1:1000, a=0.05)

## hkind代表兩個(gè)物種保守的基因在sim_data中的位置, 該例子中為第一個(gè)到第1000個(gè)是兩個(gè)物種的保守基因
## α 代表p_val的閾值

運(yùn)行完以后，factor 為：

我們可以簡(jiǎn)單的看一下SCBN的結(jié)構(gòu)

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)
}

不難發(fā)現(xiàn)它是由兩個(gè)主要的函數(shù)構(gòu)成，一個(gè)是 MediancalcNorm() 主要計(jì)算 median_val；另一個(gè)是 Iter_optimal() 主要計(jì)算 scbn_val

MediancalcNorm()

首先我們看一下 MediancalcNorm() 的功能：

MediancalcNorm <- function(orth_gene, hkind)
{
##對(duì)于物種 1 , 用物種 1 的count除以對(duì)應(yīng)基因的長(zhǎng)度, 再乘10^9 / (物種 1 所有基因count數(shù)之和
  RPKMh <- orth_gene[, 2]/orth_gene[, 1]*(10^9/sum(orth_gene[, 2]))
##對(duì)于物種 2 , 用物種 2 的count除以對(duì)應(yīng)基因的長(zhǎng)度, 再乘10^9 / (物種 2 所有基因count數(shù)之和
  RPKMm <- orth_gene[, 4]/orth_gene[, 3]*(10^9/sum(orth_gene[, 4]))
##挑選物種 1 和 2 的保守基因
  cRPKMh <- RPKMh[hkind]
  cRPKMm <- RPKMm[hkind]
## 標(biāo)準(zhǔn)化因子scale為兩個(gè)物種保守基因矯正后count值的商
  scale <- median(cRPKMh)/median(cRPKMm)
  return(scale)
}

利用p_val判斷兩個(gè)物種的差異基因

計(jì)算p_val的方法頁(yè)很簡(jiǎn)單，按照說(shuō)明書(shū)上的方法：

data(sim_data)
orth_gene <- sim_data
## hkind 代表保守基因在sim_data的位置
hkind <- 1:1000

##scale 代表上面所計(jì)算的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)

通過(guò)上述命令就可以計(jì)算出每個(gè)基因?qū)?yīng)的p_val來(lái)了
我們對(duì)應(yīng)的看一下核心函數(shù) 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數(shù)
  x <- round(x)
  ## y 代表物種 2 的count數(shù)
  y <- round(y)
  ## n1(nn) 代表物種 1 總的count數(shù)
  nn <- round(n1)
  ## n2(mm) 代表物種 2 總的count數(shù)
  mm <- round(n2)

  ## size 代表物種 1 和物種 2 對(duì)應(yīng)的每個(gè)基因(保守基因)count數(shù)之和(比方說(shuō)保守的 gene A 在物種 1 的count數(shù)加上在物種 2 中的count數(shù))
  size <- x + y
  pValue <- rep(1, length(size))

  ## 定義rate (scale為之前計(jì)算的scbn_val), lengthx為物種 1 每個(gè)保守基因的長(zhǎng)度, lengthy為物種 2 對(duì)應(yīng)每個(gè)保守基因的長(zhǎng)度, 
  ## 每個(gè)基因都有一個(gè)rate
  rate <- scale*lengthx*nn/lengthy/mm
  ## 定義概率, 每一個(gè)基因都有個(gè)prob, 這個(gè)prob
  prob <- rate/(1+rate)

#接下來(lái)需要做判斷, 當(dāng)size大于10000時(shí), 采用卡方檢驗(yàn)計(jì)算p_val
  if (any(big <- size > 10000)) {
      ibig <- (seq_along(x))[big]
    ## 每個(gè) i 代表兩個(gè)物種中每一個(gè)保守的基因
      for (i in ibig)
           pValue[i] <- chisq.test(matrix(c(x[i], y[i],
                                           nn-x[i], mm-y[i]), 2, 2))$p.value
  }

## 篩選出兩個(gè)物種對(duì)應(yīng)基因count之和大于0小于10000的基因
  size0 <- size[size > 0 & !big]
## 篩選出兩個(gè)物種對(duì)應(yīng)基因count之和大于0小于10000的基因count
  x1 <- x[size > 0 & !big]
## 篩選出兩個(gè)物種對(duì)應(yīng)基因count之和大于0小于10000的基因概率
  prob1 <- prob[size > 0 & !big]
  pValue1 <- rep(1, length(size0))
  if (length(size0))
  ## 依次遍歷每一個(gè)滿(mǎn)足條件的基因
      for(i in seq_len(length(size0))){
        ## 對(duì)每一個(gè)基因進(jìn)行二項(xiàng)分布檢驗(yàn), 計(jì)算每次實(shí)驗(yàn)結(jié)果為物種 1 count的
        ## size0[i] (size0[i]是一個(gè)數(shù)值) 為第 i 個(gè)基因兩個(gè)物種count數(shù)總和, 可以理解為二項(xiàng)分布試驗(yàn)總的量
        ## prob1[i] 為第 i 個(gè)基因二項(xiàng)分布的概率
        ## 0:size0[i] 可以理解為二項(xiàng)分布實(shí)驗(yàn)的次數(shù)為 0~size0[i] 次, 物種 1 的count數(shù)為 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 個(gè)基因的count值
          pValue1[i] <- cumsump[x1[i] + 1]
    }

  pValue[size > 0 & !big] <- pValue1
  return(pValue)
}

對(duì)于size > 10000 的基因來(lái)說(shuō)
這里作者采用的是卡方檢驗(yàn)來(lái)完成p_val的計(jì)算，卡方檢驗(yàn)的原理很簡(jiǎn)單，事實(shí)上就是列聯(lián)表檢驗(yàn)，其目的是檢驗(yàn)列聯(lián)表橫縱遍歷是否相關(guān)，如下圖：

對(duì)于size < 10000 的基因來(lái)說(shuō)
這一段比較有意思的是二項(xiàng)分布檢驗(yàn)：

二項(xiàng)分布

舉個(gè)例子

##物種 1, count數(shù)為 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數(shù)為 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
> 

我們可以看到兩種檢驗(yàn)結(jié)果成相反的結(jié)果，物種1的count較少，所以 k 取值小的時(shí)候，對(duì)應(yīng)的概率大；而物種2的count較多，所以 k 取值大的時(shí)候，對(duì)應(yīng)的概率大

不過(guò)這里有兩個(gè)疑問(wèn)：

1. 為什么作者在做卡方檢驗(yàn)的時(shí)候不對(duì)總的count做一個(gè)文庫(kù)歸一化處理
2. 在二項(xiàng)分布定義 prob(rate) 的時(shí)候?yàn)槭裁纯紤]的是基因長(zhǎng)度而不是表達(dá)量

Iter_optimal()

接著我們?cè)倏聪翴ter_optimal()的功能，這種標(biāo)準(zhǔn)化的模式是作者自己提出來(lái)的：

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)
}

其核心思想就是根據(jù)顯著性閾值（α）構(gòu)造一列數(shù)組

 fW1 <- seq(scale-0.5, scale + 0.5, 0.1)
## 或
 fW1 <- seq(0.02, scale +0.5, 0.1)

##計(jì)算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] 為不同的顯著性閾值

然后根據(jù)一系列方法去刷選，最后得到標(biāo)準(zhǔn)化因子 scbn_val