代码

我们先上代码：

library(zebrafishRNASeq)
library(sva)
library(RUVSeq)

data(zfGenes)
filter = apply(zfGenes, 1, function(x) length(x[x>5])>=2)
filtered = zfGenes[filter,]
genes = rownames(filtered)[grep("^ENS", rownames(filtered))]
controls = grepl("^ERCC", rownames(filtered))
spikes =  rownames(filtered)[grep("^ERCC", rownames(filtered))]
group = as.factor(rep(c("Ctl", "Trt"), each=3))
set = newSeqExpressionSet(as.matrix(filtered),
                          phenoData = data.frame(group, row.names=colnames(filtered)))
dat0 = counts(set)

# 设计矩阵
## 生物学差异的设计矩阵
mod1 = model.matrix(~group)
## 无差异的设计矩阵
mod0 = cbind(mod1[,1])

# DESeq2
condition <- factor(c(rep("control",3),rep("treat",3)), levels = c("control","treat"))

colData <- data.frame(row.names=colnames(dat0), condition)
dds <- DESeqDataSetFromMatrix(dat0, colData, design= ~ condition)
dds <- DESeq(dds)


dat = counts(dds, normalized=TRUE)
dds_result = data.frame(results(dds, contrast=c("condition","treat","control")))


## Estimate batch with svaseq (unsupervised)
svseq = svaseq(dat0,mod1,mod0,n.sv=2)

## 然后接上DESeq的标准流程
ddssva <- dds
ddssva$SV1 <- svseq$sv[,1]
ddssva$SV2 <- svseq$sv[,2]
design(ddssva) <- ~ SV1 + SV2 + condition

library(magrittr)
ddssva %<>% DESeq

dat_sva = counts(ddssva, normalized=TRUE)
dds_result_sva = data.frame(results(ddssva, contrast=c("condition","treat","control")))

我们注意到，用DESeq矫正批次效应中有一步，利用svaseq()来构造替代变量完成批次效应的矫正

而对于矫正前的log2(FC)

对于矫正后的log2(FC)

最终的效果是log2(FC)有一定的改变

svaseq

我们打开svaseq函数的源代码，目标集中在irwsva.build函数

n.sv = 2
dat = dat0
mod = mod1
mod0 = mod0
  
n <- ncol(dat)
m <- nrow(dat)
  
Id <- diag(n)
# resid为扣除生物学设计所带来的差异的表达矩阵
resid <- dat %*% (Id - mod %*% solve(t(mod) %*% mod) %*% t(mod))  
# 对 resid 进行SVD分解
uu <- eigen(t(resid)%*%resid)
ndf <- n - dim(mod)[2]
  
pprob <- rep(1,m)
one <- rep(1,n)
Id <- diag(n)
df1 <- dim(mod)[2] + n.sv
df0 <- dim(mod0)[2]  + n.sv

## 以下的循环代表多次筛选取均值
B = 5
for(i in 1:B){
    # 计算生物学因素对基因表达差异的影响
    mod.b <- cbind(mod,uu$vectors[,1:n.sv])
    mod0.b <- cbind(mod0,uu$vectors[,1:n.sv])
    ptmp <- f.pvalue(dat,mod.b,mod0.b)
    pprob.b <- (1-edge.lfdr(ptmp))
    
    # 计算unmodel factors对基因表达差异的影响
    mod.gam <- cbind(mod0,uu$vectors[,1:n.sv])
    mod0.gam <- cbind(mod0)
    ptmp <- f.pvalue(dat,mod.gam,mod0.gam)
    pprob.gam <- (1-edge.lfdr(ptmp))
    ## pprob.gam代表unmodel factors对基因表达差异的影响
    ## pprob.b代表生物学因素对基因表达差异的影响
    ## 1-pprob.b代表unmodel factors对基因表达差异的影响
    ## pprob为基因数目维的向量
    pprob <- pprob.gam*(1-pprob.b)
    dats <- dat*pprob
    # 数据中心化
    dats <- dats - rowMeans(dats)
    uu <- eigen(t(dats)%*%dats)
}
  
sv = svd(dats)$v[,1:n.sv, drop=FALSE]

其中：

pprob.gam代表unmodel factors对基因表达差异的影响

pprob.b代表生物学因素对基因表达差异的影响

1-pprob.b代表unmodel factors对基因表达差异的影响

值得注意的是：

resid <- dat %*% (Id - mod %*% solve(t(mod) %*% mod) %*% t(mod))

resid为扣除生物学设计所带来的差异的表达矩阵，然后对该矩阵SVD分解：

