有一个朋友问了我一个问题，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/
结论：
- 你以为的未必真的是你以为
- 相互交流可以提高水平
- 一定要学会搜索
