本文首发于“生信补给站”公众号 https://mp.weixin.qq.com/s/2HJxx0KsYhEp8MdGB1nNFQ
随着癌症基因组学的进步, 突变注释格式 (MAF) 被广泛用于存储检测到的somatic variants。The Cancer Genome Atlas 项目对30多种不同的癌症进行了测序，每种癌症类型的样本量超过200种。maftools-R包能够有效的汇总，分析和可视化MAF格式的文件。
maftools函数主要分为可视化和分析两个模块，其主要功能及简短的描述如下所示，使用时只需读取MAF文件然后降MAF对象传递给所需要的绘图或分析功能即可。
本次主要使用R-maftools包绘制组学突变结果（MAF）的oncoplot或者叫“瀑布图”。
一、 载入R包，数据
1） 载入maftools包
if (!require("BiocManager")) install.packages("BiocManager") BiocManager::install("maftools")
2） 载入数据
通过read.maf函数读入MAF文件，将各种数据（组学基因突变，拷贝数变异，临床数据，表达数据等）汇总并将其存储为MAF对象。
library(maftools) #TCGA-LAML MAF file (gz) laml.maf = system.file('extdata', 'tcga_laml.maf.gz', package = 'maftools') #clinical information (optional) laml.clin = system.file('extdata', 'tcga_laml_annot.tsv', package = 'maftools') laml = read.maf(maf = laml.maf, clinicalData = laml.clin)
MAF对象中除了上面基因突变数据和对应的临床数据外，还可以加入拷贝数变异，表达数据等其他数据类型，后面需要的时候会添加。
二、 MAF对象汇总
1） 展示MAF重点变量的summary信息
#Shows sample summry. getSampleSummary(laml) #Shows gene summary. getGeneSummary(laml) #Shows all fields in MAF getFields(laml) #shows clinical data associated with samples getClinicalData(laml) #Writes maf summary to an output file with basename laml. write.mafSummary(maf = laml, basename = 'laml')
上图为临床数据（getClinicalData）的summary结果，其余可自行输入查看结果。
2） 绘制MAF-summary图
使用 plotmafSummary 绘制 maf 文件的summary信息，如下：
#plotmafSummary
plotmafSummary(maf=laml, rmOutlier=TRUE, addStat='median', dashboard=TRUE, titvRaw=FALSE)
堆叠的 barplot展示maf 文件中每个样本中的变异数量，并添加中位线，以显示队列间的中位数变异数量。箱线图展示variant_Classification的变异类型。
三、 绘制oncoplot（瀑布）图
1 绘制基础oncoplots（瀑布图）
oncoplots或者瀑布图可以很好的展示maf文件中的变异信息，侧面条形图和顶部条形图可分别由drawRowBar和drawColBar参数控制。
#展示top20的变异genes. oncoplot(maf = laml, top = 20)
注：变异注释为Multi_Hit表示同一样本中突变多次的基因。
2 更改变异类型的颜色
#此处使用RColorBrewer的颜色，当然也可以使用任意颜色 vc_cols = RColorBrewer::brewer.pal(n = 8, name = 'Paired') names(vc_cols) = c( 'Frame_Shift_Del', 'Missense_Mutation', 'Nonsense_Mutation', 'Multi_Hit', 'Frame_Shift_Ins', 'In_Frame_Ins', 'Splice_Site', 'In_Frame_Del' ) #查看变异类型对应的颜色 print(vc_cols) #> Frame_Shift_Del Missense_Mutation Nonsense_Mutation Multi_Hit #> "#A6CEE3" "#1F78B4" "#B2DF8A" "#33A02C" #> Frame_Shift_Ins In_Frame_Ins Splice_Site In_Frame_Del #> "#FB9A99" "#E31A1C" "#FDBF6F" "#FF7F00" oncoplot(maf = laml, colors = vc_cols, top = 20)
3 添加copy number 信息
可以使用两种方式向maf文件中添加SCNA信息
- GISTIC result
- Custom copy number table
3.1 GISTIC results
GISTIC会得到很多结果，此处只需要四个文件即可：all_lesions.conf_XX.txt, amp_genes.conf_XX.txt, del_genes.conf_XX.txt, scores.gistic：
#GISTIC results all.lesions = system.file("extdata", "all_lesions.conf_99.txt", package = "maftools") amp.genes = system.file("extdata", "amp_genes.conf_99.txt", package = "maftools") del.genes = system.file("extdata", "del_genes.conf_99.txt", package = "maftools") scores.gis = system.file("extdata", "scores.gistic", package = "maftools") #Read GISTIC results along with MAF laml.plus.gistic = read.maf( maf = laml.maf, gisticAllLesionsFile = all.lesions, gisticAmpGenesFile = amp.genes, gisticDelGenesFile = del.genes, gisticScoresFile = scores.gis, isTCGA = TRUE, verbose = FALSE, clinicalData = laml.clin ) #绘制含有SCNA信息的oncolpot oncoplot(maf = laml.plus.gistic, top = 20)
3.2 Custom copy-number table
可以自定义一个copy-number结果，例如，让我们在随机的 20 个样本中创建 DNMT3A 的假 SCNA变异。
set.seed(seed = 1234) barcodes = as.character(getSampleSummary(x = laml)[,Tumor_Sample_Barcode]) #Random 20 samples dummy.samples = sample(x = barcodes, size = 20, replace = FALSE) #Genarate random CN status for above samples cn.status = sample( x = c('Amp', 'Del'), size = length(dummy.samples), replace = TRUE ) custom.cn.data = data.frame( Gene = "DNMT3A", Sample_name = dummy.samples, CN = cn.status, stringsAsFactors = FALSE ) head(custom.cn.data) #> Gene Sample_name CN #> 1 DNMT3A TCGA-AB-2898 Amp #> 2 DNMT3A TCGA-AB-2879 Amp #> 3 DNMT3A TCGA-AB-2920 Del #> 4 DNMT3A TCGA-AB-2866 Amp #> 5 DNMT3A TCGA-AB-2892 Amp #> 6 DNMT3A TCGA-AB-2863 Amp #添加SCNA信息 laml.plus.cn = read.maf(maf = laml.maf, cnTable = custom.cn.data, verbose = FALSE) oncoplot(maf = laml.plus.cn, top = 10)
4 添加 significance values
输入data.frame格式的数据，包括gene and q值两列即可：
#MutSig results laml.mutsig = system.file("extdata", "LAML_sig_genes.txt.gz", package = "maftools") oncoplot( maf = laml, mutsig = laml.mutsig, mutsigQval = 0.01, )
5 临床信息 annotations
1）使用存储在 clinical.data 中的变量进行注释 #使用FAB_classification注释 oncoplot(maf = laml, clinicalFeatures = 'FAB_classification')
2）临床信息注释的变量同样可以自定义颜色：
#更改FAB classification的颜色设置 fabcolors = RColorBrewer::brewer.pal(n = 8,name = 'Spectral') #颜色和变量的种类要一样多 names(fabcolors) = c("M0", "M1", "M2", "M3", "M4", "M5", "M6", "M7") fabcolors = list(FAB_classification = fabcolors) oncoplot( maf = laml, clinicalFeatures = 'FAB_classification', sortByAnnotation = TRUE, annotationColor = fabcolors )
3）注释多个临床数据
oncoplot( maf = laml, clinicalFeatures = c('FAB_classification','SEX','Overall_Survival_Status'), annotationColor = fabcolors)
以上就是如何使用R-maftools包完成瀑布图绘制的简单介绍，然后基本上也就完成了文献的Fig1 🤭。
开篇就说了maftools对于组学数据的分析和展示来说，算是一个宝藏“R包”，因此后续还会介绍其他的汇总，分析和可视化功能。