uu <- eigen(t(resid)%*%resid)

那么uu特征向量为

             sv1          sv2         sv3        sv4       sv5        sv6
Ctl1   0.29441072  0.402479876  0.50127306  0.4083188 0.5472996 -0.1838383
Ctl3  -0.22582801  0.012722701  0.09208784 -0.7791189 0.5472996 -0.1838383
Ctl5  -0.06858271 -0.415202578 -0.59336090  0.3708001 0.5472996 -0.1838383
Trt9  -0.50202545  0.573280615 -0.26598909  0.1234352 0.1838383  0.5472996
Trt11  0.74069432  0.007025239 -0.24253355 -0.2432421 0.1838383  0.5472996
Trt13 -0.23866886 -0.580305854  0.50852263  0.1198069 0.1838383  0.5472996

每一列代表unmodel factors，里面的元素代表不同的unmodel factors对每一个sample的影响
关于SVD分解的生物学意义可以参考：SVD分解在生物学中的意义
只不过这里对sample的影响变为了unmodel factors

然后for循环里面相当于对特征向量 uu 多次操作取平均值

pprob <- pprob.gam*(1-pprob.b)
dats <- dat*pprob

上面语句相当于得到了仅有 unmodel factors 影响下的仅有表达矩阵 dats，因此再次SVD分解，最终取得影响 sample 的特征向量（取前两个sv），而这两个sv就可以看作为 unmodel factors 对sample的影响，把这个加入 DESeq的设计矩阵中去除即可

关于f.pvalue和edge.lfdr函数

其中，函数f.pvalue：

dat = dat0

f.pvalue <- function(dat,mod,mod0){
  n <- dim(dat)[2]
  m <- dim(dat)[1]
  df1 <- dim(mod)[2]
  df0 <- dim(mod0)[2]
  p <- rep(0,m)
  Id <- diag(n)
  
  # resid为扣除因素带来的差异的表达矩阵
  resid <- dat %*% (Id - mod %*% solve(t(mod) %*% mod) %*% t(mod))
  # 计算残差矩阵每个元素平方和,作为F分布的参数
  rss1 <- rowSums(resid*resid)
  
  # resid0为未扣除因素的表达矩阵
  resid0 <- dat %*% (Id - mod0 %*% solve(t(mod0) %*% mod0) %*% t(mod0))
  # 计算残差矩阵每个元素平方和,作为F分布的参数
  rss0 <- rowSums(resid0*resid0)
 
  #构造F分布统计量,一个组别是(rss0 - rss1)/(df1-df0),另一个组别是(rss1/(n-df1)
  ## (rss0 - rss1)为参数的分布代表潜在因素所带来的差异的表达值分布
  ### 潜在因素理解为 mod(resid) 设计矩阵以外的因素；而因素理解为mod(resid) 设计矩阵设定的因素
  fstats <- ((rss0 - rss1)/(df1-df0))/(rss1/(n-df1))
  ## 由表达式来看,如果分子分母比值越接近 1,代表因素对该基因的影响越小,反之影响越大
  # 计算F分布p值
  p <-  1-pf(fstats,df1=(df1-df0),df2=(n-df1))
  ## 如果p越大,代表因素对该基因的影响越大
  ## 如果p越小,代表因素对该基因的影响越小
}

其中：

(rss0 - rss1)为参数的分布代表潜在因素所带来的差异的表达值分布

潜在因素理解为 mod(resid) 设计矩阵以外的因素；而因素理解为mod(resid) 设计矩阵设定的因素

对于 ((rss0 - rss1)/(df1-df0))/(rss1/(n-df1)) ，代表因素对该基因的影响越小,反之影响越大

对于函数edge.lfdr，它的主要功能是：

adj=1.5
eps=10^-8
lambda=0.8
  
p = ptmp
## 如果p越大,代表因素对该基因的影响越大
## 如果p越小,代表因素对该基因的影响越小
pi0 <- mean(p >= lambda)/(1 - lambda)
pi0 <- min(pi0, 1)
  
n <- length(p)
p <- pmax(p, eps)
p <- pmin(p, 1-eps)
# 将因素对该基因的影响权重累积成个正态分布
# 这个分布的均值是0
x <- qnorm(p)
myd <- density(x, adjust=adj)
mys <- smooth.spline(x=myd$x, y=myd$y)
# 标准化
y <- predict(mys, x)$y
lfdr <- pi0*dnorm(x)/y
  
# 做筛选
lfdr[lfdr > 1] <- 1
lfdr <- lfdr[order(p)]
lfdr <- mono(lfdr)
lfdr <- lfdr[rank(p)]