关于sva包

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

其中关于去除批次效应的函数有Combat和Combat_seq两个，而Combat主要面对的是带有小数的数据（比方说芯片数据），基于的是贝叶斯原理；而Combat_seq主打RNA-seq的count数据，基于负二项分布回归

Combat

测试数据的设计矩阵为：

表达矩阵为：

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

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

这样返回的就是去除批次效应后的数据了

而Combat函数的核心代码：

# 输入数据
dat=edata
batch=batch
mod = NULL    ## 忽略设计矩阵
par.prior = TRUE
prior.plots = FALSE
mean.only = FALSE
ref.batch = NULL

# 定义出批次效应的项
batch <- as.factor(batch)

# dat 为表达矩阵，s.data相当于标准化后的表达矩阵
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]] 代表所有levels = 1 的sample
    batches[[i]] <- which(batch == levels(batch)[i])
  } # list of samples in each batch  

# 计算去除批次效应的表达矩阵
bayesdata <- s.data
## batches 代表有批次效应的sample
## 所有有批次效应的sample需要扣除批次效应的影响
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) 代表批次效应所引起的基因表达量的差异，这一部分是需要被扣除的
    j <- j+1
  }

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

而这段代码的理论部分来自于文章《Adjusting batch effects in microarray expression data using empirical Bayes methods》

文章对去除批次效应分为了三步：
step_1：Standardize the data

其中：

1. Yijg 代表 gene g for sample j from batch i
2. αg 代表 overall gene expression，即 gene g 在所有sample中平均表达量
3. X 代表 a design matrix for sample conditions，代表实验设计所制定的矩阵
4. βg 代表对X矩阵的回归系数，代表不同设计矩阵各个因素下的效应值
5. X·βg 代表设计矩阵各个因素除平均表达量αg以外的额外的表达值，对于sample j （因素 j）最后加上 αg 则代表sample j 中 gene g 的表达量
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

标准化

## 定义设计矩阵
## combine batch variable and covariates
design <- cbind(batchmod,mod)
## batchmod为与batch有关的设计矩阵，mod为实验设计矩阵
  
# 均值的计算
stand.mean <- stand.mean+t(tmp %*% B.hat)
# 方差的计算
var.pooled <- ((dat-t(design %*% B.hat))^2) %*% rep(1/n.array,n.array)
## B.hat代表β.hat
## tmp为与batch相关的设计矩阵，即与，如果不考虑实验设计因素，那么tmp为0矩阵

Step_2: EB batch effect parameter estimates using parametric empirical priors
假设每一个元素所服从分布如下：

假设标准化后的表达矩阵里面的元素值 Zijg 服从正态分布，均值为 γig（γig为the additive and multiplicative batch effects of batch i for gene g），而这些参数估计需要用到贝叶斯后验来进行估计

为批次效应所引起的基因表达量的差异，这一部分是需要被扣除的

Step_3: Adjust the data for batch effects

Combat_seq

这个函数主打RNA-seq的count数据，而它的下载方式需要用：

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

# 示例数据
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)

我们参考这篇文献来分析它的原理：《ComBat-Seq: batch effect adjustment for RNA-Seq count data》
Combat_seq基于的原理是负二项分布回归

这里的负二项分布指的是基因的表达服从负二项分布，如下图

横坐标表示基因的表达量，纵坐标表示的是基因表达量在某个区间范围内的频率，那么整个count矩阵数据符合负二项分布，用具有这样数据特点的矩阵做的回归称为负二项分布回归

其中：

1. μijg 代表 gene g for sample j from batch i
2. αg 代表 overall gene expression，即 gene g 在所有sample中平均表达量
3. Xj 代表 a design matrix for sample conditions，代表实验设计所制定的矩阵
4. βg 代表对X矩阵的回归系数，代表不同设计矩阵各个因素下的效应值
5. Xj·βg 代表设计矩阵各个因素除平均表达量αg以外的额外的表达值，对于sample j （因素 j）最后加上 αg 则代表sample j 中 gene g 的表达量
6. γgi 代表批次效应 i 影响 gene g 的表达值从而产生差异的部分（均值）
7. ** ϕgi** 代表批次效应 i 影响 gene g 的表达值从而产生差异的部分（dispersion）
8. Nj 代表文库大小
   

而关于Combat_seq的核心代码：

# 读入相关的内容
## count_matrix为count矩阵
counts = count_matrix
## batch 为批次效应向量
batch=batch
group=NULL
covar_mod=NULL
full_mod=TRUE
shrink=FALSE
shrink.disp=FALSE
gene.subset.n=NULL


# 广义线性模型拟合数据，即建立真实的基因表达量
## 这里的design代表的是batch的设计矩阵
## dge_obj$counts为counts矩阵
## 参数 phi_matrix 通过edgeR包中的函数estimateGLMCommonDisp进行估计
glm_f2 <- glmFit.default(dge_obj$counts, design=design, dispersion=phi_matrix, 
                         offset=new_offset, prior.count=1e-4) 
  

# gamma_hat 为glm_f2 模型的系数，代表不同batch对每个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矩阵

然而去除批次效应所带来影响需要将原始counts值，减去对应batch所带来的影响γgi：

mu_star <- matrix(NA, nrow=nrow(counts), ncol=ncol(counts))
## mu_hat 可以理解为count值
## vec2mat(gamma_star_mat[, jj], n_batches[jj])可以理解为batch的效应值 γgi，只不过将它做成矩阵形式方便计算
for(jj in 1:n_batch){
    ## count值减去batch的效应值 γgi，最终得到去除批次效应的结果
    mu_star[, batches_ind[[jj]]] <- exp(log(mu_hat[, batches_ind[[jj]]])-vec2mat(gamma_star_mat[, jj], n_batches[jj]))
  }

最后对 mu_star 进行整数化的矫正就可以得到最终结果了：