ARTICLE · INTELLIGENCE

战地情报 · 详情页

来自尧图项目组的一线实战观察与深度解析

TCGA突变数据可视化:maftools从安装到瀑布图全流程

TCGA突变数据可视化:maftools从安装到瀑布图全流程 做肿瘤生信分析的朋友应该都有过这种体会从TCGA下载完突变数据满心期待想快速看一眼样本的整体突变情况结果光是把数据整理成能画图的格式就耗掉半天。早些年我还没有接触maftools的时候画一张瀑布图oncoplot需要自己写一长串R代码去统计每个基因的突变频率、按样本排列突变位点、手动调颜色改一次配色就要重跑一遍效率非常低。后来在一次生信培训上接触到maftools这个R包试完第一感受就是这类可视化问题终于有了一个开箱即用的标准答案。这篇内容我会以一个拿到TCGA癌症突变MAF文件的普通分析者为视角从包的安装、数据格式的认识、完整可视化流程到常见坑的排查一步步写清楚。无论你是医学生、科研助理还是刚转行做生信分析的R语言初学者只要跟着下面的代码走5分钟内输出一张可发表的瀑布图完全可行。所有代码我都测试过放在这里可以直接复制运行。1. 为什么突变数据可视化首选maftools1.1 MAF格式TCGA突变数据的通用语言在展开代码之前有必要先搞清楚maftools处理的数据到底长什么样。MAF全称是Mutation Annotation Format这是TCGA项目以及全球多个癌症基因组计划普遍采用的突变数据标准格式。你可以把它理解成一张巨大的表格每一行代表肿瘤基因组中检测到的一个突变事件每一列描述这个突变的一个属性。常见的关键列包括Hugo_Symbol基因名、Chromosome染色体号、Start_Position和End_Position突变起止位置、Variant_Classification突变类型比如Missense_Mutation错义突变、Nonsense_Mutation无义突变、Frame_Shift_Del移码缺失等以及Tumor_Sample_Barcode来自哪个样本。这种格式的巧妙之处在于它把复杂的突变信息压缩成了统一的、机器可读的表格结构。无论你是从GDC Data Portal下载的VarScan2或MuTect2流程产出还是从cBioPortal导出的精简数据最后落到本地基本都会是这个结构。maftools之所以能成为处理这类数据的霸主级工具本质上就是围绕MAF格式做了非常高密度的功能封装把上游数据清洗和下游统计可视化之间的链路打通了。1.2 maftools的设计思路从一个对象开始我第一次使用maftools时印象最深的是它的对象化设计。它不像很多R包那样要你反复传递数据框参数而是通过read.maf()把MAF文件读进来封装成一个MAF对象。后续几乎所有函数无论是画瀑布图、棒棒糖图、做共突变分析还是生存分析都直接吃这个MAF对象作为第一个参数。这个设计语言跟tidyverse里的ggplot2有些神似——先建立“数据层”再在“数据层”上叠加各种图形和分析任务。另外一个关键的设计思路是“默认值合理细节可覆盖”。maftools的函数参数非常多但大部分时候你只需要给两个参数就能跑出一张很专业的图。比如画瀑布图时它默认会按突变频率排序基因自动用不同类型的颜色区分突变种类还会在底部计算每个样本的TMB肿瘤突变负荷。这对非可视化专业出身的研究人员极其友好不需要懂得颜色搭配和排版美学默认输出的结果就已经很能打了。如果非要跟其他方案对比我自己的体会是用ggplot2手搓瀑布图适合做高度定制化的收尾阶段用maftools做90%的常规探索性分析才是效率最大化的路径。它整合了TCGA队列比较、驱动基因识别、共突变分析等功能这些东西如果全部自行实现工作量非常可观。2. 环境准备和MAF文件获取2.1 一行命令完成安装maftools位于Bioconductor仓库所以安装方式和CRAN上的普通R包略有不同。如果你之前装过BiocManager这个辅助包整个安装过程只需要一条命令。以下是我在R 4.x版本下反复验证过的安装流程直接在R控制台或RStudio中逐行运行即可if (!requireNamespace(BiocManager, quietly TRUE)) install.packages(BiocManager) BiocManager::install(maftools)安装完成后用library(maftools)加载。如果控制台没有任何报错说明环境已经就绪。需要注意的是maftools对R版本有一定要求如果你还在用R 3.6或更老的版本建议先升级R环境否则很可能碰到依赖包版本不兼容的问题。关于安装过程中的具体报错我会在第4章展开讲。2.2 从哪获取TCGA的MAF文件从我自己的使用经验来看获取TCGA突变数据最权威且最稳定的渠道是官方数据下载平台GDC Data Portal。在网站上选择你感兴趣的癌种比如肺腺癌LUAD、乳腺癌BRCA或结直肠癌COAD在Data Category中勾选Simple Nucleotide Variation在Data Format中勾选MAF然后下载得到的就是标准MAF文件。这里的文件通常比较大某些癌种的MAF文件动辄上百MB下载后不需要预处理可以直接read.maf读入。除了GDC官网还有一些备选渠道。cBioPortal提供在线导出功能可以选择某个癌症研究项目在Download面板中找到Mutation Data导出的文件同样是MAF格式但列名可能与GDC官方版略有差异读取时需要留意自动识别是否成功。UCSC Xena也提供TCGA的突变数据下载格式相对更精简一些适合后续自定义分析。如果你只是想快速测试maftools的功能不想下载整个癌种的数据可以像我一样直接用包自带的TCGA急性髓细胞白血病LAML示例数据。这个数据集虽然样本量不算大但麻雀虽小五脏俱全瀑布图、生存分析、临床注释都齐备非常适合用来跑通流程。这就是我们下面实操部分用的数据。2.3 read.maf的核心参数别在读取这一步掉链子read.maf()是maftools的门面函数大部分用户的使用场景就是一条命令读完整个文件。但我建议在读取大文件之前先搞清楚几个关键参数否则等到分析快做完才发现数据读错了返工成本很高。第一个需要注意的参数是clinicalData。TCGA的MAF文件虽然包含样本的突变信息但不包含生存时间、年龄、分期等临床信息。如果你后续要做生存分析或者按临床特征分组比较突变谱就需要在这里传入一个临床信息数据框行名是样本条码列是临床变量。maftools通过样本条码自动匹配突变和临床数据所以务必保证MAF文件中的Tumor_Sample_Barcode和临床数据行名对标一致性。第二个是isTCGA参数。当设置为TRUE时maftools会从样本条码中自动提取TCGA样本类型信息比如01代表原发肿瘤、06代表转移灶等并在后续分析中默认排除非肿瘤样本。如果你下载的是TCGA的MAF建议直接设成TRUE省去手动过滤的步骤。第三个是vc_nonSyn这是指定哪些突变类型算作非同义突变。默认情况下maftools已经内置了一套适用于绝大多数场景的列表。但如果你用的是其他项目的数据突变类型的命名风格可能略有区别这时候手动指定vc_nonSyn参数可以避免漏算或误算。还有两个小参数也值得留意。verbose默认是TRUE控制台会输出一堆运行日志初学阶段建议保留方便观察maftools内部做了哪些步骤实际上手熟练之后可以设为FALSE减少信息干扰。removeDuplicatedVariants会在读取时自动去除重复突变记录对于部分整理不规范的数据源这是很有用的保险措施。3. 5分钟核心实操全程代码与图形解读3.1 读取数据建立MAF对象现在我们直接进入正题。打开RStudio新建一个R脚本把下面代码粘贴进去。先加载maftools包然后用它内置的LAML示例数据来跑通整个流程。这里读取MAF文件的方式用的是system.file函数它返回的是包安装路径下某个文件的完整路径设计了一些新手在这里直接复制“extdata/tcga_laml.maf.gz”去读结果提示找不到文件原因就是少了前面那一大段绝对路径。library(maftools) # 读入maftools包自带的TCGA LAML示例MAF文件 laml - read.maf(maf system.file(extdata, tcga_laml.maf.gz, package maftools))如果一切顺利控制台会打印出一大串摘要信息包括样本数量、突变总数、每个样本的中位突变数等。这时候你只需要输入laml并回车maftools就会把整个MAF对象的概览信息再次汇总列出来包括突变类型分布、转换颠换比例、样本间突变负荷排名等。能看到这一页输出说明你的环境已经没问题了接下来画图顺手很多。3.2 绘制整体突变摘要图拿到MAF对象后的第一个标准动作是画一张“突变全景摘要图”。maftools提供了一个叫plotmafSummary的函数它在同一个画布上完成了两件事上半部分是突变类型的堆积柱状图展示每个样本中错义突变、无义突变等各类突变的数量下半部分是每个样本的总突变数条形图让样本的突变负荷高低一眼可见。这样一张图放在论文的附件材料里审稿人基本挑不出毛病。下面是我常用的写法。这里的rmOutlier参数表示是否移除突变数特别高的离群样本有时个别样本的突变数量会远高于其他样本导致柱状图压缩到几乎看不见差别设置成TRUE就能让图形分布更合理。addMax参数则在柱状图顶部额外标出最大值的样本名对小样本量的预实验分析特别有用。plotmafSummary( maf laml, rmOutlier TRUE, addMax TRUE )这里我建议第一次运行先不要设置任何参数直接plotmafSummary(laml)看看默认效果理解了默认输出后再逐步叠加参数去微调。很多初学者一上来就抄一大堆参数反而搞不清楚哪个参数控制哪个元素。3.3 瀑布图突变可视化的扛把子瀑布图是整个maftools的招牌功能也是肿瘤突变研究里最高频出现的图表类型。它本质上是一个热图矩阵行代表基因按突变频率从高到低排列列代表样本中间的每个色块代表该基因在该样本中存在某种特定类型的突变颜色对应突变类型。同时图上方还会自动显示每个样本的TMB值图右侧显示每个基因的突变比例。基础版瀑布图只需要一行代码oncoplot(maf laml, top 20)top 20表示只展示突变频率最高的前20个基因。在真实项目中这个数字一般设置在10到30之间比较合适。如果设置成50图形的横向信息量会过大阅读体验反而不好。如果已经有关注的候选基因可以换成genes参数固定展示某些基因比如genes c(DNMT3A, NPM1, FLT3)。waterfall图在发表时往往需要按自己的配色方案调整。这时可以传入colors参数这是一个具名向量形如oncoplot( maf laml, top 20, colors c( Missense_Mutation #4E79A7, Nonsense_Mutation #F28E2B, Frame_Shift_Del #E15759, Frame_Shift_Ins #B07AA1 ) )配色这件事我个人建议优先采用颜色区分度高的色板别用太相近的颜色否则读者根本分不清错义突变和无义突变。maftools默认的颜色其实已经很成熟除非期刊有明确的配色要求否则没必要大改。还有一个很实用的参数是removeNonMutated默认是TRUE。它会自动把没有发生任何突变的样本从图中移除让瀑布图变得更紧凑。如果分析中需要强调存在一部分全野生型样本可以设置成FALSE让这部分样本也显示在图上。我这里再补充一个实战中很有价值的组合在瀑布图下方加展示临床注释的“热轨道”。这是通过clinicalFeatures参数实现的。比如你想同时观察每个样本的FAB分型和突变谱的关系laml.clin - read.table( system.file(extdata, tcga_laml_clinical.txt, package maftools), header TRUE, sep \t ) oncoplot( maf laml, top 20, clinicalFeatures c(FAB_classification), sortByAnnotation TRUE )这段代码需要先手动读取一个临床信息文本文件它同样来自maftools包自带的示例数据。clinicalFeatures接收一个或多个列名maftools会自动在瀑布图底部画出一个对应的分组条。sortByAnnotation TRUE的意思是按临床注释分组排序样本这样一来样本在图上不再是无序排列而是按FAB分型聚在一起很容易看出不同亚型间的突变模式差异。3.4 单基因棒棒糖图看驱动突变的位置瀑布图解决的是“哪些基因经常突变”的问题但如果你想知道某个基因的突变究竟发生在蛋白的哪个功能域瀑布图就无能为力了。这时候要请出棒棒糖图。所谓lollipop图横坐标是蛋白序列位置纵坐标是突变在该位置上出现的样本数每个“棒棒糖”代表一个突变位点。maftools还能自动标注出已知的蛋白结构域位置方便判断突变是否倾向于落在特定功能区域。lollipopPlot( maf laml, gene DNMT3A, showDomainLabel TRUE, cex 0.5 )这段代码画的基因是DNMT3A一个在LAML中高频突变的表观遗传学调控基因。运行后你会看到图上会标注出PWWP结构域和甲基转移酶结构域大概落在蛋白的哪个区段同时用不同颜色标出错义突变和截短突变。别小看这个功能它能在几秒内帮你判断某个基因的突变是否具有“热点分布”特征——如果很多突变都集中在同一个结构域往往提示这个结构域的功能异常是疾病发生的关键。3.5 森林图与共突变分析两个基因的关系斯不容忽视在肿瘤基因组探索中除了单基因频率大家也很关心基因之间的共发生或互斥关系。比如A基因突变和B基因突变究竟是一起出现还是互相排斥对于理解信号通路之间的代偿机制很有参考价值。maftools提供了somaticInteractions函数它基于配对样本的突变状态做Fisher精确检验输出一张对称矩阵热图红色代表共发生绿色代表互斥。somaticInteractions( maf laml, top 25, pvalue c(0.05, 0.01) )这里的top 25表示取突变频率最高的25个基因做两两组合检验。矩阵图上方会出现一个调节显示阈值的滑块在实际发表截图时可以手动选择合适的显著性展示。我第一次跑这个函数的时候有点困惑为什么热图只显示半三角后来才意识到这个矩阵本身就是对称的展示上三角就够用了全矩阵反而显得冗余。此外还有forestPlot用于绘制基因突变频率的森林图显示突变频率及其置信区间。这个图形的用处更多在于补充材料展示所有基因的突变频率分布情况帮助审稿人快速了解数据整体水平。3.6 生存分析与TCGA队列比较把突变和临床关联起来maftools不只是画图工具它还整合了生存分析功能。比如你想验证某个基因是否具有预后价值可以按该基因是否突变把样本分成两组然后做Kaplan-Meier分析。maftools封装好了计算和绘图一行代码输出一个带有log-rank检验p值的生存曲线图mafSurvival( maf laml, genes DNMT3A, time days_to_last_follow_up, Status Overall_Survival_Status, isTCGA TRUE )这里time和Status分别对应临床数据中的随访时间列和生存状态列。运行后会返回生存曲线的ggplot对象你可以用ggplot2主题函数进一步修改配色、字体等细节。有一点要注意这个函数对临床数据的格式有要求time列必须是数值型Status列必须是0和1的数值编码如果从Excel复制数据时混入了字符串会直接报错。TCGA队列比较是另一个我很喜欢的功能。tcgaCompare函数会把当前数据的TMB和各癌种的TCGA公共队列放在一起做对比在同一个箱线图上展示你的队列处在什么水平。这在肿瘤免疫相关分析中特别常见因为TMB高低直接和免疫治疗效果相关。tcgaCompare(maf laml)运行后应该可以看到LAML样本的TMB分布和泛癌种比较的结果。它会自动下载一个内置参考队列文件如果网络环境一般这一步可能会卡一会儿。如果下载失败也可以从maftools的GitHub仓库手动下载后通过tcgaCompare(maf laml, tcgaCache 路径)指定本地缓存文件这个技巧我后边细说。3.7 药物互作预测和突变签名进阶加分项如果分析做到后面还有精力maftools还有一批比较进阶的功能值得尝试。drugInteractions可以根据已突变基因在Drug Gene Interaction Database中检索可能相关的药物输出一张药物与基因互作的表格。它还能自动对每一个药物列出对应的互作基因数方便快速筛选有潜在临床指导意义的候选药物。drugInteractions(maf laml)这个函数的输出结果是一个数据框需要先赋值给一个对象再用View查看直接运行的话控制台输出会有截断。如果你对肿瘤突变签名mutational signature感兴趣还可以用signatureEnrichment函数估计样本中是否存在COSMIC已知的突变签名并和参考队列做富集比较。不过坦率说这两块功能的参数复杂度和运行时间都明显高于前面的基础可视化对非生信背景的朋友来说前期可以先跑通不必死磕。4. 常见问题与排查技巧实录4.1 安装环节卡住十有八九是依赖包版本问题maftools安装失败的案例我见过太多了最常见的报错是“there is no package called ‘xxx’”其中xxx往往是data.table、RColorBrewer等基础依赖。这种情况的根源通常不是网络问题而是当前环境中这些依赖包的版本过旧和maftools的最新版本不兼容。我自己的解决套路是先执行BiocManager::install(maftools, update TRUE, ask FALSE)强制让BiocManager去检查并更新所有相关依赖如果还是失败再单独安装报错信息里提到的那个包然后重新加载maftools。另一个非常隐蔽的坑是某些服务器环境存在多个R版本切换R环境后旧环境里安装的包在新环境下完全不可见导致library(maftools)直接报错。遇到这种情况先用sessionInfo()确认当前R版本和.libPaths()指向的包安装路径是否正常再决定是重装还是调整环境变量。4.2 数据读入不是所有MAF都叫TCGA标准MAF虽然MAF有官方标准但实际操作中读入其他项目产生的MAF文件时经常会遇到一些“同分异构”问题。最常见的是列名大小写不对比如hugo_symbol而不是Hugo_Symbolmaftools默认是严格匹配大写列名的大小写不一致会导致关键列识别正确但注释列大面积缺失。另一个典型问题是非TCGA数据源的样本条码格式五花八门直接用isTCGA TRUE会引发样本类型推断错误。这时候建议设isTCGA FALSE通过临床数据或样本名自己构建分组。如果你读入的文件报“Required fields are missing”之类的错误建议先不要直接用read.maf而是先用read.table或fread把原始文件读进来用colnames()检查一下表头确认哪些列名和标准MAF不一致再用rename或data.table的setnames修正。这属于笨办法但在面对陌生数据源时老实检查表头的效率远高于一遍遍尝试读入。4.3 画图环节内存、字体与输出格式maftools在读取包含几万行突变的大MAF文件时内存占用会明显偏高。如果你用的是老一点的笔记本电脑在读取前建议先用fread或read.table配合nrows参数抽样预览文件结构确认无误后完整读入时关闭其他占内存的应用。如果行为个别的多等一会儿可以用options(maftools.verbose FALSE)关掉输出日志节省一点I/O时间。字体问题是另一个多次困扰我的细节。在中文Windows系统上某些maftools自动生成的图在导出PDF后中文注释可能出现方块乱码。我的习惯是凡是奔向投稿的图全部通过pdf()或png()函数显式控制输出尺寸和字体尽量不依赖RStudio自带的绘图窗口直接另存为图片。比如pdf(oncoplot.pdf, width 10, height 8) oncoplot(maf laml, top 20) dev.off()这样可以保证图片的像素和排版完全可控避免后期再返工。为了帮你快速排查我把上面讲到的常见问题整理成了一张速查表症状可能原因解决思路安装时报错“no package called”依赖包版本过旧或缺失运行BiocManager::install(maftools, update TRUE, ask FALSE)强制更新依赖读入后样本数明显偏少样本条码格式不标准或isTCGA设置错误设isTCGA FALSE检查Tumor_Sample_Barcode列瀑布图基因名显示重叠top值太大导致行间距过小减小top值或调整fontSize参数生存分析报错time列有非数值临床数据列格式错误用as.numeric()转换时间列确认Status为0/1编码tcgaCompare下载卡住网络问题导致无法下载内置参考文件手动下载参考队列用tcgaCache参数指定本地路径输出图片在Word中模糊直接截图保存使用pdf()或png(res 300)控制导出分辨率4.4 让百色图在报告里更耐看的几个微调画图这种事的最终效果很多时候不是功能不够而是微调不得要领。瀑布图如果顶部样本条码挤成一团把barcode_mar参数调大一些比如barcode_mar 3000能明显缓解标签重叠。如果发现单个基因的突变类型太过单一导致图右侧图例几乎只有两种颜色可以在oncoplot里加drawRowBar FALSE先把行柱状图关掉让整体构图更干净。如果是做组会汇报或者论文初稿我建议输出前把图整体放大到合适尺寸。同一个oncoplot在RStudio窗口里看效果不错导出PNG后往往因为分辨率不够显得模糊。我的习惯是设置png(oncoplot.png, width 3000, height 2200, res 300)这样插入Word或PPT里至少是印刷清晰度。对于方框图、箱线图这些矢量图形则推荐保存为PDF后续在Illustrator里统一排版也不会失真。5. 一些适合继续深挖的扩展方向写完上面的基础流程如果你还有余力可以试着把maftools和其他R生态工具组合使用。比如把突变数据和表达数据联起来做关联分析或者借助complexheatmap包在oncoplot的基础上自定义更复杂的注释热图。不过从投入产出比来看我仍然建议先熟练掌握maftools自带的核心环节等真正理解了MAF数据结构和瀑布图的视觉逻辑后再考虑深度定制。从我这个过来人的角度看maftools最大价值在于非常明显降低了肿瘤突变数据的入门门槛。我以前手写瀑布图的时候一边整理数据一边怀疑人生如今用maftools几分钟就能出图而且代码稳定可复现很多早期的分析脚本到现在还在用。只要你手里有一份标准MAF文件按照上面的代码走一遍基本就能覆盖90%的日常探索需求。最后再说一个我的个人习惯每接到一个新的癌种数据我第一步永远是plotmafSummary oncoplot(top 20)快速建立全局认知看着瀑布图从空荡荡到布满色块整个样本队列的突变画像就迅速立体起来了。这种“先总览后深挖”的工作流比一开始就针对单个基因做细致分析要高效得多。希望这份上手流程能帮你少踩一些坑省下更多时间去做真正有价值的生物学解读。
RELATED READING

延伸阅读

更多一线实战笔记与深度复盘,助您持续精进