DESeq2的baseMean和log2FoldChange是如何得到的？

有一个朋友问了我一个问题，DESeq2的baseMean是如何计算？我最初都是认为baseMean计算的是对照组的样本标准化counts的均值。由于我在分析结果里还会提供所有样本的标准化的counts，所以这个baseMean我也没有太过在意。在他提出问题后，我做了如下的探索，以R包airway为例

首先提取出差异表达的baseMean

library(airway)
data("airway")
se <- airway

library("DESeq2")
ddsSE <- DESeqDataSet(se, design = ~ cell + dex)
ddsSE

dds <- ddsSE

keep <- rowSums(counts(dds)) >= 10
dds <- dds[keep,]

dds$condition <- colData(dds)$dex
dds$condition <- relevel(dds$condition, ref = "untrt")
dds <- DESeq(dds)
res <- results(dds)
res

部分结果

以ENSG00000000003为例，提取该基因组在所有样本的表达量

normlized_counts <- counts(dds, normalized=TRUE)

untrt <- normlized_counts['ENSG00000000003', c(1, 3, 5, 7)]
trt <-  normlized_counts['ENSG00000000003', c(2, 4, 6, 8)]

注: 这里的对照组和实验组是间隔存放。

两者的均值如下，

> mean(untrt)
[1] 801.4966
> mean(trt)
[1] 615.6991

不难发现结果都和原来的708有差异。但是好像将两者进行平均后，结果就是708。

> mean(c(untrt, trt))
[1] 708.5979

结果验证了我的猜想。当然我还用关键字 "deseq2 basemean calculation"进行了搜索，链接如下

https://support.bioconductor.org/p/75244/

而log2FoldChange计算也不是简单的用标准化的counts进行计算，因为计算的时候需要考虑零值以及其他效应，参考链接为https://support.bioconductor.org/p/77021/

结论：

1. 你以为的未必真的是你以为
2. 相互交流可以提高水平
3. 一定要学会搜索