2023年11月27日发(作者：)

TCGA数据库MAF⽂件的分析以及可视化（maftools包的使⽤）

这两天“王院长”帮别⼈分析snv数据，顺便给我推荐了⼀个R包：maftools。我看了⼀下官⽹，感觉是个⾮常强⼤的R

包。那必须要学习起来~实际上⽹上已经有很多教程了，但都是⽤maftools⾃带的⽰例数据来演⽰的。这⾥我⾃⼰去

TCGA⽹站上随便挑了⼀个maf⽂件来做练习（⽤⽰例数据来练习就不知道哪⾥会报错、哪⼀步需要注意了，事实证明在

练习过程中，有许多步骤都需要根据你的MAF⾥包含的信息进⾏参数修改，否则是要报错的。如果只是使⽤⽰例数据，

这些需要注意的地⽅都不会被发现的）。maftools的官⽅⽹站：here

从TCGA⽹站上随便选了⼀个SNV的MAF⽂件，并且把对应的临床数据下载下来（这第⼀步就是需要注意的地⽅，但是如果你⽤的是maftools⾃

带的数据跑流程，你不会知道第⼀步就是⼀个坑）:

在excel⾥把临床数据的列名改⼀下，否则后续分析会报错

读取MAF⽂件

这⼀步构建maf对象同样⾮常重要，也是⼀个需要注意的点。

#这⾥需要注意的是：下⾯⾥必须加上isTCGA = TRUE

#否则这⼀步会报错，有兴趣的同学可以去掉这⼀个参数设置运⾏⼀下，看看是什么报错。

#另外，如果你看其他教程⾥使⽤的是maftools⾃带的数据，这⼀步的参数也是没有isTCGA=TRUE的，这就是为什么我没有使⽤⾃带数据来练习，因为你不知道在分析你⾃⼰数据

> maf <- ("",clinicalData = "",isTCGA = TRUE)

> maf

An object of class MAF

> getGeneSummary(maf)

Hugo_Symbol Frame_Shift_Del Frame_Shift_Ins In_Frame_Del In_Frame_Ins Missense_Mutation Nonsense_Mutation Nonstop_Mutation

1: TP53 49 21 8 3 220 68 0

2: TTN 6 5 0 1 310 18 0

3: FAT1 25 23 3 0 26 61 0

4: CDKN2A 15 3 3 0 20 51 0

5: MUC16 2 2 0 3 109 11 0

---

15432: ZNHIT1 0 0 0 0 1 0 0

15433: ZNHIT2 0 0 0 0 1 0 0

15434: ZP1 0 0 0 0 1 0 0

15435: ZRSR1 0 0 0 0 1 0 0

15436: ZSCAN30 0 0 0 0 1 0 0

Splice_Site Translation_Start_Site total MutatedSamples AlteredSamples

在总结图⾥，列出了在这个MAF⽂件⾥所有的variant分类，对于这个样品来说，错义突变的⽐例是最⾼的；variant的类

型是SNP最多；SNP的突变⽅式是C到T的⽐例是最⾼的；平均每个样品⾥的variants数量是92个；以及top10的突变基

因。）

（2）Oncoplot图

oncoplot图也叫做瀑布图：

> oncoplot(maf = maf, top = 10)

上⾯图⾥⿊⾊的Multi-hit的意思是该基因在同⼀个样品⾥有多于⼀个的突变。

当你的样品数很多的时候，⽐如像我这个MAF⽂件有508个样品，你还可以添加⼀个参数，让样品之间没有灰⾊的框框，看起来会更好看⼀些

（参考： 肿瘤变异数据分析和可视化⼯具maftools：突变数据下载和可视化）

> oncoplot(maf = maf, top = 10,borderCol=NULL)#样品不显⽰边界

如果你不喜欢图⾥的颜⾊，也可以根据这篇⽂章⾥的代码来更改颜⾊：maftools | 从头开始绘制发表级oncoplot（瀑布

图）

对于oncoplot这个函数，你还可以在图上添加其他的信息，可以参考官⽹：here

（3）转换和颠换的可视化

titv

