本文首發於公眾號“生信補給站”，https://mp.weixin.qq.com/s/WG4JHs9RSm5IEJiiGEzDkg
之前介紹了使用maftools | 從頭開始繪制發表級oncoplot（瀑布圖） R-maftools包繪制組學突變結果（MAF）的oncoplot或者叫“瀑布圖”，以及一些細節的更改和注釋。
本文繼續介紹maftools對於MAF文件的其他應用，為更易理解和重現，本次使用TCGA下載的LIHC數據。
一 數據部分
setwd("C:\\Users\\Maojie\\Desktop\\maftools-V2\\")
library(maftools)
laml.maf = read.csv("TCGA.LIHC.mutect.maf.csv",header=TRUE)
#本次只展示maf的一些統計繪圖，只讀入組學數據，不添加臨床數據
laml = read.maf(maf = laml.maf)
#查看數據的基本情況
laml
An object of class MAF
ID summary Mean Median
1: NCBI_Build 1 NA NA
2: Center 1 NA NA
3: Samples 364 NA NA
4: nGenes 12704 NA NA
5: Frame_Shift_Del 1413 3.893 3
6: Frame_Shift_Ins 551 1.518 1
7: In_Frame_Del 277 0.763 0
8: In_Frame_Ins 112 0.309 0
9: Missense_Mutation 28304 77.972 63
10: Nonsense_Mutation 1883 5.187 4
11: Nonstop_Mutation 45 0.124 0
12: Splice_Site 1051 2.895 2
13: Translation_Start_Site 65 0.179 0
14: total 33701 92.840 75
可以將MAF文件的gene ，sample的 summary 的信息，輸出到laml前綴的summary文件
write.mafSummary(maf = laml, basename = 'laml')
laml_geneSummary.txt
laml_sampleSummary.txt
二 繪圖部分
1,首先繪制MAF文件的整體結果圖
plotmafSummary(maf = laml, rmOutlier = TRUE, addStat = 'median', dashboard = TRUE, titvRaw = FALSE)
2,oncoplot圖
#oncoplot for top ten mutated genes.
oncoplot(maf = laml, top = 20)
添加SCNA信息，添加P值信息，添加臨床注釋信息，更改顏色等可參考 鏈接 。。
3 Oncostrip
可以使用 oncostrip 函數展示特定基因在樣本中的突變情況，此處查看肝癌中關注較多的'TP53','CTNNB1', 'ARID1A'三個基因，如下：
oncostrip(maf = laml, genes = c('TP53','CTNNB1', 'ARID1A'))
4 Transition ， Transversions
titv函數將SNP分類為Transitions_vs_Transversions，並以各種方式返回匯總表的列表。匯總數據也可以顯示為一個箱線圖，顯示六種不同轉換的總體分布，並作為堆積條形圖顯示每個樣本中的轉換比例。
laml.titv = titv(maf = laml, plot = FALSE, useSyn = TRUE)
#plot titv summary
plotTiTv(res = laml.titv)
5 Rainfall plots
使用rainfallPlot參數繪制rainfall plots，展示超突變的基因組區域。detectChangePoints設置為TRUE，rainfall plots可以突出顯示潛在變化的區域.
brca <- system.file("extdata", "brca.maf.gz", package = "maftools")
brca = read.maf(maf = brca, verbose = FALSE)
rainfallPlot(maf = laml, detectChangePoints = TRUE, pointSize = 0.6)
6 Compare mutation load against TCGA cohorts
通過tcgaComapre函數實現laml（自有群體）與TCGA中已有的33個癌種隊列的突變負載情況的比較。
#cohortName 給輸入的隊列命名
laml.mutload = tcgaCompare(maf = laml, cohortName = 'LIHC-2')
7 Genecloud
使用 geneCloud參數繪制基因雲，每個基因的大小與它突變的樣本總數成正比。
geneCloud(input = laml, minMut = 15)
8 Somatic Interactions
癌症中的許多引起疾病的基因共同發生或在其突變模式中顯示出強烈的排他性。可以使用somaticInteractions函數使用配對Fisher 's精確檢驗來分析突變基因之間的的co-occurring 或者exclusiveness。
#exclusive/co-occurance event analysis on top 10 mutated genes.
Interact <- somaticInteractions(maf = laml, top = 25, pvalue = c(0.05, 0.1))
#提取P值結果
Interact$gene_sets
gene_set pvalue
1: CTNNB1, AXIN1, TP53 0.0001486912
2: CTNNB1, TP53, ARID1A 0.0018338597
3: AXIN1, TP53, APOB 0.0087076043
4: CSMD3, AXIN1, ALB 0.0130219628
5: AXIN1, TP53, ALB 0.0173199619
6: CTNNB1, AXIN1, ARID1A 0.0363739468
9 Comparing two cohorts (MAFs)
由於癌症的突變模式各不相同，因此可是 mafComapre參數比較兩個不同隊列的差異突變基因，檢驗方式為fisher檢驗。
#輸入另一個 MAF 文件
Our_maf <- read.csv("Our_maf.csv",header=TRUE)
our_maf = read.maf(maf = Our_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 = laml, m2 = our_maf, m1Name = 'LIHC', m2Name = 'OUR', minMut = 5)
print(pt.vs.rt)
-
result部分會有每個基因分別在兩個隊列中的個數以及P值和置信區間等信息。
-
SampleSummary 會有兩個隊列的樣本數。
1) Forest plots
比較結果繪制森林圖
forestPlot(mafCompareRes = pt.vs.rt, pVal = 0.01, color = c('royalblue', 'maroon'), geneFontSize = 0.8)
10 Oncogenic 信號通路
`OncogenicPathways 功能查看顯著富集通路
OncogenicPathways(maf = laml)
#會輸出統計結果
Pathway alteration fractions
Pathway N n_affected_genes fraction_affected
1: RTK-RAS 85 68 0.8000000
2: WNT 68 55 0.8088235
3: NOTCH 71 52 0.7323944
4: Hippo 38 30 0.7894737
5: PI3K 29 24 0.8275862
可以對上面富集的通路中選擇感興趣的進行完成的突變展示：
PlotOncogenicPathways(maf = laml, pathways = "PI3K")
好了，以上就是使用maftools包對MAF格式的組學數據的匯總，分析，可視化。
后台回復“maf文件”即可獲得示例的maf文件和代碼
【覺得不錯，右下角點擊賞個“在看”，轉發就是贊賞，謝謝！】