作为天天跟转录组数据打交道的人我对一个场景特别熟悉跑完DESeq2拿到一个几万行的results表格然后对着log2FoldChange和padj两列发呆。表格确实在那儿基因也确实筛出来了但你要是直接把这张表扔给导师或者写进文章里大概率会被问一句“你的差异基因到底长什么样整体趋势是什么”这时候火山图和热图就是最直观的答案。但这俩图恰恰是很多初学者卡壳的地方——RStudio里报错一堆调出来的图要么配色辣眼睛要么聚类乱成一团要么导出之后分辨率低得没法看。这篇文章就是来解决这个问题的。我会从读懂DESeq2结果表开始把火山图和热图的原理讲清楚再给一个可以直接跑的R脚本最后把我在实际使用中踩过的坑一并列出来。目标是让你照着做5分钟出一套能放进文章里的图。1. 先搞懂DESeq2结果表再谈画图——每一列到底在说什么很多人的误区是一上来就画图但图只是结果的呈现你要是看不懂结果表画出来的图也就只能自欺欺人。DESeq2跑完之后用results()函数提取出来的表格每一行是一个基因每一列是一个统计量。我在教学和实际项目中见过太多人只看log2FoldChange和padj两列其他列一概不管这其实是个坏习惯。1.1 结果表的六列数据分别代表什么标准的结果表长这样列名含义怎么理解baseMean该基因在所有样本中的平均归一化计数衡量基因总体表达量高低数值太低说明这个基因本身就不怎么表达差异分析的可信度要打问号log2FoldChange组间表达差异的对数倍数变化正数代表在处理组/实验组中上调负数代表下调绝对值越大差异幅度越大lfcSElog2FoldChange的标准误反映log2FC估计的稳定性这个值越大说明这个基因的差异越不可靠statWald统计量等于log2FC除以lfcSE是计算p值的中间量pvalue未校正的P值单基因层面的显著性做多重检验校正之前的值padj校正后的P值推荐使用的显著性指标DESeq2默认用BH方法校正也就是FDR我建议你形成肌肉记忆判断差异基因看padj和log2FoldChange判断数据质量看baseMean和lfcSE。一个基因就算log2FC算出来是10如果baseMean只有0.5lfcSE大到3以上那这个“差异”很可能是低表达量基因的噪声导致的筛选时要格外小心。1.2 padj和pvalue怎么选——为什么我不建议用原始p值这个问题在答疑群里几乎每周都会被问一次。标准答案是用padj。差异表达分析是几万次统计检验同时进行你要是不做多重检验校正哪怕所有基因都没有真实差异也会有一堆假阳性冒出来——p 0.05这个阈值在几千次检验里光靠随机误差就能刷出几百个“显著”基因。DESeq2默认的校正方法是Benjamini-Hochberg也就是控制FDR。你只需要知道一个核心逻辑padj的意思是“在所有被我判定为显著的基因里预期有多少比例是假阳性”。所以你用padj 0.05做阈值意思是允许5%的假阳性率这在转录组领域是大家普遍接受的默认标准。在提交文章时审稿人通常也默认你用的是padj如果你用原始p值材料方法部分就得专门说明理由否则容易被质疑。1.3 差异基因筛选的经典阈值组合筛选差异基因最常用的阈值组合是|log2FoldChange| 1对应上下调2倍。取log2是为了让上调2倍log2FC 1和下调2倍log2FC -1在数值上对称而且倍数变化在表达量层面是乘法关系取对数后才适合做线性统计建模。padj 0.055%的FDR阈值。baseMean 10过滤掉极低表达的基因。这一步不是DESeq2必需的但很多文章会写主要是为了避免低表达量基因的倍数变化虚高。满足条件的基因就是你的差异表达基因DEGs后面所有图、所有富集分析、所有后续实验的基因清单都是从这一批基因里来的。1.4 DESeq2的独立过滤机制——为什么会有基因padj显示为NA新手第一次看结果表往往会慌怎么这么多行padj是NA是不是算错了不是。这是DESeq2的一种主动过滤机制。它会把baseMean低到几乎没有统计检验能力的基因自动把padj标记为NA不参与多重检验校正。原理是这些基因本身表达量极低就算算出p值也基本不可能显著把它们排除掉可以降低多重检验的总次数提高剩余基因的检验效能。所以看到NA千万别以为是bug也不用特意去补全。分析流程里处理这个很简单筛选差异基因时先用!is.na(res$padj)把NA行剔除再做阈值判断就行。脚本里我会写清楚。2. 火山图横轴纵轴怎么选、阈值怎么定、颜色怎么分火山图的大名估计没人不知道但很多人的理解停留在“左边红、右边红、中间灰”的层面。你要是想画出真正能发表、能让审稿人一眼看懂的水平得先搞清楚它背后的数学逻辑。2.1 火山图各轴的含义与坐标计算火山图本质上是一张散点图横轴log2FoldChange。0的位置代表没有差异右边是上调左边是下调。横轴范围过大的情况通常是某些低表达基因在组间出现或消失导致log2FC极端绘图时可用xlim控制范围。纵轴-log10(padj)也就是把padj取负对数。为什么这么干因为padj是0到1之间的数越接近0越显著在坐标轴上会全部挤在底部。取-log10之后padj越小-log10(padj)越大显著基因就被“顶”到图的上方。padj 0.05时-log10(0.05) ≈ 1.3这就是图上常见水平虚线的位置。整张图看起来像火山就是因为在横轴两端和纵轴上方显著基因密集分布形似火山口喷发的样子。2.2 上、下调和无差异的判定逻辑我画图时习惯把基因分成三类类别条件图上表现显著上调log2FoldChange 1 且 padj 0.05右上区域红色显著下调log2FoldChange -1 且 padj 0.05左上区域蓝色无显著差异其他情况中间灰色区域这个配色方案红/蓝/灰几乎成了转录组文章的事实标准红色代表在外条件刺激下或疾病组中表达升高蓝色代表下降灰色代表没有统计学差异。2.3 为什么有的文章用pvalue而不是padj画火山图这个问题很容易在查阅文献时遇到明明我说用padj怎么有些高分文章纵轴写的是-log10(pvalue)原因是这样Early access时代的很多经典转录组文章确实用pvalue因为pvalue没有经过校正显著基因数量更多图看起来更“热闹”适合展示宽泛的表达变化趋势。但是现在的主流趋势包括大多数期刊审稿人的预期已经转向padj。你要是投的期刊没特别说明建议还是用padj审稿人不容易挑毛病。我自己的习惯是正文里放padj的火山图但如果补充材料里需要展示“更灵敏”的变化趋势可以附一张pvalue版本的图两者不冲突。3. 5分钟出图可直接跑的火山图热图完整R脚本下面这个脚本是我从实际项目中整理出来的删掉了跟项目相关的业务代码只保留通用逻辑。你只需要改文件路径、样本分组、比较组名称就能直接用。3.1 脚本结构概览脚本整体分五步加载数据 → 构建DESeq2对象 → 提取差异结果 → 画火山图 → 画热图。我用的是内置的airway数据集做演示这样你不需要准备任何真实数据就能把全流程跑通看到效果后再替换成自己的数据排错成本最低。3.2 火山图R代码含详细注释# 加载必要的R包 library(DESeq2) library(ggplot2) library(dplyr) library(pheatmap) # 示例数据airway 数据集来自 RNA-Seq 公共数据 # 替换成自己的数据时只需要保证countData和colData格式一致即可 library(airway) data(airway) se - airway # 构建 DESeq2 分析对象 dds - DESeqDataSet(se, design ~ cell dex) # 这一步是核心分析跑完就得到了所有基因的差异统计结果 dds - DESeq(dds) # 提取差异结果contrast指定比较哪个分组 # 这里比较的是 dex 处理组trt相比未处理组untrt res - results(dds, contrast c(dex, trt, untrt)) res - as.data.frame(res) # 先过滤掉 padj 为 NA 的基因原因前面讲过 res - res[!is.na(res$padj), ] # 加一列标记基因是上调、下调还是不显著 res$change - NS res$change[res$log2FoldChange 1 res$padj 0.05] - UP res$change[res$log2FoldChange -1 res$padj 0.05] - DOWN # 统计一下三类基因的数量后面标到图里用 table(res$change) # 定义火山图的配色红色上调、蓝色下调、灰色不显著 my_colors - c(UP #D32F2F, DOWN #1976D2, NS #BDBDBD) # 开始画火山图 ggplot(res, aes(x log2FoldChange, y -log10(padj), color change)) geom_point(size 1.2, alpha 0.6) scale_color_manual(values my_colors) # 画横轴和纵轴的阈值虚线位置就是log2FC±1和padj0.05 geom_vline(xintercept c(-1, 1), linetype dashed, color #444444, linewidth 0.5) geom_hline(yintercept -log10(0.05), linetype dashed, color #444444, linewidth 0.5) labs( x expression(log[2]~Fold~Change), y expression(-log[10]~adjusted~P~value) ) theme_classic(base_size 14) theme( legend.position top, legend.title element_blank(), axis.title element_text(size 14), axis.text element_text(size 12) )这段代码跑完之后你会在RStudio的Plots窗口看到一张基本的火山图。先别急着导出我后面会专门讲如何调到“发表级”。3.3 热图R代码含详细注释热图部分的思路和火山图不一样火山图展示的是所有基因的全局差异分布热图则只需要展示筛选出的差异基因在样本间的具体表达模式。我通常是取top显著基因来画这样图上的基因名标注清晰聚类块也一目了然。# 按显著性排序取前30个差异最显著的基因 sig_genes - res %% arrange(padj) %% head(30) # 提取这30个基因的归一化计数矩阵 # 注意这里用vst方差稳定变换后的数据做热图比直接用counts更合理 # 因为counts的均值-方差关系会导致高表达基因的数值范围过大把热图的颜色梯度拉偏 vsd - vst(dds, blind TRUE) mat - assay(vsd)[rownames(vsd) %in% rownames(sig_genes), ] # 做行基因方向的Z-score归一化 # 把每个基因在所有样本中的表达量转换为以0为中心的标准分数 # 这样处理的目的是不同基因的基础表达量差异很大有的几百有的几万 # 不归一化的话高表达基因会把整个颜色映射主导掉低表达基因的模式根本看不出来 mat_scale - t(scale(t(mat))) # 画热图 pheatmap( mat_scale, cluster_cols TRUE, cluster_rows TRUE, show_rownames TRUE, show_colnames TRUE, fontsize_row 9, fontsize_col 10, color colorRampPalette(c(#1976D2, white, #D32F2F))(100), border_color NA, main Top 30 Significant DEGs )3.4 脚本中几个容易被忽略但很关键的设计第一点热图数据用vst而不是原始counts。这是我反复强调的。原始counts数据不同基因的表达量跨度太大动不动就是几十到几万倍的差距如果直接用原始值做热图低表达基因的颜色会被高表达基因“吞掉”。vst的作用就是把数据做方差稳定变换让不同表达量水平的基因在数值上可比。第二点Z-score归一化是按行基因进行的。scale(t(mat))的t转置你可能第一次看会觉得绕实际逻辑是矩阵原本是“行基因列样本”但scale()函数是按列操作的所以要先把矩阵转置成“行样本列基因”scale完再转置回来。归一化后每个基因在样本间的均值为0标准差为1颜色红蓝代表相对高低而不是绝对表达量。第三点盲法变换blind TRUE。vst的blind参数设置为TRUE意味着变换过程中不考虑样本的分组信息。这对热图展示是公平的因为如果变换过程利用了分组信息得到的表达矩阵会自带“偏向”后续聚类的结果就会有循环论证的嫌疑。4. 脚本里没写但实际必踩的坑——从安装到报错全排查脚本能跑通是理想情况实际过程中我见过太多人在前三步就被卡住。这一节我把高频问题按出现顺序整理一遍你遇到报错时对照着查就行。4.1 环境安装BiocManager装不上DESeq2怎么办DESeq2是Bioconductor包不能用install.packages()直接装。正确姿势是if (!requireNamespace(BiocManager, quietly TRUE)) install.packages(BiocManager) BiocManager::install(DESeq2)常见的坑有两个R版本太老。DESeq2的更新频率不算特别高但新版本的Bioconductor往往会要求新的R版本。我以前在服务器上遇到过R 3.5时代装的DESeq2后来想更新结果直接报错说不兼容。解决办法是去 R官网镜像 装新版本的R或者用conda建一个独立的R环境别跟系统R混在一起。网速问题导致下载中断。BiocManager默认从Bioconductor官方服务器下载国内访问偶尔会超时。如果你遇到Timeout of 60 seconds was reached之类的报错可以设置镜像options(BioC_mirror https://mirrors.tuna.tsinghua.edu.cn/bioconductor) options(CRAN https://mirrors.tuna.tsinghua.edu.cn/CRAN) BiocManager::install(DESeq2)装好之后验证一下library(DESeq2)不报错就说明环境没问题。4.2 DESeqDataSetFromMatrix和DESeqDataSet的区别——何时用哪个我的示例代码用的是DESeqDataSet因为airway数据集是SummarizedExperiment对象可以直接用。但实际项目中你手里的数据大概率是“一个counts表格 一个样本信息表”这时候要用DESeqDataSetFromMatrixcount_matrix - read.csv(counts.csv, row.names 1, check.names FALSE) col_data - read.csv(colData.csv, row.names 1) dds - DESeqDataSetFromMatrix( countData count_matrix, colData col_data, design ~ condition )这里最容易犯的错误是counts表格的格式行名必须是基因名列名必须和colData的行名严格一一对应。我之前接过一个求助对方报错说“nrow of Coldata does not equal nrow of countData”后来发现是列名没对齐——counts表里有一列样本在colData里找不到或者两个文件的样本名大小写不一致。DESeq2在这块非常严格容不得半点马虎。4.3 一个让无数人崩溃的报错all samples have the same condition如果你的设计公式是~ condition但condition只有一种取值——比如所有样本都标记成“control”——DESeq2会直接在DESeq()这一步报错。这个道理很直白没有组间差异可比较自然无法做差异分析。这个报错的常见场景是样本信息表读入后condition列被当成数值型0和1没有转成因子。DESeq2设计公式要求分组变量是因子类型。解决办法很简单col_data$condition - factor(col_data$condition, levels c(control, treatment))levels的顺序决定了比较的基准组第一个level是分母后面的是分子。如果你发现跑出来的log2FoldChange正负方向和预期相反就是levels顺序反了。4.4 火山图坐标轴范围太宽或太窄怎么处理实际数据里经常遇到一个情况部分基因的log2FC高达10甚至15或者padj小到1e-300导致整张图的点全部挤在中间看起来很“空”。这纯粹是坐标轴范围的问题。解决办法是加上坐标轴限制ggplot(...) xlim(-8, 8) ylim(0, 80)但注意xlim和ylim会直接裁剪数据超出范围的基因会被丢到图外。这对火山图的展示影响不大主要差异基因通常落在坐标范围内但如果你担心审稿人问“那些被裁掉的基因去哪了”可以在材料方法里写一句“火山图展示|log2FC| 8的基因超出范围的数据点不计入作图”。另一种更优雅的方案是用ggrepel包给最显著的基因加标签同时用坐标轴限制让图更紧凑library(ggrepel) # 先筛选出要标注的基因padj最小的前10个 top_genes - res %% arrange(padj) %% head(10) ggplot(res, aes(x log2FoldChange, y -log10(padj), color change)) geom_point(size 1.2, alpha 0.6) geom_text_repel( data top_genes, aes(label rownames(top_genes)), size 3, max.overlaps 20 )加了标签之后图和文章的“讲故事”能力会强很多读者一眼就能看到关键的基因是哪几个。4.5 热图的聚类结果乱成一团怎么办热图里最让人头疼的问题就是聚类结果看起来毫无规律样本没有按照分组聚在一起基因也是东一块西一块。这个问题的根源往往是数据质量问题不是画图参数的问题。我排查的顺序是先看样本是否按照处理/对照分组聚在一起。如果样本没有按组聚类先检查是不是样本标签搞反了再检查批次效应——你从公共数据库下载的数据尤其容易出现这种情况不同批次测序的样本混在一起组间差异很容易被批次差异掩盖。基因聚类乱的情况通常是差异基因本身没有协调一致的表达模式。这种情况在生物学上可能也合理——有些基因确实只在某个时间点或特定条件下协同变化。画图时可以用cluster_rows FALSE先取消基因聚类按p值排序从上往下排列出来的图也很有说服力很多人反而更喜欢这种干干净净的表达矩阵展示。4.6 中文显示成方块怎么解决R语言在Linux服务器上画图中文经常显示成方框框。尤其你的基因名是中文注释比如某些非模式物种用中文描述或者图的标题里有中文导出PDF后中文字体全部丢失。两个办法换字体。在画图前设置字体为系统中文字体# Windows windowsFonts(Arial windowsFont(Arial)) # macOS/Linux 可以用 showtext 包来加载系统字体 library(showtext) font_add(Arial, arial.ttf) showtext_auto()最省事的办法图中不要出现中文。发表级图片通常要求英文基因名比如TP53本来就是英文坐标轴和标题也用英文完全没必要在图上写中文。这也是行业惯例我基本不在图里放中文从根源上规避了这个坑。4.7 导出图片清晰度不够——300 dpi只是底线RStudio的Plots窗口直接Export保存的图片默认分辨率往往是屏幕级别拉到word里一放大就糊。发表级图片至少需要300 dpi而且最好是矢量图。我的导出标准姿势是ggsave( volcano_plot.png, plot last_plot(), width 8, height 6, dpi 300 )如果你投的期刊要求矢量图绝大多数期刊的图片都要求矢量图用PDF格式ggsave( volcano_plot.pdf, plot last_plot(), width 8, height 6, useDingbats FALSE )PDF格式的好处是无限放大不模糊而且后期用AIAdobe Illustrator编辑非常方便可以改字体、改颜色、加标注。很多杂志的最终出版图都是在PDF基础上微调的。pheatmap的导出稍微特殊一点要用pdf()或者png()设备包裹pdf(heatmap.pdf, width 8, height 8) pheatmap(mat_scale, ...) dev.off()5. 从“能出图”到“能发表”——最终细节打磨清单最后这一部分我把从“能跑出图”到“图能进文章”的细节差异一次性说透。这些细节单独看都不起眼但合在一起就是“论文配图”和“练习作业”的分水岭。5.1 配色不止是好看——色盲友好设计值得了解传统的红绿配色上调红色、下调绿色在色盲人群中容易撞色。我在一次学术会议上听报告时有听众当场提了这个问题——做生物医学研究的谁都不想自己的图被部分读者错误解读。所以现在很多期刊也在逐步推荐色盲友好的配色方案。红蓝配色红色Vs蓝色就是其中一种对红绿色盲人群也相对友好。我的脚本里用的#D32F2F深红和#1976D2深蓝就是从Material Design色板里选的在屏幕上和打印出来都清晰可辨。如果你还需要其他颜色可以参考ColorBrewer的Set1、Set2方案或者用RColorBrewer包的brewer.pal()函数直接取色。5.2 字号和线条粗细——期刊审稿人最在意但从不写在意见里的东西我编辑过几篇论文的配图最大的感受是很多图不是内容不行而是细节不到位。具体来说就是字体太小、线条太细、图例位置遮挡数据点。发表级图片的硬性标准图片最小字号不小于7pt通常要求8pt及以上轴标题字号要大于轴刻度字号主次分明散点大小1.5左右透明度和点大小要配合让密集区域的点不至于重叠到看不出密度阈值虚线用虚线线型线宽0.5到1之间既清晰又不抢数据点的视觉权重图例放在图的右上角或顶部但不能遮挡主要数据区域你需要记住一个原则你在屏幕上看到的图和打印出来的效果不完全一样。屏幕看起来舒服的位置打印出来可能偏小或偏暗。所以导出的PDF最后要放大到100%检查一遍再用AI打开改一下细节才算完成。5.3 图注怎么和图片配合——审稿人真正在看什么这个点很微妙。很多人以为图注是随便写写其实审稿人看图的顺序往往是先看图再看图注最后才看正文。图注写得清不清楚直接决定审稿人能不能快速理解你的核心发现。火山图的图注至少要包含以下信息用什么软件和包做的差异分析DESeq2比较的是哪两个组dex处理组 vs 未处理组阈值标准|log2FC| 1padj 0.05每个颜色代表什么红色上调、蓝色下调、灰色不显著一共鉴定出多少个差异基因上调X个下调Y个热图的图注需要额外说明数据经过什么变换vst归一化颜色代表的是什么Z-score归一化后的相对表达量row和column的聚类依据欧氏距离 层次聚类这些信息不用全部写进图里但必须出现在图注里。审稿人拿着你的图对照图注能完整还原你的分析流程这才是合格的文章配图。5.4 两个实用小技巧scale颜色自定义和批量出图最后分享两个实用技巧。第一个技巧热图的颜色梯度可以用自定义的色板。我的默认配色是蓝色到白色到红色但不同期刊的审美风格不一样。如果你想调整成绿色到黑色到红色的风格可以这样my_breaks - seq(-2, 2, length.out 101) my_colors - colorRampPalette(c(#2C7BB6, #FFFFFF, #D7191C))(100) pheatmap( mat_scale, breaks my_breaks, color my_colors, ... )第二个技巧批量出图。如果你有多个比较组比如“处理A vs 对照”、“处理B vs 对照”、“处理C vs 对照”手动一个一个跑太浪费时间。用一个循环就能搞定treatments - c(A, B, C) for (trt in treatments) { res_i - results(dds, contrast c(treatment, trt, control)) res_i - as.data.frame(res_i) res_i - res_i[!is.na(res_i$padj), ] res_i$change - NS res_i$change[res_i$log2FoldChange 1 res_i$padj 0.05] - UP res_i$change[res_i$log2FoldChange -1 res_i$padj 0.05] - DOWN p - ggplot(res_i, aes(x log2FoldChange, y -log10(padj), color change)) geom_point(size 1.2, alpha 0.6) scale_color_manual(values my_colors) geom_vline(xintercept c(-1, 1), linetype dashed) geom_hline(yintercept -log10(0.05), linetype dashed) labs( x expression(log[2]~Fold~Change), y expression(-log[10]~adjusted~P~value), title paste0(Treatment , trt, vs Control) ) theme_classic(base_size 14) ggsave(paste0(volcano_, trt, .png), p, width 8, height 6, dpi 300) }这套流程跑下来几组比较的图一次全出文件名清晰后续整理还省时间。5.5 一个关于pheatmap出图的提醒先画图到屏幕再导出pheatmap有个小特点它不经过ggplot的对象系统所以不能像ggsave那样直接保存。很多初学者直接在RStudio里用Export保存得到的PNG分辨率不够。正确的做法是用pdf()或者png()包裹我前面已经写过了。另外一个和热图相关的细节如果你的基因名比较长在热图里会重叠显示看起来一团黑。这时候有两个办法一是用fontsize_row 7压缩字号二是把show_rownames改为只显示部分标签比如每隔5个显示一个。但说实话真正发表级的热图基因名不一定要全部标注——你可以只标注你重点关注的几个基因其余不标这样图面干净得多。审稿人真要看所有基因名补充表格里有。我和DESeq2打交道五六年从最初对着报错信息发懵到现在一套脚本跑完差异分析和可视化最深的体会是这类工具用多了你会发现真正拉开差距的从来不是代码本身而是你对数据、对统计原理、对图表表达逻辑的理解。脚本5分钟就能跑完但能跑出什么东西、跑出的东西能不能说服审稿人靠的是前期的积累和细节的把控。这篇里的脚本和建议都是我在实际项目中验证过的。你要是按着操作一遍应该能在半小时内跑通全流程——如果中间有报错翻翻我这篇里的排查部分大概率能找到答案。真要还有搞不定的地方多半是你的数据格式特殊到时候对照DESeq2的官方文档或者differential expression的教程逐行排查即可。