DESeq2的shrink(以apeglm为例)

DESeq2对logFC会有一个shrink的过程，我们以apeglm为例：
step 1：创建deseq对象

# 创建deseq对象
library(DESeq2)

set.seed(1)
dds <- makeExampleDESeqDataSet(n=500,betaSD=1)
dds <- DESeq(dds,betaPrior = F)
res <- results(dds)

step 2：apeglm分析
接下来对lfcshrink的源码进行分析

dds=dds
res=res
type="apeglm"
lfcThreshold=0
svalue=FALSE
returnList=FALSE
format=c("DataFrame","GRanges","GRangesList")
saveCols=NULL
apeAdapt=TRUE
apeMethod="nbinomCR"
parallel=FALSE
quiet=FALSE
coefNum = 2

# 获得利用标准化因子进行标准化后的矩阵   
Y <- counts(dds)
# 定义设计矩阵
design <- model.matrix(design(dds), data=colData(dds))

# 计算dds中基因的离散 α
disps <- dispersions(dds)

# 计算dds中每个sample的标准化因子
offset <- matrix(log(sizeFactors(dds)),
                   nrow=nrow(dds), ncol=ncol(dds), byrow=TRUE)

# 设置权重矩阵
weights <- matrix(1, nrow=nrow(dds), ncol=ncol(dds))
mle <- log(2) * cbind(res$log2FoldChange, res$lfcSE)
log.lik <- NULL
apeT <- NULL

# 进行apeglm回归
fit <- apeglm::apeglm(Y=Y, x=design, log.lik=log.lik, param=disps,
                        coef=coefNum, mle=mle, threshold=apeT,
                        weights=weights, offset=offset,
                        method=apeMethod)
conv <- fit$diag[,"conv"]
# 定义 logFC
res$log2FoldChange <- log2(exp(1)) * fit$map[,coefNum]
res$lfcSE <- log2(exp(1)) * fit$sd[,coefNum]

其中，条件比较为：

设计矩阵design为：

将condition A作为截距项（基准线），如果对表达量取对数，那么线性模型做变形可以理解为：
对于某基因：

其中：

1. E_{condition B} 代表条件B下某基因的表达量；E_{condition A} 代表条件A下某基因的表达量
2. β代表效应值
3. X代表设计矩阵（design）

由上图，由于存在条件B在设计矩阵中表示为 1，所以condition B vs condition A 的logFC值为 β

因此只需要估计出β值即可得到对数变化值，而apeglm采用的是：

贝叶斯方式估计β值（MAP矩阵，maximum a posteriori），MAP矩阵的值即为β值，因此shrink后的logFC为：

res$log2FoldChange <- log2(exp(1)) * fit$map[,coefNum]

fit$map[,coefNum]：