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

这两天“王院长”帮别人分析snv数据，顺便给我推荐了一个R包：maftools。我看了一下官网，感觉是个非常强大的R包。那必须要学习起来~实际上网上已经有很多教程了，但都是用maftools自带的示例数据来演示的。这里我自己去TCGA网站上随便挑了一个maf文件来做练习（用示例数据来练习就不知道哪里会报错、哪一步需要注意了，事实证明在练习过程中，有许多步骤都需要根据你的MAF里包含的信息进行参数修改，否则是要报错的。如果只是使用示例数据，这些需要注意的地方都不会被发现的）。maftools的官方网站：here从TCGA网站上随便选了一个SNV的MAF文件，并且把对应的临床数据下载下来（这第一步就是需要注意的地方，但是如果你用的是maftools自带的数据跑流程，你不会知道第一步就是一个坑）:

在excel里把临床数据的列名改一下，否则后续分析会报错读取MAF文件这一步构建maf对象同样非常重要，也是一个需要注意的点。#这里需要注意的是：下面read.maf里必须加上isTCGA = TRUE#否则这一步会报错，有兴趣的同学可以去掉这一个参数设置运行一下，看看是什么报错。#另外，如果你看其他教程里使用的是maftools自带的数据，这一步的参数也是没有isTCGA=TRUE的，这就是为什么我没有使用自带数据来练习，因为你不知道在分析你自己数据的时候会出现什么报错。> maf <- read.maf("TCGA.HNSC.mutect.somatic.maf.gz",clinicalData = "clinical.tsv",isTCGA = TRUE)-Reading|--------------------------------------------------||==================================================|-Validating-Silent variants: 35007 -Summarizing--Possible FLAGS among top ten genes: TTN MUC16 SYNE1-Processing clinical data-Finished in 10.3s elapsed (9.030s cpu)现在读取文件后，可以看看maf这个对象的简单的总结：> mafAn object of class MAF ID summary Mean Median 1: NCBI_Build GRCh38 NA NA 2: Center BI NA NA 3: Samples 508 NA NA 4: nGenes 15436 NA NA 5: Frame_Shift_Del 2055 4.045 2 6: Frame_Shift_Ins 1382 2.720 1 7: In_Frame_Del 328 0.646 0 8: In_Frame_Ins 251 0.494 0 9: Missense_Mutation 56861 111.931 7810: Nonsense_Mutation 4796 9.441 611: Nonstop_Mutation 63 0.124 012: Splice_Site 1482 2.917 213: Translation_Start_Site 84 0.165 014: total 67302 132.484 92查看每个样品的总结情况：> getSampleSummary(maf) Tumor_Sample_Barcode Frame_Shift_Del Frame_Shift_Ins In_Frame_Del In_Frame_Ins Missense_Mutation Nonsense_Mutation 1: TCGA-F7-A624 248 57 7 1 2174 63 2: TCGA-CV-7568 4 1 1 1 1645 105 3: TCGA-D6-6516 5 1 2 0 1331 78 4: TCGA-CR-7402 28 9 1 1 696 59 5: TCGA-CN-4723 8 4 0 0 564 73 --- 504: TCGA-CR-7392 0 0 0 0 11 0505: TCGA-D6-6827 1 1 2 0 7 0506: TCGA-CR-7391 0 0 0 0 4 0507: TCGA-CN-6017 1 0 0 0 1 0508: TCGA-CR-7393 1 0 0 0 0 0 Nonstop_Mutation Splice_Site Translation_Start_Site total 1: 2 60 1 2613 2: 0 34 3 1794 3: 0 24 0 1441 4: 0 27 1 822 5: 1 10 1 661 --- 504: 0 0 0 11505: 0 0 0 11506: 0 0 0 4507: 0 0 0 2508: 0 0 0 1查看基因的情况：> 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 015433: ZNHIT2 0 0 0 0 1 0 015434: ZP1 0 0 0 0 1 0 015435: ZRSR1 0 0 0 0 1 0 015436: ZSCAN30 0 0 0 0 1 0 0 Splice_Site Translation_Start_Site total MutatedSamples AlteredSamples 1: 28 1 398 340 340 2: 6 0 346 194 194 3: 5 0 143 114 114 4: 12 0 104 101 101 5: 2 0 129 97 97 --- 15432: 0 0 1 1 115433: 0 0 1 1 115434: 0 0 1 1 115435: 0 0 1 1 115436: 0 0 1 1 1除此之外，你还可以使用下面的代码查看其它信息：#shows clinical data associated with samples查看临床信息> getClinicalData(maf)#Shows all fields in MAF查看maf文件里所有的可用信息> getFields(maf)可视化（1）画MAF文件的summary图> plotmafSummary(maf = maf, rmOutlier = TRUE, addStat = 'median', dashboard = TRUE, titvRaw = FALSE)在总结图里，列出了在这个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。这个总结性的数据也可以被可视化，并且可以展示每一个样品里每一种类型的比例：> hnsc.titv = titv(maf = maf, plot = FALSE, useSyn = TRUE)#plot titv summary> plotTiTv(res = hnsc.titv)（3）氨基酸改变的Lollipop plots图（棒棒糖图）如果你想画氨基酸change的棒棒糖图，首先要确保你的MAF文件里有相应的信息。然而MAF文件没有明确的标准对于氨基酸changes这一项信息的命名，在不同的study里可能有着不同的命名方法。默认情况下，棒棒糖图是根据MAF文件里AAChanges这一列画的。如果你的MAF文件里没有AAChanges这一项，那么会返回一个warning message提示你你的MAF文件里都有什么信息。比如我使用的这个MAF文件里，如果根据官网上写的参数设置，R会返回给我一个warning message:#在官网的例子里，他们用的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" [13] "Tumor_Seq_Allele2" "dbSNP_RS" [15] "dbSNP_Val_Status" "Tumor_Sample_Barcode" [17] "Matched_Norm_Sample_Barcode" "Match_Norm_Seq_Allele1" [19] "Match_Norm_Seq_Allele2" "Tumor_Validation_Allele1" [21] "Tumor_Validation_Allele2" "Match_Norm_Validation_Allele1" [23] "Match_Norm_Validation_Allele2" "Verification_Status" [25] "Validation_Status" "Mutation_Status" [27] "Sequencing_Phase" "Sequence_Source" [29] "Validation_Method" "Score" [31] "BAM_File" "Sequencer" [33] "Tumor_Sample_UUID" "Matched_Norm_Sample_UUID" [35] "HGVSc" "HGVSp" [37] "HGVSp_Short" "Transcript_ID" [39] "Exon_Number" "t_depth" [41] "t_ref_count" "t_alt_count" [43] "n_depth" "n_ref_count" [45] "n_alt_count" "all_effects" [47] "Allele" "Gene" [49] "Feature" "Feature_type" [51] "One_Consequence" "Consequence" [53] "cDNA_position" "CDS_position" [55] "Protein_position" "Amino_acids" [57] "Codons" "Existing_variation" [59] "ALLELE_NUM" "DISTANCE" [61] "TRANSCRIPT_STRAND" "SYMBOL" [63] "SYMBOL_SOURCE" "HGNC_ID" [65] "BIOTYPE" "CANONICAL" [67] "CCDS" "ENSP" [69] "SWISSPROT" "TREMBL" [71] "UNIPARC" "RefSeq" [73] "SIFT" "PolyPhen" [75] "EXON" "INTRON" [77] "DOMAINS" "GMAF" [79] "AFR_MAF" "AMR_MAF" [81] "ASN_MAF" "EAS_MAF" [83] "EUR_MAF" "SAS_MAF" [85] "AA_MAF" "EA_MAF" [87] "CLIN_SIG" "SOMATIC" [89] "PUBMED" "MOTIF_NAME" [91] "MOTIF_POS" "HIGH_INF_POS" [93] "MOTIF_SCORE_CHANGE" "IMPACT" [95] "PICK" "VARIANT_CLASS" [97] "TSL" "HGVS_OFFSET" [99] "PHENO" "MINIMISED" [101] "ExAC_AF" "ExAC_AF_Adj" [103] "ExAC_AF_AFR" "ExAC_AF_AMR" [105] "ExAC_AF_EAS" "ExAC_AF_FIN" [107] "ExAC_AF_NFE" "ExAC_AF_OTH" [109] "ExAC_AF_SAS" "GENE_PHENO" [111] "FILTER" "CONTEXT" [113] "src_vcf_id" "tumor_bam_uuid" [115] "normal_bam_uuid" "case_id" [117] "GDC_FILTER" "COSMIC" [119] "MC3_Overlap" "GDC_Validation_Status" Error in lollipopPlot(maf = maf, gene = "TTN", AACol = "Protein_Change", : Column Protein_Change not found.那么根据我的MAF文件的情况，则要选择“HGVSp_Short”这一项来进行画图：> 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，我们可以使用rainfallPlot来画图。如果detectChangePoints参数设置为TRUE，rainfallPlot图还可以突出显示突变间距离潜在变化所在的区域。NOTE：因为我使用的数据里没有检测到changepoints，所以我就借用一下官网的图> rainfallPlot(maf = maf, detectChangePoints = TRUE, pointSize = 0.6)（5）与TCGA cohort比较突变负荷将自己的MAF与自带的33个TCGA cohorts进行比较：> hnsc.mutload = 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:由于我的数据里没有i_TumorVAF_WU这一项信息，所以这里的代码会有所不同。如果你使用的MAF文件里也没有相关信息也没有关系，将参数设置成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 37297: FAM135B KMT2D 1.000000e+00 0.9455631 389 7 72 40298: CASP8 USH2A 1.000000e+00 0.9510660 401 6 55 46299: PKHD1L1 CASP8 1.000000e+00 0.9716819 411 5 47 45300: XIRP2 CASP8 1.000000e+00 0.9961390 412 5 47 44 Event pair event_ratio 1: Co_Occurence CDKN2A, TP53 91/259 2: Co_Occurence DMD, PAPPA2 15/57 3: Co_Occurence RYR2, TTN 38/171 4: Co_Occurence FLG, MUC16 28/107 5: Co_Occurence CSMD3, PCLO 29/107 --- 296: Mutually_Exclusive PAPPA2, PIK3CA 7/116297: Mutually_Exclusive FAM135B, KMT2D 7/112298: Mutually_Exclusive CASP8, USH2A 6/101299: Mutually_Exclusive CASP8, PKHD1L1 5/92300: Mutually_Exclusive CASP8, XIRP2 5/91（8）鉴定肿瘤driver 基因maftools有一个功能oncodrive，它可以鉴定一个MAF文件里的肿瘤driver基因。这个功能根据oncodriveCLUST算法进行运算。这个概念是基于大多数致癌基因的变异在少数特定位点(即热点)富集。> hnsc.sig = oncodrive(maf = maf, AACol = 'HGVSp_Short', minMut = 5, pvalMethod = 'zscore')Estimating background scores from synonymous variants.. |===========================================================================| 100%Not enough genes to build background. Using predefined values. (Mean = 0.279; SD = 0.13)Estimating cluster scores from non-syn variants.. |===========================================================================| 100%Comapring with background model and estimating p-values..Done !> head(hnsc.sig) Hugo_Symbol Frame_Shift_Del Frame_Shift_Ins In_Frame_Del In_Frame_Ins1: HRAS 0 1 0 02: MAPK1 0 0 0 03: RAC1 0 0 0 04: HIST1H4K 0 0 0 05: HIST1H3C 0 0 0 06: MOS 0 0 0 0 Missense_Mutation Nonsense_Mutation Nonstop_Mutation Splice_Site1: 28 1 0 02: 7 1 0 03: 12 0 0 04: 5 1 0 05: 6 0 0 06: 6 0 0 0 Translation_Start_Site total MutatedSamples AlteredSamples clusters1: 0 30 29 29 22: 0 8 8 8 13: 0 12 12 12 34: 0 6 6 6 25: 0 6 6 6 16: 0 6 6 6 2 muts_in_clusters clusterScores protLen zscore pval fdr1: 26 0.8666667 189 4.520513 3.084500e-06 0.00067190962: 7 0.87500003604.584615 2.274114e-06 0.00067190963: 11 0.8627961 211 4.490739 3.548818e-06 0.00067190964: 5 0.7845178 103 3.888598 5.041239e-05 0.00715855915: 5 0.7500000 136 3.623077 1.455596e-04 0.01653557056: 4 0.6666667 346 2.982051 1.431620e-03 0.1355267092 fract_muts_in_clusters1: 0.86666672: 0.87500003: 0.91666674: 0.83333335: 0.83333336: 0.6666667> plotOncodrive(res = hnsc.sig, fdrCutOff = 0.01, useFraction = TRUE,labelSize = 1)plotOncodrive将结果绘制成散点图，点的大小与基因中发现的簇数成比例。x轴显示了在这些clusters中观察到的突变数量(或突变的部分)。（9）添加pfam domain区域maftools附带pfamDomains功能，可以在氨基酸changes里添加pfam域信息。根据受影响的结构域，还总结了氨基酸结构域的变化。这有助于了解特定癌症群体中哪些domain最容易受到影响。> hnsc.pfam = pfamDomains(maf = maf, AACol = 'HGVSp_Short', top = 10)#Protein summary (Printing first 7 columns for display convenience)> hnsc.pfam$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 --- 65993: ZZZ3 90 Missense_Mutation 1 6 0.1666666765994: pk 669 Missense_Mutation 1 4 0.2500000065995: pk 448 Missense_Mutation 1 4 0.2500000065996: pk 466 Missense_Mutation 1 4 0.2500000065997: pk 353 Missense_Mutation 1 4 0.25000000 DomainLabel 1: PI3Ka_I 2: ANK 3: PI3Ka_I 4: PI3Kc_IA_alpha 5: P53 --- 65993: <NA>65994: <NA>65995: <NA>65996: <NA>65997: <NA>#用下面的代码查看计算出的pfam区域的信息#Domain summary (Printing first 3 columns for display convenience)> hnsc.pfam$domainSummary[,1:3, with = FALSE] DomainLabel nMuts nGenes 1: 7tm_1 1490 495 2: COG5048 969 343 3: Cadherin_repeat 814 106 4: FN3 544 115 5: I-set 380 86 --- 5725: zf-UBR 1 15726: zf-ZPR1 1 15727: zf-piccolo 1 15728: zf-rbx1 1 15729: zinc finger 1 1（10）临床富集分析clinicalEnrichment功能可以把任何临床特征与相关的样品进行富集分析。 比如我想把临床分期与突变基因之间的关系表示出来：> adcc_clinical_stage.ce = 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_stage.ce$groupwise_comparision[p_value < 0.001] Hugo_Symbol Group1 Group2 n_mutated_group1 n_mutated_group21: ANO9 '-- Rest 3 of 14 4 of 4942: ADAMTSL1 Stage III Rest 8 of 104 4 of 4043: INPP4B Stage IVB Rest 3 of 9 8 of 4994: FMN2 '-- Rest 5 of 14 26 of 4945: ITGB1 Stage III Rest 7 of 104 3 of 404 p_value OR_low OR_high fdr1: 0.0005491163 0 0.1740052 12: 0.0005801867 0 0.3843828 13: 0.0005938635 0 0.1801677 14: 0.0008275002 0 0.3323002 15: 0.0008798140 0 0.3826840 1比如上面的ADAMTSL1 这个基因，主要在stageIII这个分期有突变的富集，同样的，INPP4B基因的突变主要在Stage IVB这个分期中出现。因为在临床样品的信息里，有些病人的信息并不完整，所以有一个群体是--，这里就会出现有的基因在“--”群体里有突变的富集。> plotEnrichmentResults(enrich_res = adcc_clinical_stage.ce, pVal = 0.001)（11）药物-基因相互作用drugInteractions 功能是用来检查药物-基因相互作用，以及基因的成药性，这个分析是由Drug Gene Interaction database汇编而成的。> dgi = drugInteractions(maf = maf, fontSize = 0.75)上面的图显示了在这个MAF文件中涉及到的前5个基因，以及潜在的可给药的基因类别。也可以提取药物-基因相互作用的信息（基因名称、相互作用类型、药物名称）。例如，下面是已知/报道的药物与RYR2相互作用的结果。NOTE：这个药物-基因信息只用于研究使用，并不具有临床指导意义> RYR2.dgi = drugInteractions(genes = "RYR2", drugs = TRUE)Number of claimed drugs for given genes: Gene N1: RYR2 9> RYR2.dgi[,.(Gene, interaction_types, drug_name, drug_claim_name)] Gene interaction_types drug_name1: RYR2 activator,channel blocker RYANODINE2: RYR2 activator SURAMIN3: RYR2 antagonist MAGNESIUM4: RYR2 activator ADENOSINE TRIPHOSPHATE5: RYR2 activator CAFFEINE6: RYR2 activator,antagonist cA27: RYR2 antagonist 8: RYR2 channel blocker PROCAINE9: RYR2 channel blocker RUTHENIUM RED drug_claim_name1: 43032: 17283: 7084: 17135: 4076: 7077: DAP0004998: 42919: 2432（12）Oncogenic信号通路富集OncogenicPathways功能可以查看样品里关于已知Oncogenic信号通路的富集情况：> OncogenicPathways(maf = maf) Pathway N n_affected_genes fraction_affected 1: NRF2 3 3 1.0000000 2: TP53 6 6 1.0000000 3: TGF-Beta 7 7 1.0000000 4: Cell_Cycle 15 10 0.6666667 5: MYC 13 11 0.8461538 6: PI3K 29 26 0.8965517 7: Hippo 38 36 0.9473684 8: NOTCH 71 58 0.8169014 9: WNT 68 61 0.897058810: RTK-RAS 85 80 0.9411765 Mutated_samples Fraction_mutated_samples 1: 57 0.11220472 2: 349 0.68700787 3: 47 0.09251969 4: 127 0.25000000 5: 37 0.07283465 6: 163 0.32086614 7: 279 0.54921260 8: 274 0.53937008 9: 171 0.3366141710: 278 0.54724409也可以对某一个完整的信号通路进行可视化：> PlotOncogenicPathways(maf = maf, pathways = "PI3K")肿瘤抑制基因标记为红色，oncogene标记为蓝色。横轴上的点是样品数量，但不是全部的508个样品，只是PI3K信号通路有突变的样品。（13）肿瘤异质性和MATH分数肿瘤通常是异质性的，即由多个克隆组成。这种异质性可以通过聚类变异等位基因频率来推断。inferHeterogeneity函数使用vaf信息聚类变异(使用mclust)，以推断克隆。默认情况下，inferHeterogeneity函数查找MAF文件里包含的vaf信息的（t_vaf）。但是，如果你的MAF文件里的vaf信息名称不是t_vaf，可以使用参数vafCol手动指定。在我选的这个MAF文件里并没有vaf信息，参数设置vafCol = NULL，它会自动计算。#从MAF文件里随便选取一个样品来查看肿瘤异质性> tcga.cn.4741.het = 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(tcga.cn.4741.het$clusterMeans) Tumor_Sample_Barcode cluster meanVaf1: TCGA-CN-4741 2 0.38460202: TCGA-CN-4741 1 0.08782953: TCGA-CN-4741 3 0.77905454: TCGA-CN-4741 outlier 0.2222222> plotClusters(clusters = tcga.cn.4741.het)从上图中可以清楚的看到有三个克隆，平均变异等位基因频率约为77.9%(主要克隆)，38.46%和8.8%。尽管突变等位基因频率的聚类给我们一个关于异质性的概念，它也可能是衡量异质性的程度的一个数值。MATH分数(在上面的图中是54.585)是一种简单的肿瘤内异质性的定量测量，它计算出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("BSgenome.Hsapiens.UCSC.hg38")> library(BSgenome.Hsapiens.UCSC.hg38, quietly = TRUE)> hnsc.tnm = trinucleotideMatrix(maf = maf, ref_genome = "BSgenome.Hsapiens.UCSC.hg38")-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富集的样品在突变模式上的差异。plotApobecDiff函数采用trinucleotideMatrix计算APOBEC富集分数，将样本分为APOBEC富集和非APOBEC富集。分组后，比较这两组，以确定改变的基因的差异。> plotApobecDiff(tnm = hnsc.tnm, maf = maf, pVal = 0.01)Showing top 10 of 24 differentially mutated genes$results Hugo_Symbol Enriched nonEnriched pval or ci.up ci.low adjPval 1: NSD1 13 44 0.0005989521 0.3383091 0.6608655 0.16251777 1 2: CASP8 34 18 0.0011488374 2.6532642 5.1495987 1.40831724 1 3: SMTN 8 0 0.0013041485 Inf Inf 2.22107156 1 4: FSTL5 2 18 0.0019547843 0.1341226 0.5700077 0.01491342 1 5: TGFBR2 16 5 0.0029623612 4.3012453 15.2586418 1.47518180 1 --- 10105: ZSCAN4 3 3 1.0000000000 1.2766417 9.6232806 0.16933136 110106: ZSCAN5B 3 4 1.0000000000 0.9546251 5.7067946 0.13838894 110107: ZSWIM2 2 3 1.0000000000 0.8479356 7.4690287 0.07022655 110108: ZWILCH 1 2 1.0000000000 0.6356815 12.2842428 0.01072353 110109: ZZZ3 3 3 1.0000000000 1.2766417 9.6232806 0.16933136 1$SampleSummary Cohort SampleSize Mean Median1: Enriched 223 133.915 1082: nonEnriched 284 131.824 81特征分析特征分析包括以下几步：（1）estimateSignatures：在一系列的值上运行NMF并测量拟合优度。（2）plotCophenetic：绘制elblow图并帮助你决定特征的最佳数量。最佳可能的特征是共同相关性显著下降的值。（3）extractSignatures：使用非负矩阵分解将矩阵分解为n个特征。根据上述两个步骤选择n的值。如果你可以自己估计n的值，可以跳过以上两步。（4）compareSignatures：从上一步骤中提取的特征可以与COSMIC数据库中已知的特征进行比较，并计算余弦相似度以确定最佳匹配。（5）plotSignatures：可视化> library('NMF')> hnsc.sign = estimateSignatures(mat = hnsc.tnm, nTry = 10) #这里我用的nTry是10，官网用的是6-Running NMF for 10 ranksCompute NMF rank= 2 ... + measures ... OKCompute NMF rank= 3 ... + measures ... OKCompute NMF rank= 4 ... + measures ... OKCompute NMF rank= 5 ... + measures ... OKCompute NMF rank= 6 ... + measures ... OKCompute NMF rank= 7 ... + measures ... OKCompute NMF rank= 8 ... + measures ... OKCompute NMF rank= 9 ... + measures ... OKCompute NMF rank= 10 ... + measures ... OK-Finished in 00:06:53 elapsed (29.9s cpu)> plotCophenetic(res = hnsc.sign)最佳可能值是y轴上相关值显著下降之前的值。那么上图的话，有两个显著下降，一个是2，另一个是在6。这里官网的示例数据里只有一个显著下降的值，并且官网在这里写到：对于突变负荷较高的实体肿瘤，如果有足够的样本数量，会有更多的特征。我所用的HNSC肿瘤突变负荷比官网使用的LAML突变负荷要高，所以我选择n=6。确定好n的值，就可以正式分析了。> hnsc.sig = extractSignatures(mat = hnsc.tnm, n = 6)-Running NMF for factorization rank: 6-Finished in9.390s elapsed (7.940s cpu)> hnsc.og30.cosm = compareSignatures(nmfRes = hnsc.sig, sig_db = "legacy")-Comparing against COSMIC signatures--------------------------------------Found Signature_1 most similar to COSMIC_13 Aetiology: APOBEC Cytidine Deaminase (C>G) [cosine-similarity: 0.832]--Found Signature_2 most similar to COSMIC_1 Aetiology: spontaneous deamination of 5-methylcytosine [cosine-similarity: 0.937]--Found Signature_3 most similar to COSMIC_4 Aetiology: exposure to tobacco (smoking) mutagens [cosine-similarity: 0.946]--Found Signature_4 most similar to COSMIC_7 Aetiology: UV exposure [cosine-similarity: 0.938]--Found Signature_5 most similar to COSMIC_6 Aetiology: defective DNA mismatch repair [cosine-similarity: 0.798]--Found Signature_6 most similar to COSMIC_5 Aetiology: Unknown [cosine-similarity: 0.805]------------------------------------> hnsc.v3.cosm = compareSignatures(nmfRes = hnsc.sig, sig_db = "SBS")-Comparing against COSMIC signatures--------------------------------------Found Signature_1 most similar to SBS2 Aetiology: APOBEC Cytidine Deaminase (C>T) [cosine-similarity: 0.714]--Found Signature_2 most similar to SBS1 Aetiology: spontaneous or enzymatic deamination of 5-methylcytosine [cosine-similarity: 0.933]--Found Signature_3 most similar to SBS4 Aetiology: exposure to tobacco (smoking) mutagens [cosine-similarity: 0.919]--Found Signature_4 most similar to SBS7b Aetiology: UV exposure [cosine-similarity: 0.979]--Found Signature_5 most similar to SBS6 Aetiology: defective DNA mismatch repair [cosine-similarity: 0.719]--Found Signature_6 most similar to SBS5 Aetiology: Unknown [cosine-similarity: 0.78]------------------------------------> library('pheatmap')> pheatmap::pheatmap(mat = hnsc.og30.cosm$cosine_similarities, cluster_rows = FALSE, main = "cosine similarity against validated signatures")> maftools::plotSignatures(nmfRes = hnsc.sig, title_size = 1, sig_db = "SBS")特征富集分析上面鉴定出的特征可以进一步分配到样品上，并且使用signatureEnrichment功能进行富集分析，鉴定每一个特征的突变富集。> hnsc.se <- signatureEnrichment(maf = maf, sig_res = hnsc.sig)Running k-means for signature assignment..Performing pairwise and groupwise comparisions..Sample size per factor in Signature:Signature_1 Signature_2 194 313 Estimating mutation load and signature exposures..这一步有一点问题，上一步计算特征数量的时候，如果我用的n=6得出的hnsc.sig，在这里就会报错。我查阅了网上，也有人和我有一样的报错情况，他是把n值降低，这一步就不会报错了（signatureEnrichment not assigning all samples to signatures）,所以我试了n=2又重新运行了一遍，这里就没有出现报错信息，在这个网址里也提到，如果他使用n>3的数值也会报错，不知道是不是个bug。当然上面的图也要重新做。这里我就不再贴一遍图了。>plotEnrichmentResults(enrich_res = hnsc.se, pVal = 0.005)