如何将edgeR生成的TMM标准化RNA-seq文件转换为差异表达基因列表

前提说明 #

你已获得edgeR生成的TMM标准化表达矩阵，行对应基因、列对应样本，首列为基因名，包含27个A组样本、18个B组样本，需计算B组相对A组的差异表达基因，采用BH法多重检验校正，筛选阈值FDR<0.05。

完整代码 #

# 1. 加载所需R包，如未安装先执行安装命令
# install.packages("BiocManager")
# BiocManager::install(c("edgeR", "limma"))
library(limma)
library(edgeR)

# 2. 读入TMM标准化表达文件，注意修改文件路径为你的本地路径
expr_df <- read.csv("你的TMM标准化文件路径.csv", row.names = 1, header = TRUE, check.names = FALSE)
# 若为txt格式用：expr_df <- read.table("你的TMM标准化文件路径.txt", row.names = 1, header = TRUE, sep = "\t", check.names = FALSE)

# 3. 构建分组信息，**必须确保样本顺序和表达矩阵的列顺序完全对应**
# 如下假设表达矩阵前27列为A组样本，后18列为B组样本，可根据实际顺序调整
group <- factor(c(rep("A", 27), rep("B", 18)), levels = c("A", "B"))

# 4. 构建设计矩阵与对比矩阵，指定比较组为B vs A
design <- model.matrix(~0 + group)
colnames(design) <- levels(group)
contrast_mat <- makeContrasts(B_vs_A = B - A, levels = design)

# 5. 差异表达分析
fit <- lmFit(expr_df, design)
fit_contrast <- contrasts.fit(fit, contrast_mat)
fit_ebayes <- eBayes(fit_contrast, trend = TRUE)

# 6. 提取全基因差异分析结果，校正方法指定为Benjamini-Hochberg
all_deg_res <- topTable(fit_ebayes, coef = "B_vs_A", adjust = "BH", number = Inf, sort.by = "none")
# 给结果添加基因名列
all_deg_res$gene <- rownames(all_deg_res)
# 调整列顺序为：基因名、logFC、p值、校正后p值（FDR）
all_deg_res <- all_deg_res[, c("gene", "logFC", "P.Value", "adj.P.Val")]

# 7. 筛选FDR<0.05的显著差异表达基因
sig_deg_res <- all_deg_res[all_deg_res$adj.P.Val < 0.05, ]

# 8. 输出结果文件
write.table(all_deg_res, "全基因差异分析结果.txt", sep = "\t", row.names = FALSE, quote = FALSE)
write.table(sig_deg_res, "显著差异表达基因结果_FDR小于0.05.txt", sep = "\t", row.names = FALSE, quote = FALSE)

注意事项 #

• 请务必核对分组的样本顺序和表达矩阵的列顺序一致，否则会得到完全错误的分析结果
• 若你的输入文件为原始计数矩阵而非经过TMM标准化的log转换表达值，需要使用edgeR标准差异分析流程，上述代码不适用
• 可根据研究需求自行调整FDR筛选阈值，也可额外添加log2FC的阈值（比如|logFC|>1）进一步缩小差异基因范围

内容的提问来源于stack exchange，提问作者hamid