如果你做了全基因組、全外顯子組或者targeted sequencing,，拿到了一堆基因變異位點(diǎn)（SNP）,，下一步你該分析點(diǎn)什么？

如果還沒(méi)有頭緒,，你可以認(rèn)識(shí)一下這個(gè)神器：maftools,，不僅馬上有了分析思路，還能收獲一堆結(jié)果～

這是一個(gè)R包,，先來(lái)看一下這個(gè)包都能干點(diǎn)什么：

具體用法可以看官網(wǎng)說(shuō)明書：
http:///packages/release/bioc/manuals/maftools/man/maftools.pdf

這個(gè)包好在哪里呢,？

首先，MAF是非常常見(jiàn)的描述基因突變的文件形式,，只要拿到這么一個(gè)文件,，就可以做一系列的突變分析。

第二,，最基本的和稍微進(jìn)階一點(diǎn)的突變分析套路,，maftools已經(jīng)給你安排的明明白白的了，這些結(jié)果也夠研究一陣的了,。

那么什么是MAF(Mutation Annotation Format ),，可以看下這個(gè)鏈接https://docs.gdc./Data/File_Formats/MAF_Format/，里面除了說(shuō)了什么MAF外,，還對(duì)MAF格式的每個(gè)列名做了解釋說(shuō)明,，一般來(lái)說(shuō)我們從TCGA下載是突變文件都是MAF格式的。

如果是VCF格式的突變分析結(jié)果,，那么可以用vcf2maf工具從vcf轉(zhuǎn)化為MAF格式,，至于為什么我想說(shuō)要轉(zhuǎn)為MAF格式呢，因?yàn)橛袀€(gè)maftools這個(gè)R包,，其基礎(chǔ)MAF格式可以做很多分析以及可視化工作,，特別方便！

在R中安裝maftools包

有了包,，下一步就是讀入數(shù)據(jù)開(kāi)始分析啦～

先看看你需要給它什么樣的數(shù)據(jù)？

必須給的文件包括（Required input files）

• MAF文件 - 可以將MAF壓縮成 .gz結(jié)尾的文件作為input,，也可以不壓縮,。

• 如果你還需要臨床分析,，再給一個(gè)臨床數(shù)據(jù)文件（可選）。

• 如果有拷貝數(shù)數(shù)據(jù)也可以給,，maftools接受GISTIC的輸出文件（這一類文件也可選）,。

讀入MAF文件

library(maftools)

這里我用一個(gè)示例數(shù)據(jù)，這個(gè)示例數(shù)據(jù)只要值安裝了maftools就可以用,，讀入文件用 read.maf這個(gè)函數(shù)

讀入MAF后，可以總結(jié)一下這個(gè)文件中樣本（samples）的情況：

# Shows sample summary.getSampleSummary(laml)# 結(jié)果為每個(gè)樣本的變異類型總結(jié)：

總結(jié)一下這個(gè)文件中基因（Gene）的情況：

# Shows gene summary.getGeneSummary(laml)# 結(jié)果為每個(gè)基因的變異類型總結(jié)：

對(duì)MAF還不是很了解的同學(xué),，可以看看這里,，這個(gè)命令展示了MAF文件中所有列的名稱?？梢粤私庖幌翸AF都包含哪些內(nèi)容,。

# 展示一下MAF所有列都是什么# Shows all fields in MAFgetFields(laml)

如果有臨床信息，可以看一下臨床信息什么樣子：

# shows clinical data associated with samplesgetClinicalData(laml)

4.1 MAF概覽

plotmafSummary函數(shù)將每個(gè)樣本中的變體數(shù)顯示為堆積條形圖,，將變體類型顯示為 Variant_Classification 匯總的箱形圖,。我們可以在堆積的條形圖中添加平均線或中線，以顯示整個(gè)群組中變體的平均值/中值數(shù),。

4.2 突變?nèi)皥D（Oncoplots）

#這里只展示了突變檢出率top10的基因,。參數(shù)可以自己調(diào)整#oncoplot for top ten mutated genes.oncoplot(maf = laml, top = 10)

