關(guān)于sva包

sva的用戶說明：http://www.bioconductor.org/packages/release/bioc/vignettes/sva/inst/doc/sva.pdf

其中關(guān)于去除批次效應(yīng)的函數(shù)有Combat和Combat_seq兩個(gè)歉提，而Combat主要面對的是帶有小數(shù)的數(shù)據(jù)（比方說芯片數(shù)據(jù)）区转，基于的是貝葉斯原理；而Combat_seq主打RNA-seq的count數(shù)據(jù)恋拷，基于負(fù)二項(xiàng)分布回歸

Combat

測試數(shù)據(jù)的設(shè)計(jì)矩陣為：

表達(dá)矩陣為：

# 提取batch的信息
batch = pheno$batch
# 建立對照組模型湘捎，這不也可以省略
modcombat = model.matrix(~1, data=pheno)

# 去除批次效應(yīng)
combat_edata1 = ComBat(dat=edata, batch=batch, mod=NULL, par.prior=TRUE,  prior.plots=FALSE)

這樣返回的就是去除批次效應(yīng)后的數(shù)據(jù)了

而Combat函數(shù)的核心代碼：

# 輸入數(shù)據(jù)
dat=edata
batch=batch
mod = NULL    ## 忽略設(shè)計(jì)矩陣
par.prior = TRUE
prior.plots = FALSE
mean.only = FALSE
ref.batch = NULL

# 定義出批次效應(yīng)的項(xiàng)
batch <- as.factor(batch)

# dat 為表達(dá)矩陣，s.data相當(dāng)于標(biāo)準(zhǔn)化后的表達(dá)矩陣
s.data <- (dat-stand.mean)/(sqrt(var.pooled) %*% t(rep(1,n.array)))

# 找到相同levels的sample舷胜，并歸類
n.batch <- nlevels(batch)
  batches <- list()
  for (i in 1:n.batch) {
    ## 把相同 levels 的 sample 歸為一類活翩，即 batches[[1]] 代表所有l(wèi)evels = 1 的sample
    batches[[i]] <- which(batch == levels(batch)[i])
  } # list of samples in each batch  

# 計(jì)算去除批次效應(yīng)的表達(dá)矩陣
bayesdata <- s.data
## batches 代表有批次效應(yīng)的sample
## 所有有批次效應(yīng)的sample需要扣除批次效應(yīng)的影響
for (i in batches){
    bayesdata[,i] <- (bayesdata[,i]-t(batch.design[i,]%*%gamma.star))/(sqrt(delta.star[j,])%*%t(rep(1,n.batches[j]))) 
    ### t(batch.design[i,]%*%gamma.star) 代表批次效應(yīng)所引起的基因表達(dá)量的差異材泄，這一部分是需要被扣除的
    j <- j+1
  }

# 最后乘上系數(shù)
 bayesdata <- (bayesdata*(sqrt(var.pooled)%*%t(rep(1,n.array))))+stand.mean # FIXME

而這段代碼的理論部分來自于文章《Adjusting batch effects in microarray expression data using empirical Bayes methods》

文章對去除批次效應(yīng)分為了三步：
step_1：Standardize the data

其中：

1. Yijg 代表 gene g for sample j from batch i
2. αg 代表 overall gene expression，即 gene g 在所有sample中平均表達(dá)量
3. X 代表 a design matrix for sample conditions峦树，代表實(shí)驗(yàn)設(shè)計(jì)所制定的矩陣
4. βg 代表對X矩陣的回歸系數(shù)，代表不同設(shè)計(jì)矩陣各個(gè)因素下的效應(yīng)值
5. X·βg 代表設(shè)計(jì)矩陣各個(gè)因素除平均表達(dá)量αg以外的額外的表達(dá)值急灭，對于sample j （因素 j）最后加上 αg 則代表sample j 中 gene g 的表達(dá)量
6. γig 代表 the additive and multiplicative batch effects of batch i for gene g
7. δig 代表 the additive and multiplicative batch effects of batch i for gene g

標(biāo)準(zhǔn)化

## 定義設(shè)計(jì)矩陣
## combine batch variable and covariates
design <- cbind(batchmod,mod)
## batchmod為與batch有關(guān)的設(shè)計(jì)矩陣埋凯，mod為實(shí)驗(yàn)設(shè)計(jì)矩陣
  
# 均值的計(jì)算
stand.mean <- stand.mean+t(tmp %*% B.hat)
# 方差的計(jì)算
var.pooled <- ((dat-t(design %*% B.hat))^2) %*% rep(1/n.array,n.array)
## B.hat代表β.hat
## tmp為與batch相關(guān)的設(shè)計(jì)矩陣，即與，如果不考慮實(shí)驗(yàn)設(shè)計(jì)因素换怖，那么tmp為0矩陣

Step_2: EB batch effect parameter estimates using parametric empirical priors
假設(shè)每一個(gè)元素所服從分布如下：