功能可以把SNP分类为转换和颠换，并且返回⼀个总结性的table。这个总结性的数据也可以被可视化，并且可以展⽰每⼀个样品⾥每⼀种类

型的⽐例：

> = titv(maf = maf, plot = FALSE, useSyn = TRUE)

#plot titv summary

> plotTiTv(res = )

（3）氨基酸改变的Lollipop plots图（棒棒糖图）

如果你想画氨基酸change的棒棒糖图，⾸先要确保你的MAF⽂件⾥有相应的信息。然⽽MAF⽂件没有明确的标准对于氨基酸changes这⼀项信

息的命名，在不同的study⾥可能有着不同的命名⽅法。默认情况下，棒棒糖图是根据MAF⽂件⾥AAChanges这⼀列画的。如果你的MAF⽂件

⾥没有AAChanges这⼀项，那么会返回⼀个warning message提⽰你你的MAF⽂件⾥都有什么信息。⽐如我使⽤的这个MAF⽂件⾥，如果根

#在官⽹的例⼦⾥，他们⽤的MAF⽂件⾥关于氨基酸change的命名是“Protein_change”，⽽在我的数据⾥不是这个命名，所以出现以下报错：

> lollipopPlot(maf = maf, gene = 'TTN', AACol = 'Protein_Change', showMutationRate = TRUE)

Available fields:

[1] "Hugo_Symbol" "Entrez_Gene_Id"

[3] "Center" "NCBI_Build"

[5] "Chromosome" "Start_Position"

[7] "End_Position" "Strand"

[9] "Variant_Classification" "Variant_Type"

[11] "Reference_Allele" "Tumor_Seq_Allele1"

> lollipopPlot(maf = maf, gene = 'FAT1', AACol = 'HGVSp_Short', showMutationRate = TRUE)

这⾥要说明的是：Lollipop plots图默认使⽤的是基因的最长的那个转录本

也可以针对某⼀个突变位点进⾏标注：

> lollipopPlot(maf = maf, gene = 'FAT1', showDomainLabel = FALSE, labelPos = 4136)

（4）Rainfall plots图

肿瘤基因组，特别是实体肿瘤，具有局部超突变（hyper-mutations）。这种超突变基因组区域可以通过在线性基因组尺度上绘制变异间距离来

进⾏可视化。这种图通常称为Rainfall plots，我们可以使⽤来画图。如果detectChangePoints参数设置为TRUE，rainfallPlot图还

rainfallPlot

可以突出显⽰突变间距离潜在变化所在的区域。NOTE：因为我使⽤的数据⾥没有检测到changepoints，所以我就借⽤⼀下官⽹的图

（5）与TCGA cohort⽐较突变负荷

将⾃⼰的MAF与⾃带的33个TCGA cohorts进⾏⽐较：

> d = tcgaCompare(maf = maf, cohortName = 'Example-hnsc')

Performing pairwise t-test for differences in mutation burden..

Warning message:

In FUN(X[[i]], ...) : Removed 0 samples with zero mutations.

（6）变异等位基因频率（VAF）可视化

NOTE:由于我的数据⾥没有这⼀项信息，所以这⾥的代码会有所不同。如果你使⽤的MAF⽂件⾥也没有相关信息也没有关系，将

i_TumorVAF_WU

参数设置成vafCol = NULL，会⾃动计算并进⾏可视化：

plotVaf

> plotVaf(maf = maf, vafCol = NULL)

（7）Somatic Interactions

对top基因的突变事件画相关性图：

> somaticInteractions(maf = maf, top = 25, pvalue = c(0.05, 0.01))

gene1 gene2 pValue oddsRatio 00 11 01 10

1: CDKN2A TP53 3.652900e-09 5.7583311 158 91 249 10

2: DMD PAPPA2 2.851697e-07 7.9922491 436 15 29 28

3: TTN RYR2 3.054658e-07 4.8395698 299 38 15 156

4: FLG MUC16 2.233155e-06 3.9689241 373 28 69 38

5: PCLO CSMD3 2.291170e-06 3.8123129 372 29 60 47

---

