RNA-Seq 差异表达分析通常需要三类输入：基因表达计数矩阵、样本分组信息和基因注释。DESeq2 与 edgeR 都以原始整数 counts 为输入，通过负二项分布建模分析组间差异，但在标准化、离散度估计和统计检验上采用了不同实现。
一、准备数据
library(dplyr)
library(tibble)
library(DESeq2)
library(edgeR)
counts <- read.csv("airway_counts.csv", check.names = FALSE)
metadata <- read.csv("airway_metadata.csv", row.names = 1)
anno <- read.csv("annotables_grch38.csv")表达矩阵应当每行对应一个基因、每列对应一个样本，metadata 每行对应一个样本。正式分析前首先确认二者顺序一致：
stopifnot(all(colnames(counts)[-1] == rownames(metadata)))
stopifnot(all(counts[-1] == round(as.matrix(counts[-1]))))DESeq2 和 edgeR 需要原始整数计数，不能直接使用 TPM、FPKM、CPM 或经过缩放的数据。原始代码中的文件名是 airway_scaledcounts.csv，因此必须先确认文件内容；如果确实是 scaled counts，应回到比对或定量结果重新取得 raw counts。
基因注释可以通过基因 ID 合并，但不建议简单使用 distinct(symbol) 丢弃重复 symbol，因为多个 Ensembl ID 可能映射到同一个 symbol。更稳妥的方式是保留稳定基因 ID 作为行名，将 symbol 作为结果注释；如果业务确实需要按 symbol 合并，应明确采用求和或其他规则。
二、使用 DESeq2
count_matrix <- counts |>
column_to_rownames("gene_id") |>
as.matrix()
metadata$dex <- relevel(factor(metadata$dex), ref = "control")
dds <- DESeqDataSetFromMatrix(
countData = count_matrix,
colData = metadata,
design = ~ dex
)
dds <- DESeq(dds)
resultsNames(dds)
res <- results(dds, contrast = c("dex", "treated", "control"))
res <- as.data.frame(res)
res$gene_id <- rownames(res)DESeq() 已经依次完成 size factor、离散度和模型估计，一般不需要在它之后再次调用 estimateSizeFactors()。如果需要查看标准化计数，可以使用：
normalized_counts <- counts(dds, normalized = TRUE)三组实验仍然使用同一设计公式，只需要重新构建包含正确分组的 dds，再分别提取对比结果：
metadata$dex <- factor(metadata$dex, levels = c("a", "b", "c"))
dds3 <- DESeqDataSetFromMatrix(
countData = count_matrix,
colData = metadata,
design = ~ dex
)
dds3 <- DESeq(dds3)
res_b_vs_a <- results(dds3, contrast = c("dex", "b", "a"))
res_c_vs_a <- results(dds3, contrast = c("dex", "c", "a"))只修改外部 metadata$dex 不会自动改变已经创建的 dds，因此分组改变后应重新创建 DESeqDataSet，或者显式修改 colData(dds) 后重新拟合。
三、使用 edgeR
group <- factor(metadata$dex, levels = c("control", "treated"))
dge <- DGEList(counts = count_matrix, group = group)
dge <- calcNormFactors(dge)
design <- model.matrix(~ 0 + group)
colnames(design) <- levels(group)
dge <- estimateDisp(dge, design, robust = TRUE)
plotMDS(dge)
plotBCV(dge)
fit <- glmQLFit(dge, design, robust = TRUE)
contrast <- makeContrasts(treated - control, levels = design)
test <- glmQLFTest(fit, contrast = contrast)
deg <- topTags(test, n = Inf)$table
write.csv(deg, "edger-treated-vs-control.csv")plotMDS() 用于观察样本间整体距离，plotBCV() 用于查看生物学变异。makeContrasts(treated - control) 明确表示 treated 相对 control 的变化方向；反过来写会使 logFC 符号相反。
四、结果解释与检查
DESeq2 常用 log2FoldChange、pvalue 和 padj，edgeR 常用 logFC、PValue 和 FDR。两种工具的显著基因不一定完全相同，比较结果时要保持输入基因集合、样本分组和 contrast 方向一致。
分析完成后至少检查：
- counts 是否为非负整数。
- counts 列名是否与 metadata 行名完全对应。
- 参考组和对比方向是否符合问题定义。
- 是否存在极低表达基因，是否需要预过滤。
- 多重检验校正后是否仍然显著。
- 样本在 MDS/PCA 中是否按实验组分离，是否存在离群点。
原始代码展示了两套工具的完整轮廓，但数据对象在修改分组后没有重新创建、重复 symbol 被直接删除、contrast 方向和拼写也存在问题。将数据对应关系固定下来，再分别运行两套模型，结果才具有可解释性。