假設(shè)標(biāo)準(zhǔn)化后的表達(dá)矩陣?yán)锩娴脑刂?Zijg 服從正態(tài)分布彻坛，均值為 γig（γig為the additive and multiplicative batch effects of batch i for gene g），而這些參數(shù)估計(jì)需要用到貝葉斯后驗(yàn)來進(jìn)行估計(jì)

為批次效應(yīng)所引起的基因表達(dá)量的差異扛施，這一部分是需要被扣除的

Step_3: Adjust the data for batch effects

Combat_seq

這個(gè)函數(shù)主打RNA-seq的count數(shù)據(jù)昌阿，而它的下載方式需要用：

# 下載方式
devtools::install_github("zhangyuqing/sva-devel")

# 示例數(shù)據(jù)
count_matrix <- matrix(rnbinom(400, size=10, prob=0.1),nrow=50, ncol=8)
batch <- c(rep(1, 4), rep(2, 4))
adjusted <- ComBat_seq(count_matrix, batch=batch, group=NULL)

我們參考這篇文獻(xiàn)來分析它的原理：《ComBat-Seq: batch effect adjustment for RNA-Seq count data》
Combat_seq基于的原理是負(fù)二項(xiàng)分布回歸

這里的負(fù)二項(xiàng)分布指的是基因的表達(dá)服從負(fù)二項(xiàng)分布，如下圖

橫坐標(biāo)表示基因的表達(dá)量灶轰，縱坐標(biāo)表示的是基因表達(dá)量在某個(gè)區(qū)間范圍內(nèi)的頻率刷钢，那么整個(gè)count矩陣數(shù)據(jù)符合負(fù)二項(xiàng)分布，用具有這樣數(shù)據(jù)特點(diǎn)的矩陣做的回歸稱為負(fù)二項(xiàng)分布回歸

其中：

1. μijg 代表 gene g for sample j from batch i
2. αg 代表 overall gene expression，即 gene g 在所有sample中平均表達(dá)量
3. Xj 代表 a design matrix for sample conditions敞嗡，代表實(shí)驗(yàn)設(shè)計(jì)所制定的矩陣
4. βg 代表對X矩陣的回歸系數(shù)航背，代表不同設(shè)計(jì)矩陣各個(gè)因素下的效應(yīng)值
5. Xj·βg 代表設(shè)計(jì)矩陣各個(gè)因素除平均表達(dá)量αg以外的額外的表達(dá)值，對于sample j （因素 j）最后加上 αg 則代表sample j 中 gene g 的表達(dá)量
6. γgi 代表批次效應(yīng) i 影響 gene g 的表達(dá)值從而產(chǎn)生差異的部分（均值）
7. ** ?gi** 代表批次效應(yīng) i 影響 gene g 的表達(dá)值從而產(chǎn)生差異的部分（dispersion）
8. Nj 代表文庫大小
   

而關(guān)于Combat_seq的核心代碼：

# 讀入相關(guān)的內(nèi)容
## count_matrix為count矩陣
counts = count_matrix
## batch 為批次效應(yīng)向量
batch=batch
group=NULL
covar_mod=NULL
full_mod=TRUE
shrink=FALSE
shrink.disp=FALSE
gene.subset.n=NULL


# 廣義線性模型擬合數(shù)據(jù)箕肃，即建立真實(shí)的基因表達(dá)量
## 這里的design代表的是batch的設(shè)計(jì)矩陣
## dge_obj$counts為counts矩陣
## 參數(shù) phi_matrix 通過edgeR包中的函數(shù)estimateGLMCommonDisp進(jìn)行估計(jì)
glm_f2 <- glmFit.default(dge_obj$counts, design=design, dispersion=phi_matrix, 
                         offset=new_offset, prior.count=1e-4) 
  

# gamma_hat 為glm_f2 模型的系數(shù)今魔，代表不同batch對每個(gè)sample的影響
gamma_hat <- glm_f2$coefficients[, 1:n_batch]
# mu_hat為利用模型擬合的count值
mu_hat <- glm_f2$fitted.values
phi_hat <- do.call(cbind, genewise_disp_lst)

design矩陣

然而去除批次效應(yīng)所帶來影響需要將原始counts值，減去對應(yīng)batch所帶來的影響γgi：

mu_star <- matrix(NA, nrow=nrow(counts), ncol=ncol(counts))
## mu_hat 可以理解為count值
## vec2mat(gamma_star_mat[, jj], n_batches[jj])可以理解為batch的效應(yīng)值 γgi叫倍，只不過將它做成矩陣形式方便計(jì)算
for(jj in 1:n_batch){
    ## count值減去batch的效應(yīng)值 γgi，最終得到去除批次效應(yīng)的結(jié)果
    mu_star[, batches_ind[[jj]]] <- exp(log(mu_hat[, batches_ind[[jj]]])-vec2mat(gamma_star_mat[, jj], n_batches[jj]))
  }

最后對 mu_star 進(jìn)行整數(shù)化的矯正就可以得到最終結(jié)果了：