差异分析对比：DESeq2/edgeR/limma-voom¶

# === DESeq2 差异分析 ===
library(DESeq2)

# 从计数矩阵创建
dds <- DESeqDataSetFromMatrix(
  countData = count_matrix,               # 基因×样本的计数矩阵
  colData = sample_info,                  # 样本信息表（含分组）
  design = ~ condition                    # 实验设计公式
)

# 过滤低表达基因
keep <- rowSums(counts(dds) >= 10) >= 3   # 至少3个样本中计数≥10
dds <- dds[keep, ]

# 运行差异分析
dds <- DESeq(dds)                          # 一步完成：标准化+离散度+检验
res <- results(dds,
  contrast = c("condition", "Treatment", "Control"),  # 比较组
  alpha = 0.05                             # FDR阈值
)

# log2FC收缩（推荐，减少噪声）
res_shrink <- lfcShrink(dds,
  coef = "condition_Treatment_vs_Control", # 系数名
  type = "apeglm"                          # ★推荐的收缩方法
)

# 查看结果
summary(res)                               # 汇总统计
sig_genes <- subset(res, padj < 0.05 & abs(log2FoldChange) > 1)  # 筛选显著

# === edgeR 差异分析 ===
library(edgeR)

# 创建DGEList对象
dge <- DGEList(
  counts = count_matrix,                  # 计数矩阵
  group = sample_info$condition            # 分组信息
)

# 过滤低表达
keep <- filterByExpr(dge, group = dge$samples$group)  # 自动过滤
dge <- dge[keep, , keep.lib.sizes=FALSE]

# TMM标准化
dge <- calcNormFactors(dge, method="TMM")  # 计算标准化因子

# 设计矩阵
design <- model.matrix(~0 + condition, data=sample_info)  # 无截距设计
colnames(design) <- gsub("condition", "", colnames(design))

# 估计离散度
dge <- estimateDisp(dge, design)           # 估计离散度（common+trended+tagwise）

# 差异检验（准似然F检验，推荐）
fit <- glmQLFit(dge, design)               # 拟合模型
contrast <- makeContrasts(Treatment - Control, levels=design)  # 定义对比
qlf <- glmQLFTest(fit, contrast=contrast)  # F检验

# 查看结果
topTags(qlf, n=20)                          # 前20个差异基因
sig_edgeR <- topTags(qlf, n=Inf, p.value=0.05)$table  # 所有显著基因

# === limma-voom 差异分析 ===
library(limma)
library(edgeR)  # 需要edgeR做TMM标准化

# 创建DGEList并标准化
dge <- DGEList(counts = count_matrix)
dge <- calcNormFactors(dge, method="TMM")

# 过滤低表达
keep <- filterByExpr(dge, group = sample_info$condition)
dge <- dge[keep, , keep.lib.sizes=FALSE]

# 设计矩阵
design <- model.matrix(~0 + condition, data=sample_info)
colnames(design) <- gsub("condition", "", colnames(design))

# ★ voom转换（关键步骤）
v <- voom(dge, design, plot=TRUE)           # 估计均值-方差关系
# 或用voomWithQualityWeights（处理异质样本）
# v <- voomWithQualityWeights(dge, design, plot=TRUE)

# 拟合线性模型
fit <- lmFit(v, design)                    # 线性拟合
contrast <- makeContrasts(Treatment - Control, levels=design)
fit2 <- contrasts.fit(fit, contrast)        # 应用对比
fit2 <- eBayes(fit2)                        # 经验贝叶斯调节

# 查看结果
topTable(fit2, n=20)                        # 前20个
sig_limma <- topTable(fit2, n=Inf, p.value=0.05)  # 所有显著

# === 取三种方法的交集（更可靠） ===
library(VennDiagram)

# 提取各方法的显著基因（FDR<0.05, |logFC|>1）
deseq2_genes <- rownames(subset(res, padj<0.05 & abs(log2FoldChange)>1))
edger_genes <- rownames(subset(sig_edgeR, FDR<0.05 & abs(logFC)>1))
limma_genes <- rownames(subset(sig_limma, adj.P.Val<0.05 & abs(logFC)>1))

# 画韦恩图
venn.diagram(
  list(DESeq2=deseq2_genes, edgeR=edger_genes, limma=limma_genes),
  filename="DEG_venn.png",
  fill=c("red", "blue", "green"),
  alpha=0.5
)

# 取交集（最可靠的差异基因）
robust_DEGs <- Reduce(intersect, list(deseq2_genes, edger_genes, limma_genes))
cat("三方法交集基因数:", length(robust_DEGs), "\n")

# === 通用可视化 ===
library(EnhancedVolcano)

# 火山图（DESeq2结果）
EnhancedVolcano(res,
  lab = rownames(res),                     # 基因名
  x = 'log2FoldChange',                   # x轴
  y = 'padj',                             # y轴
  pCutoff = 0.05,                          # p值阈值
  FCcutoff = 1,                            # logFC阈值
  title = 'Treatment vs Control'
)

# MA图
plotMA(res, ylim=c(-5,5), main="MA Plot")

# 热图（显著基因）
library(pheatmap)
sig_matrix <- assay(vst(dds))[rownames(sig_genes), ]  # VST标准化
pheatmap(sig_matrix, scale="row",                       # 行标准化
         show_rownames=FALSE,
         annotation_col=sample_info["condition"])

# === 差异分析速查 ===
# DESeq2
dds <- DESeqDataSetFromMatrix(counts, colData, ~condition)
dds <- DESeq(dds); res <- results(dds, contrast=c("condition","T","C"))

# edgeR
dge <- DGEList(counts); dge <- calcNormFactors(dge)
dge <- estimateDisp(dge, design); fit <- glmQLFit(dge, design)
qlf <- glmQLFTest(fit, contrast=con)

# limma-voom
v <- voom(dge, design); fit <- lmFit(v, design)
fit2 <- eBayes(contrasts.fit(fit, con)); topTable(fit2)

# 选择: 小样本→DESeq2 | 大样本→limma | 低丰度→edgeR | 发表→三者交集