如果只想挑幾個(gè)感興趣的基因看突變譜也可以：

4.3 Transition and Transversions.

titv函數(shù)將SNP分類為轉(zhuǎn)換和顛換,，并以各種方式返回匯總表的列表。匯總數(shù)據(jù)也可以顯示為一個(gè)箱線圖,，顯示六種不同轉(zhuǎn)換的總體分布，并作為堆積條形圖顯示每個(gè)樣本中的轉(zhuǎn)換比例,。

laml.titv = titv(maf = laml, plot = FALSE, useSyn = TRUE)#plot titv summaryplotTiTv(res = laml.titv)

4.4 氨基酸改變的棒棒糖圖??

棒棒糖圖顯示了蛋白質(zhì)結(jié)構(gòu)上的突變點(diǎn),。許多癌基因，某些位點(diǎn)比其他位點(diǎn)更容易突變,，變異頻率更高,，這些位置被認(rèn)為是突變熱點(diǎn)(hot-spots),。棒棒糖圖可以用于顯示它們以及其他突變。我們可以使用函數(shù)lollipopPlot繪制這樣的圖,。這個(gè)功能要求我們?cè)趍af文件中有氨基酸改變信息,。然而，MAF文件沒(méi)有關(guān)于命名氨基酸變化字段的明確指南,，不同的研究具有不同的氨基酸變化的字段（或列）名稱,。默認(rèn)情況下，lollipopPlot查找列AAChange,，如果在MAF文件中找不到它,，則會(huì)打印所有可用字段并顯示警告消息。對(duì)于以下示例,，MAF文件包含字段/列名稱“Protein_Change”下的氨基酸變化,。我們將使用參數(shù)AACol手動(dòng)指定它。此函數(shù)還將繪圖作為ggplot對(duì)象返回,，如果需要,，用戶稍后可以修改該對(duì)象。

4.5 瀑布圖（Rainfall plots）

癌癥基因組,，特別是實(shí)體瘤會(huì)有些位點(diǎn)存在超突變（hyper-mutations）,。rainfallPlot 這個(gè)函數(shù)可以將這一現(xiàn)象通過(guò)瀑布圖體現(xiàn)出來(lái)。

brca <- system.file('extdata', 'brca.maf.gz', package = 'maftools')brca = read.maf(maf = brca, verbose = FALSE)

4.6 Plotting VAF

如果你的MAF文件有VAF信息的話,，可以簡(jiǎn)單畫個(gè)箱線圖看看所有樣本在不同gene上的VAF分布,，vafCol就是指定MAF列中的VAF的列名。

plotVaf(maf = laml, vafCol = 'i_TumorVAF_WU')

5.1 基因之間的伴發(fā)互斥關(guān)系

深綠色的為樣本中傾向共同發(fā)生突變的基因；而深黃色為互斥基因,，也就是在一個(gè)樣本中,，“有我沒(méi)他,，有他沒(méi)我”的基因。

##      gene1  gene2       pValue oddsRatio  00 11 01 10              Event
##   1: ASXL1  RUNX1 0.0001541586 55.215541 176  4 12  1       Co_Occurence
##   2:  IDH2  RUNX1 0.0002809928  9.590877 164  7  9 13       Co_Occurence
##   3:  IDH2  ASXL1 0.0004030636 41.077327 172  4  1 16       Co_Occurence
##   4:  FLT3   NPM1 0.0009929836  3.763161 125 17 16 35       Co_Occurence
##   5:  SMC3 DNMT3A 0.0010451985 20.177713 144  6 42  1       Co_Occurence
##              pair event_ratio
##   1: ASXL1, RUNX1        4/13
##   2:  IDH2, RUNX1        7/22
##   3:  ASXL1, IDH2        4/17
##   4:   FLT3, NPM1       17/51
##   5: DNMT3A, SMC3        6/43
##  ---
## 296: ASXL1, PLCE1         0/9
## 297: FAM5C, RAD21        0/10
## 298: FAM5C, PLCE1         0/9
## 299: PLCE1, RAD21         0/9
## 300:  EZH2, PLCE1         0/

5.2 生存分析 （Survival analysis）

