maftools | 從頭開始繪制發表級oncoplot（瀑布圖）

本文首發於微信公眾號 **“ 生信補給站 ”** ，期待您的關注！！！

原文鏈接：https://mp.weixin.qq.com/s/G-0PtaoO6bYuhx_D_Rlrlw

對於組學數據的分析和展示來說，maftools算是一個寶藏“R包”，可用於MAF格式的組學數據的匯總，分析和可視化展示。

​ 隨着癌症基因組學的進步, 突變注釋格式+Specification) (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對象。

#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包”，因此后續還會介紹其他的匯總，分析和可視化功能。

【關注“生信補給站”公眾號，更多的精彩內容】

​