296: PAPPA2 PIK3CA 1.000000e+00 0.9221460 385 7 79 37

maftools有⼀个功能，它可以鉴定⼀个MAF⽂件⾥的肿瘤driver基因。这个功能根据oncodriveCLUST算法进⾏运算。这个概念是基

oncodrive

于⼤多数致癌基因的变异在少数特定位点(即热点)富集。

> = oncodrive(maf = maf, AACol = 'HGVSp_Short', minMut = 5, pvalMethod = 'zscore')

plotOncodrive将结果绘制成散点图，点的⼤⼩与基因中发现的簇数成⽐例。x轴显⽰了在这些clusters中观察到的突变数量(或突变的部分)。

（9）添加pfam domain区域

maftools附带pfamDomains功能，可以在氨基酸changes⾥添加pfam域信息。根据受影响的结构域，还总结了氨基酸结构域的变化。这有助

于了解特定癌症群体中哪些domain最容易受到影响。

> = pfamDomains(maf = maf, AACol = 'HGVSp_Short', top = 10)

#Protein summary (Printing first 7 columns for display convenience)

> $proteinSummary[,1:7, with = FALSE]

HGNC AAPos Variant_Classification N total fraction

1: PIK3CA 545 Missense_Mutation 24 88 0.27272727

2: CDKN2A 80 Nonsense_Mutation 20 104 0.19230769

3: PIK3CA 542 Missense_Mutation 18 88 0.20454545

4: PIK3CA 1047 Missense_Mutation 15 88 0.17045455

5: TP53 273 Missense_Mutation 13 398 0.03266332

---

> adcc_clinical_ = clinicalEnrichment(maf = maf, clinicalFeature = 'ajcc_clinical_stage')

Sample size per factor in ajcc_clinical_stage:

'-- Stage I Stage II Stage III Stage IVA Stage IVB Stage IVC

14 20 95 104 260 9 6

> adcc_clinical_$groupwise_comparision[p_value < 0.001]

Hugo_Symbol Group1 Group2 n_mutated_group1 n_mutated_group2

1: ANO9 '-- Rest 3 of 14 4 of 494

2: ADAMTSL1 Stage III Rest 8 of 104 4 of 404

3: INPP4B Stage IVB Rest 3 of 9 8 of 499

上⾯的图显⽰了在这个MAF⽂件中涉及到的前5个基因，以及潜在的可给药的基因类别。也可以提取药物-基因相互作⽤的信息（基因名称、相互

作⽤类型、药物名称）。例如，下⾯是已知/报道的药物与RYR2相互作⽤的结果。NOTE：这个药物-基因信息只⽤于研究使⽤，并不具有临床指

导意义

> = drugInteractions(genes = "RYR2", drugs = TRUE)

Number of claimed drugs for given genes:

Gene N

1: RYR2 9

> [,.(Gene, interaction_types, drug_name, drug_claim_name)]

Gene interaction_types drug_name

> OncogenicPathways(maf = maf)

Pathway N n_affected_genes fraction_affected

肿瘤抑制基因标记为红⾊，oncogene标记为蓝⾊。横轴上的点是样品数量，但不是全部的508个样品，只是PI3K信号通

路有突变的样品。

（13）肿瘤异质性和MATH分数

肿瘤通常是异质性的，即由多个克隆组成。这种异质性可以通过聚类变异等位基因频率来推断。函数使⽤vaf信息聚类变异(使⽤

inferHeterogeneity

mclust)，以推断克隆。默认情况下，函数查找MAF⽂件⾥包含的vaf信息的（t_vaf）。但是，如果你的MAF⽂件⾥的vaf信息

inferHeterogeneity

名称不是t_vaf，可以使⽤参数vafCol⼿动指定。在我选的这个MAF⽂件⾥并没有vaf信息，参数设置vafCol = NULL，它会⾃动计算。

#从MAF⽂件⾥随便选取⼀个样品来查看肿瘤异质性

> = inferHeterogeneity(maf = maf, tsb = 'TCGA-CN-4741', vafCol = NULL)

t_vaf field is missing, but found t_ref_count & t_alt_count columns. Estimating vaf..

