如何将EdgeR的TMM标准化计数作为输入用于DESeq2分析？
我有多个不同实验条件下的RNA-seq样本，完成测序并比对到参考基因组后，合并原始计数得到如下数据框：
> df_merge T0 DJ21 DJ24 DJ29 DJ32 Rec2 Rec6 Rec9 G10 421 200 350 288 284 198 314 165 G1000 17208 10608 11720 11421 10142 10768 10331 6121 G10000 37 16 19 21 28 12 9 4 G10002 45 13 44 27 12 35 74 14 G10003 136 79 162 429 184 112 192 162 G10004 54 162 73 169 102 300 429 180 G10006 1 0 1 0 0 0 0 0 G10007 3 4 7 2 1 1 1 0 G1001 9030 8366 10608 13604 9808 10654 11663 7985 ... ... ... ... ... ... ... ... ...
我选择用edgeR执行TMM标准化（DESeq2本身不内置支持该方法），使用如下脚本完成计算：
## Normalisation by the TMM method (Trimmed Mean of M-value) dge <- DGEList(df_merge) # 从计数数据创建DGEList对象 dge2 <- calcNormFactors(dge, method = "TMM") # 计算TMM标准化因子
计算后得到的标准化因子如下：
> dge2$samples group lib.size norm.factors T0 1 129884277 1.1108130 DJ21 1 110429304 0.9453988 DJ24 1 126410256 1.0297216 DJ29 1 123008035 1.0553169 DJ32 1 118968544 0.9927826 Rec2 1 119000510 0.9465131 Rec6 1 114775318 1.0053686 Rec9 1 90693946 0.9275454
我最初尝试用上述标准化因子处理原始计数，得到标准化后的伪计数：
# 用cpm函数计算标准化伪计数并存入数据框 pseudo_TMM <- log2(cpm(dge2) + 1) df_TMM <- melt(pseudo_TMM, id = rownames(raw_counts_wn)) names(df_TMM)[1:2] <- c ("id", "sample") df_TMM$method <- rep("TMM", nrow(df_TMM))
最终得到的TMM标准化后计数如下：
> pseudo_TMM T0 DJ21 DJ24 DJ29 DJ32 Rec2 Rec6 Rec9 G10 1.970115581 1.54384913 1.88316953 1.68642670 1.76745996 1.46356074 1.89575666 1.56628879 G1000 6.910138402 6.68101996 6.50839579 6.47542172 6.44077248 6.59395683 6.50032388 6.20481983 G10000 0.329354263 0.20571418 0.19656414 0.21632677 0.30692404 0.14605339 0.10835095 0.06701850 G10002 0.391657436 0.16931112 0.42010652 0.27261134 0.13960084 0.39037793 0.71483462 0.22209164 G10003 0.958011321 0.81287356 1.16642722 2.10593537 1.35494357 0.99592405 1.41354030 1.54881003 G10004 0.458675608 1.35147467 0.64230087 1.20281148 0.89809414 1.87320592 2.23810756 1.65064058 G10006 0.009964976 0.00000000 0.01104103 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 G10007 0.029690785 0.05424318 0.07556948 0.02205789 0.01216343 0.01275200 0.01244875 0.00000000 G1001 5.990679797 6.34224022 6.36623615 6.72515956 6.39302663 6.57876150 6.67346174 6.58377191 ... ... ... ... ... ... ... ... ...
我平时做差异基因表达分析习惯用DESeq2的DESeqDataSetFromHTSeqCount()和DESeq()函数，该流程默认执行RLE标准化。现在我希望直接使用已经完成TMM标准化的数据开展DESeq2差异分析，已知可以通过DESeqDataSetFromMatrix()函数从矩阵构建DESeqDataSet对象，求可行的操作方案。
首先明确一个核心避坑点：绝对不能把你当前计算得到的log2转换后TMM伪计数（pseudo_TMM）直接输入DESeq2。DESeq2的差异分析模型基于负二项分布构建，要求输入必须是未经转换的原始整数计数，输入log转换后的连续型数值会直接违背模型的基本假设，得到的差异结果完全不可靠。
你不需要把标准化后的计数喂给DESeq2，只需要把edgeR计算得到的TMM标准化因子手动传入DESeq2，让软件跳过自带的RLE标准化因子估计步骤即可，后续的离散度估计、模型拟合、差异检验流程和你平时的操作完全一致，具体步骤如下：
- 第一步：从edgeR的计算结果中转换得到DESeq2可用的标准化因子
edgeR算出来的norm.factors需要和原始文库大小相乘得到有效库大小，再做中心化处理（让所有因子的几何均值为1，和DESeq2内部sizeFactor的计算逻辑对齐），代码如下：# 计算TMM对应的有效文库大小 tmm_eff_lib <- dge2$samples$lib.size * dge2$samples$norm.factors # 中心化得到符合DESeq2要求的sizeFactors tmm_size_factors <- tmm_eff_lib / exp(mean(log(tmm_eff_lib))) - 第二步：用原始整数计数矩阵构建DESeqDataSet对象
这里必须用你最开始得到的未做任何转换的原始计数表df_merge，不能用cpm或者log转换后的表，同时按你的实验设计准备样本分组信息：# 请将下面的condition替换为你实际的实验分组 coldata <- data.frame( row.names = colnames(df_merge), condition = factor(c("T0", "DJ", "DJ", "DJ", "DJ", "Rec", "Rec", "Rec")) ) dds <- DESeqDataSetFromMatrix( countData = df_merge, colData = coldata, design = ~ condition ) - 第三步：手动替换DESeq2的标准化因子为TMM计算结果
这一步是核心，赋值后DESeq2就不会再运行默认的RLE标准化计算，直接使用你传入的TMM因子做文库校正：sizeFactors(dds) <- tmm_size_factors - 第四步：按常规流程运行差异分析即可
dds <- DESeq(dds) # 后续提取差异结果、做下游分析的操作和你之前用的标准DESeq2流程完全一致 # 示例：提取DJ组相对T0组的差异结果 res_dj_vs_t0 <- results(dds, contrast = c("condition", "DJ", "T0"))
补充说明：如果你需要得到TMM标准化后的计数用于可视化（比如PCA、热图），可以在给
dds赋值完sizeFactors之后，用counts(dds, normalized=TRUE)提取线性尺度的标准化计数，不要直接用之前算的log2(cpm+1)做差异分析输入，这类转换后的值仅适合用于探索性可视化。
内容的提问来源于stack exchange，提问作者Paul Sourbé