1) Mutation in any given genes

2）Predict genesets associated with survival

鑒別引起更差生存的基因集

#Using top 20 mutated genes to identify a set of genes (of size 2) to predict poor prognostic groupsprog_geneset = survGroup(maf = laml, top = 20, geneSetSize = 2, time = 'days_to_last_followup', Status = 'Overall_Survival_Status', verbose = FALSE)

##     Gene_combination P_value    hr  WT Mutant
##  1:      FLT3_DNMT3A 0.00104 2.510 164     18
##  2:      DNMT3A_SMC3 0.04880 2.220 176      6
##  3:      DNMT3A_NPM1 0.07190 1.720 166     16
##  4:      DNMT3A_TET2 0.19600 1.780 176      6
##  5:        FLT3_TET2 0.20700 1.860 177      5
##  6:        NPM1_IDH1 0.21900 0.495 176      6
##  7:      DNMT3A_IDH1 0.29300 1.500 173      9
##  8:       IDH2_RUNX1 0.31800 1.580 176      6
##  9:        FLT3_NPM1 0.53600 1.210 165     17
## 10:      DNMT3A_IDH2 0.68000 0.747 178      4
## 11:      DNMT3A_NRAS 0.99200 0.986 178      4

Above results show a combination (N = 2) of genes which are associated with poor survival (P < 0.05). We can draw KM curve for above results with the function mafSurvGroup

5.3 兩個(gè)隊(duì)列的比較

不同癌癥的突變譜是不同的,，我們可以通過(guò)比較兩個(gè)隊(duì)列來(lái)鑒別不同的突變基因,。比如，Madan et. al 文章中顯示,，復(fù)發(fā)的急性早幼粒細(xì)胞白血?。╮elapsed APL）有PML和RARA突變的趨勢(shì)，而原發(fā)急性早幼粒細(xì)胞白血?。╬rimary APL）則沒(méi)有,。

這種兩個(gè)隊(duì)列間的差異可以用mafCompare函數(shù)來(lái)鑒別，函數(shù)在兩個(gè)隊(duì)列多有基因中用fisher test鑒別差異化突變基因,。

#讀入Primary APL MAFprimary.apl = system.file('extdata', 'APL_primary.maf.gz', package = 'maftools')primary.apl = read.maf(maf = primary.apl)#讀入Relapse APL MAFrelapse.apl = system.file('extdata', 'APL_relapse.maf.gz', package = 'maftools')relapse.apl = read.maf(maf = relapse.apl）#兩個(gè)MAF比較#Considering only genes which are mutated in at-least in 5 samples in one of the cohort to avoid bias due to genes mutated in single sample.pt.vs.rt <- mafCompare(m1 = primary.apl, m2 = relapse.apl, m1Name = 'Primary', m2Name = 'Relapse', minMut = 5)print(pt.vs.rt)

5.3.1 Forest plots

forestPlot(mafCompareRes = pt.vs.rt, pVal = 0.1, color = c('royalblue', 'maroon'), geneFontSize = 0.8)

5.3.2 Co-onco plots

另外一種可視化辦法，就是將基因在兩個(gè)不同隊(duì)列中的oncoplot并排畫出來(lái),?？梢酝瑫r(shí)比對(duì)同一基因在不同隊(duì)列中的突變頻率、突變類型的不同,。

不得不說(shuō),，這個(gè)包做的挺全面的，實(shí)在是太貼心了,。

5.3.3 Lollipop plot-2

還記得第一部分可視化里面介紹的棒棒糖圖嗎,？

這里有了兩個(gè)隊(duì)列，還可以來(lái)一個(gè)兩個(gè)狀態(tài)比較的棒棒糖圖??

lollipopPlot2(m1 = primary.apl, m2 = relapse.apl, gene = 'PML', AACol1 = 'amino_acid_change', AACol2 = 'amino_acid_change', m1_name = 'Primary', m2_name = 'Relapse')

上半部分展示了原發(fā)中PML蛋白變異位點(diǎn),，下半部分展示了復(fù)發(fā)中PML蛋白變異位點(diǎn),。