Processing TCGA-CN-4741..

> print($clusterMeans)

Tumor_Sample_Barcode cluster meanVaf

1: TCGA-CN-4741 2 0.3846020

2: TCGA-CN-4741 1 0.0878295

⼀种简单的肿瘤内异质性的定量测量，它计算出vaf分布的宽度。MATH的分数越⾼，代表着越差的outcome（可以理解为预后效果）。所以

MATH分数也可以作为⽣存分析的代理变量。

（14）Mutational Signatures

每⼀种癌症在其发展过程中都会留下以特定的核苷酸替换为特征的标志。Alexandrov等⼈从7000多个癌症样本中发现了这种突变特征

（Mutational Signatures）。这些特征可以通过分解核苷酸置换矩阵来提取，并根据突变碱基周围的直接碱基将其划分为96个置换类别。提取

的特征也可以与那些已经经过验证的特征进⾏⽐较。

特征分析的第⼀步是获得突变碱基周围的碱基，形成突变矩阵。注意:早期版本的maftools需要⼀个fasta⽂件作为输⼊。从1.8.0版本开

始，BSgenome对象被⽤于更快的序列提取。

#这⾥需要注意的是，你要知道你的样品是根据哪⼀个基因组⽐对来的，使⽤的包也会不⼀样，

#这⾥官⽹⽤的是hg19，然⽽我这⾥如果⽤hg19就会报错

#当你⽤maf查看你的MAF⽂件信息summary的时候，如果NCBI_build那⼀项是hg19，这⾥就要⽤hg19;我的MAF⽂件是⽤hg38来build的，所以这⾥也要⽤相应的基因组

> BiocManager::install("38")

> library(38, quietly = TRUE)

> = trinucleotideMatrix(maf = maf, ref_genome = "38")

-Extracting 5' and 3' adjacent bases

-Extracting +/- 20bp around mutated bases for background C>T estimation

-Estimating APOBEC enrichment scores

--Performing one-way Fisher's test for APOBEC enrichment

---APOBEC related mutations are enriched in 43.984 % of samples (APOBEC enrichment score > 2 ; 223 of 507 samples)

-Creating mutation matrix

--matrix of dimension 507x96

上⾯的分析过程主要是分两步：

（1）计算APOBEC富集分数

（2）准备特征分析的突变矩阵

APOBEC诱导的突变在实体肿瘤中更为常见，主要与TCW motif中发⽣的C>T转换事件有关。使⽤Roberts等⼈所描述的⽅法计算上述命令中的

APOBEC富集分数。简单地说，在⼀个给定的样本中，将TCW motif中发⽣的C>T突变与所有C>T突变的富集情况的⽐例与背景胞嘧啶和发⽣在

突变碱基20bp内的TCWs进⾏⽐较。

我们还可以分析APOBEC富集和⾮APOBEC富集的样品在突变模式上的差异。函数采⽤trinucleotideMatrix计算APOBEC富集分

plotApobecDiff

数，将样本分为APOBEC富集和⾮APOBEC富集。分组后，⽐较这两组，以确定改变的基因的差异。

> plotApobecDiff(tnm = , maf = maf, pVal = 0.01)

Showing top 10 of 24 differentially mutated genes

> library('NMF')

> = estimateSignatures(mat = , nTry = 10) #这⾥我⽤的nTry是10，官⽹⽤的是6

-Running NMF for 10 ranks

Compute NMF rank= 2 ... + measures ... OK

Compute NMF rank= 3 ... + measures ... OK

Compute NMF rank= 4 ... + measures ... OK

Compute NMF rank= 5 ... + measures ... OK

Compute NMF rank= 6 ... + measures ... OK

> = extractSignatures(mat = , n = 6)

特征富集分析

上⾯鉴定出的特征可以进⼀步分配到样品上，并且使⽤signatureEnrichment功能进⾏富集分析，鉴定每⼀个特征的突变富集。

> <- signatureEnrichment(maf = maf, sig_res = )

Running k-means for signature assignment..

Performing pairwise and groupwise comparisions..

>plotEnrichmentResults(enrich_res = , pVal = 0.005)