R语言实战：用edgeR和limma绘制专业级RNA-seq差异表达MA图
在生物信息学分析中，RNA-seq数据的可视化是理解差异表达基因分布的关键环节。MA图作为一种经典的可视化工具，能够直观展示基因表达变化与平均表达水平的关系，帮助研究者快速识别显著差异表达的基因。本文将深入探讨如何利用R语言中的edgeR和limma包，从原始计数数据出发，一步步完成差异分析并生成发表级的MA图。
1. 准备工作与环境配置
在开始绘制MA图之前，需要确保已经完成了RNA-seq数据的预处理和差异表达分析。假设我们已经获得了原始的基因计数矩阵，并进行了基本的质量控制和标准化处理。
首先，我们需要安装并加载必要的R包：
# 安装Bioconductor管理器（如未安装） if (!require("BiocManager", quietly = TRUE)) install.packages("BiocManager") # 安装edgeR和limma包 if (!require("edgeR", quietly = TRUE)) BiocManager::install("edgeR") if (!require("limma", quietly = TRUE)) BiocManager::install("limma") # 加载必要的包 library(edgeR) library(limma) library(ggplot2) # 用于后续图形美化对于RNA-seq数据，我们通常使用edgeR或DESeq2进行差异表达分析。本文将以edgeR为例，展示完整的分析流程。假设我们有一个包含两组样本（例如对照组和实验组）的基因计数矩阵，数据存储在一个名为"count_data"的数据框中，行名为基因ID，列名为样本名，同时有一个分组向量"group"指示每个样本属于哪一组。
2. 差异表达分析流程
2.1 数据准备与预处理
首先，我们需要创建一个DGEList对象，这是edgeR处理数据的基本结构：
# 创建DGEList对象 dge <- DGEList(counts = count_data, group = group) # 过滤低表达基因（CPM>1在至少两个样本中） keep <- rowSums(cpm(dge) > 1) >= 2 dge <- dge[keep, , keep.lib.sizes = FALSE] # 计算标准化因子（TMM方法） dge <- calcNormFactors(dge)2.2 差异表达分析
接下来，我们进行差异表达分析：
# 设计矩阵 design <- model.matrix(~group) # 估计离散度 dge <- estimateDisp(dge, design) # 拟合广义线性模型 fit <- glmQLFit(dge, design) # 进行似然比检验 qlf <- glmQLFTest(fit, coef = 2) # 获取差异表达结果 results <- topTags(qlf, n = Inf)3. 使用edgeR的plotSmear绘制基础MA图
edgeR提供了专门的plotSmear函数来可视化差异表达结果：
# 确定显著差异表达的基因（FDR<0.05，|logFC|>1） de_genes <- decideTestsDGE(qlf, p.value = 0.05, lfc = 1) # 绘制MA图 plotSmear(qlf, de.tags = rownames(qlf)[as.logical(de_genes)], main = "MA Plot (edgeR)", xlab = "Average logCPM", ylab = "logFC", cex = 0.5) abline(h = c(-1, 0, 1), col = c("blue", "black", "blue"), lty = c(1, 2, 1))提示：plotSmear函数会自动处理logFC和平均表达量的计算，用户只需提供差异分析结果对象。
plotSmear函数的主要参数包括：
de.tags：显著差异表达基因的名称smearWidth：控制点分布的宽度col：点的颜色设置pch：点的形状
4. 使用limma的plotMA绘制高级MA图
limma包也提供了plotMA函数，可以生成更加美观的MA图：
# 准备数据 ma_data <- data.frame( logFC = qlf$table$logFC, logCPM = qlf$table$logCPM, FDR = qlf$table$FDR ) # 标记显著基因 ma_data$status <- "NotSig" ma_data$status[ma_data$FDR < 0.05 & ma_data$logFC > 1] <- "Up" ma_data$status[ma_data$FDR < 0.05 & ma_data$logFC < -1] <- "Down" ma_data$status <- factor(ma_data$status, levels = c("Up", "Down", "NotSig")) # 绘制MA图 plotMA(ma_data, status = ma_data$status, main = "MA Plot (limma)", xlab = "Average logCPM", ylab = "logFC", col = c("red", "blue", "gray"), cex = 0.6) abline(h = c(-1, 0, 1), col = c("blue", "black", "red"), lty = c(1, 2, 1))limma的plotMA函数提供了更多自定义选项：
status：基因的分类（上调、下调、不显著）col：不同类别基因的颜色hl.cex：显著基因的点大小xlim/ylim：坐标轴范围
5. 使用ggplot2创建发表级MA图
为了获得更加精美的发表级图形，我们可以使用ggplot2进行自定义绘制：
library(ggplot2) library(ggrepel) # 用于避免标签重叠 # 准备数据 plot_data <- data.frame( gene = rownames(qlf$table), logFC = qlf$table$logFC, logCPM = qlf$table$logCPM, FDR = qlf$table$FDR ) # 标记显著基因和top基因 plot_data$significant <- "No" plot_data$significant[plot_data$FDR < 0.05 & abs(plot_data$logFC) > 1] <- "Yes" # 选择top10差异表达基因进行标注 top_genes <- plot_data[order(plot_data$FDR), ] top_genes <- top_genes[top_genes$significant == "Yes", ] top_genes <- head(top_genes, 10) # 创建ggplot对象 ma_plot <- ggplot(plot_data, aes(x = logCPM, y = logFC, color = significant)) + geom_point(alpha = 0.6, size = 1.5) + scale_color_manual(values = c("gray", "red")) + geom_hline(yintercept = c(-1, 1), linetype = "dashed", color = "blue") + geom_hline(yintercept = 0, color = "black") + geom_text_repel(data = top_genes, aes(label = gene), size = 3, box.padding = 0.5, max.overlaps = Inf) + labs(title = "Differential Expression MA Plot", x = "Average logCPM", y = "log2 Fold Change", color = "Significant") + theme_minimal() + theme(legend.position = "right", plot.title = element_text(hjust = 0.5, face = "bold")) # 显示图形 print(ma_plot) # 保存高分辨率图片 ggsave("MA_plot.png", plot = ma_plot, width = 8, height = 6, dpi = 300)使用ggplot2绘图的优势在于：
- 完全可定制的图形元素（颜色、大小、形状等）
- 轻松添加基因标签而不重叠
- 支持多种输出格式和高分辨率导出
- 可以轻松添加额外的图层和注释
6. 高级定制与问题解决
6.1 颜色与主题定制
在ggplot2中，我们可以轻松修改图形的整体外观：
# 自定义颜色方案 ma_plot + scale_color_manual(values = c("gray", "#1E88E5"), labels = c("Not significant", "FDR < 0.05 & |logFC| > 1")) + theme_bw() + theme(panel.grid.major = element_line(color = "gray90"), panel.grid.minor = element_blank(), legend.title = element_blank())6.2 处理大规模数据集
当处理包含数万个基因的数据集时，绘图可能会变得缓慢。可以考虑以下优化策略：
# 只绘制显著基因和随机子集的不显著基因 set.seed(123) sig_genes <- which(plot_data$significant == "Yes") non_sig_genes <- which(plot_data$significant == "No") sampled_genes <- sample(non_sig_genes, min(5000, length(non_sig_genes))) plot_data_filtered <- plot_data[c(sig_genes, sampled_genes), ] ggplot(plot_data_filtered, aes(x = logCPM, y = logFC, color = significant)) + geom_point(alpha = 0.6) + # 其余参数保持不变6.3 添加额外信息
我们可以在MA图上添加额外的统计信息：
# 计算并显示显著基因数量 num_up <- sum(plot_data$significant == "Yes" & plot_data$logFC > 0) num_down <- sum(plot_data$significant == "Yes" & plot_data$logFC < 0) ma_plot + annotate("text", x = max(plot_data$logCPM)*0.9, y = max(plot_data$logFC)*0.9, label = paste("Up:", num_up, "\nDown:", num_down), hjust = 1, vjust = 1, size = 4)6.4 交互式MA图
对于探索性分析，可以使用plotly创建交互式MA图：
library(plotly) plot_data$text <- paste("Gene:", plot_data$gene, "\nlogFC:", round(plot_data$logFC, 2), "\nFDR:", format(plot_data$FDR, scientific = TRUE, digits = 2)) p <- ggplot(plot_data, aes(x = logCPM, y = logFC, color = significant, text = text)) + geom_point(alpha = 0.6) + scale_color_manual(values = c("gray", "red")) + labs(x = "Average logCPM", y = "log2 Fold Change") ggplotly(p, tooltip = "text") %>% layout(hoverlabel = list(bgcolor = "white"))7. 输出与报告整合
7.1 导出高质量图片
对于发表级图片，需要注意以下参数：
# PDF格式（矢量图，适合编辑） ggsave("MA_plot.pdf", plot = ma_plot, width = 8, height = 6, device = cairo_pdf) # 使用cairo PDF设备确保字体嵌入 # TIFF格式（高分辨率位图，适合期刊投稿） ggsave("MA_plot.tiff", plot = ma_plot, width = 8, height = 6, dpi = 600, compression = "lzw")7.2 在R Markdown中整合分析
可以将整个分析流程整合到R Markdown文档中，实现可重复研究：
```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE, fig.width = 8, fig.height = 6) ``` ## RNA-seq差异表达分析 ### 数据预处理 ```{r preprocessing} # 数据预处理代码 dge <- DGEList(counts = count_data, group = group) keep <- rowSums(cpm(dge) > 1) >= 2 dge <- dge[keep, , keep.lib.sizes = FALSE] dge <- calcNormFactors(dge) ``` ### 差异表达分析 ```{r diff_exp} # 差异分析代码 design <- model.matrix(~group) dge <- estimateDisp(dge, design) fit <- glmQLFit(dge, design) qlf <- glmQLFTest(fit, coef = 2) ``` ### MA图可视化 ```{r ma_plot, fig.cap="MA plot showing differentially expressed genes"} # MA图代码 ggplot(plot_data, aes(x = logCPM, y = logFC, color = significant)) + geom_point(alpha = 0.6) + scale_color_manual(values = c("gray", "red")) + labs(x = "Average logCPM", y = "log2 Fold Change") ```7.3 结果解读与报告
在MA图中，我们通常关注：
- 对称性：好的MA图应该大致对称于logFC=0线
- 趋势线：可以使用局部回归(loess)线检查是否存在表达量依赖的偏差
- 显著基因分布：检查上下调基因是否均匀分布
可以添加趋势线来评估数据质量：
ma_plot + geom_smooth(method = "loess", se = FALSE, color = "darkgreen", size = 0.5)在实际项目中，我经常发现初学者容易忽略MA图的对称性检查。一个常见的错误是只关注显著基因的数量，而忽视了整体分布模式。通过添加趋势线，可以快速识别是否存在低表达基因的假阳性问题，这在数据分析质量评估中非常有用。
