{'fraction_non_alphanumeric': 0.032124645892351276, 'fraction_numerical': 0.014504249291784703, 'mean_word_length': 5.370871683811586, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 7, 'https://': 4, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Maftools: efficient and comprehensive analysis of somatic variants in cancer

Numerous large-scale genomic studies of matched tumor-normal samples have established the somatic landscapes of most cancer types. However, the downstream analysis of data from somatic mutations entails a number of computational and statistical approaches, requiring usage of independent software and numerous tools. Here, we describe an R Bioconductor package, Maftools, which offers a multitude of analysis and visualization modules that are commonly used in cancer genomic studies, including driver gene identification, pathway, signature, enrichment, and association analyses. Maftools only requires somatic variants in Mutation Annotation Format (MAF) and is independent of larger alignment files. With the implementation of well-established statistical and computational methods, Maftools facilitates data-driven research and comparative analysis to discover novel results from publicly available data sets. In the present study, using three of the well-annotated cohorts from The Cancer Genome Atlas (TCGA), we describe the application of Maftools to reproduce known results. More importantly, we show that Maftools can also be used to uncover novel findings through integrative analysis.

Numerous large-scale genomic studies of matched tumor-normal samples have established the somatic landscapes of most cancer types. However, the downstream analysis of data from somatic mutations entails a number of computational and statistical approaches, requiring usage of independent software and numerous tools. Here, we describe an R Bioconductor package, Maftools, which offers a multitude of analysis and visualization modules that are commonly used in cancer genomic studies, including driver gene identification, pathway, signature, enrichment, and association analyses. Maftools only requires somatic variants in Mutation Annotation Format (MAF) and is independent of larger alignment files. With the implementation of well-established statistical and computational methods, Maftools facilitates data-driven research and comparative analysis to discover novel results from publicly available data sets. In the present study, using three of the wellannotated cohorts from The Cancer Genome Atlas (TCGA), we describe the application of Maftools to reproduce known results. More importantly, we show that Maftools can also be used to uncover novel findings through integrative analysis.
[Supplemental material is available for this article.] With advances in cancer genomics and reduction in costs, wholegenome sequencing (WGS) and whole-exome sequencing (WXS) of large cohorts of cancer samples have become the mainstream way of determining genetic abnormalities associated with cancer (Mardis and Wilson 2009;Vogelstein et al. 2013;Wheeler and Wang 2013). Along with the many ongoing studies, a plethora of published large-scale data sets offer opportunities for reanalysis to advance our understanding of cancer genome and biology. Such cohort-based large-scale characterizations often produce large amounts of data in the form of somatic variants containing single-nucleotide variants (SNV) and small insertion/deletions (indels). Somatic variants provide baseline data for many analyses such as driver gene detection, pathway analysis, mutational signatures, and estimation of tumor heterogeneity, to name a few (Alexandrov et al. 2013a;Lawrence et al. 2013b). However, these downstream analyses of mutational data often entail many computational/statistical approaches, which are laborious, time-consuming, and cumbersome. On the other hand, visualization of these complex and heterogeneous data plays key roles in genomic studies, with researchers finding it difficult to generate complicated publication-quality images, such as oncoplots and lollipop plots. Although a number of tools and software exist for each of these tasks, they all require specific input data format . Although tools such as Mutational Significance in Cancer (MuSiC) offer multiple analysis domains in a single software package, they require large alignment files and computational resources, hindering their usage in analyzing public data sets (Dees et al. 2012).
To address these issues, we developed a user-friendly, R Bioconductor package, which we termed "Maftools." Maftools offers a multitude of analysis and visualization modules while only requiring a single input text file containing somatic variants in MAF format. MAF is a standard tab delimited text file format introduced by TCGA for storing and distributing somatic variants, containing complete somatic landscape of the cohort. MAF also has several advantages over variant call format (VCF) in storing annotations of hundreds of samples while maintaining readability and portability. With MAF files as a standard input, functions are implemented in Maftools to perform many commonly used statistical and computational analyses in cancer genome studies, including but not limited to driver gene detection and analysis of pathways, de novo signatures, and clinical parameter enrichment. Maftools also provides options to integrate and analyze copy number variation (CNV) data generated by programs such as genomic identification of significant targets in cancer (GISTIC) and circular binary segmentation (CBS) algorithms (Olshen et al. 2004;Mermel et al. 2011). In addition, Maftools can be used to perform variant annotations and other common tasks such as data format conversion and subset operations. Usage of Maftools is straightforward with self-explanatory functions and is implemented as an open source R package available through the Bioconductor project. Maftools is independent of alignment files, facilitating analysis of public data sets through data-mining approaches. Lastly, Maftools provides various plotting functions to help researchers generate intricate publication-quality images.
To demonstrate the application and performance of Maftools, we used three different TCGA cohorts: esophageal carcinoma (ESCA; WXS; N = 186), acute myeloid leukemia (AML; WXS; N = 192), and breast invasive carcinoma (BRCA; WGS; N = 96), which all contain both somatic variants as well as copy number variations (The Cancer Genome Atlas Network 2012b. ESCA is unique in its mutational pattern and its two molecularly distinct subtypes-esophageal adenocarcinoma (EAC) and esophageal squamous cell carcinoma (ESCC)-as we and others have characterized previously (Lin et al. 2014(Lin et al. , 2018aHao et al. 2016). AML was selected because it has several well-known dysregulated pathways and clinical associations. WGS data from BRCA is ideal for the demonstration of rainfall plots and identification of hypermutated genomic regions (known as "Kataegis") (Alexandrov et al. 2013b).

Maftools package overview
As described above, the Maftools package was developed to bring the majority of standard analysis and visualization modules into a single channel through implementation of well-established statistical and computational approaches, while only requiring a single and unaltered input data format.
Functions of Maftools are divided into three main modules ( Fig. 1): analysis, visualization, and annotation. Each of these modules and key functions are described below.

Visualization
Clear and concise visualization of large-scale genomic data is a key step in displaying critical information in an effective, precise, and easy-to-comprehend manner. The visualization module in the Maftools package offers various plotting functions to generate cus-tomizable and feature-rich publication-quality plots from both somatic mutation and CNV data sets. Some of the visualizations offered by Maftools include oncoplots (also known as waterfall plots, coMut plots, and oncoprints) to display somatic landscape, lollipop plots (also known as stick plots, needle plots, and stem plots) to illustrate the distribution of variants on a linearized protein structure, summary plots to summarize mutation load, transition/transversions plot, and rainfall plot to visualize Kataegis phenotype.
Using ESCA and BRCA cohorts, we show some of the key visualizations generated using Maftools. Figure 2A, known as oncoplots, displays significantly mutated genes (FDR < 0.1) identified by MutSigCV algorithm in ESCA cohort (Lawrence et al. 2013b). Genes are sorted by mutational frequency, and samples are sorted and ordered according to tumor histology, thereby differentiating the mutational spectrum between and within subtypes of ESCA.
Options are available to include annotations as bottom annotation bars to display clinical parameters. The transition and transversion plot (Fig. 2B) summarizes SNVs into six categories. Figure 2C shows a lollipop plot for a highly mutated gene, TP53, in ESCA. To keep up with the consistency of the plot design, lollipop plots generated by Maftools follow the same visual aesthetics of the commonly used online tools available as a part of cBioPortal (Cerami et al. 2012). Supplemental Figure S1 lists the rest of the visualization options available for plotting mutational data.
Multiple studies have reported hypermutated genomic regions in several cancer types (Nik-Zainal et al. 2012;Alexandrov et al. 2013a;D'Antonio et al. 2016). These genomic regions referred to as "Kataegis," presumed to be the result of aberrant activity of apolipoprotein B mRNA editing enzyme, catalytic polypeptidelike (APOBEC) cytidine deaminases, tend to be enriched in C > T and C > G substitutions occurring within TpCpN trinucleotide contexts (Lada et al. 2012 Maftools can also be used to visualize and summarize copy number data generated by GISTIC or segmentation files generated by CBS algorithms. Copy number data can be easily integrated along with mutation data, further facilitating integrative analysis. Supplemental Figure S2 displays plotting options available to visualize CNV data.

Mutational signatures and enrichment analysis
Every cancer, as it progresses, leaves behind a characteristic mutational pattern that can reveal its underlying mutagenic processes. Alexandrov et al. (2013a) have shown that such mutagenic processes can be identified by utilizing dimensional  Table headers describe available modules, namely, Visualization, Set operation, Variant annotations, and Analysis. A small description and the corresponding function (bold italics) in Maftools package are provided for every module. The typical workflow begins with MAF object creation either by reading an MAF file or by converting existing annotations to an MAF object, which is later passed to a desired function as an input (arrows). The visualization module includes functions to generate publication ready plots from the input MAF object, whereas the analysis module offers functions to perform commonly performed analyses in cancer genomics. The variant annotation module performs variant annotations using oncotator API and format conversion of annotations generated by programs such as ANNOVAR. The set operation includes function for subsetting MAF object based on user-defined queries. reduction techniques such as non-negative matrix factorization (NMF). Specifically, this method classifies SNVs into six different transition and transversion events, each of which is further classified into 16 subtypes based on immediate 5 ′ and 3 ′ bases surrounding the mutated base (Nik-Zainal et al. 2012;Alexandrov et al. 2013a). Typical de novo signature analysis includes frequency matrix generation and NMF decomposition. We implemented two functions, namely, trinucleotideMatrix and extractSignatures, to streamline the complete process. We further extended this analysis by signtaureEnrichment function to perform sample stratification, signature enrichment, and association analysis.
Furthermore, we implemented a method to perform sample classification and analysis of signature enrichment, wherein samples were assigned to identified signatures using k-means clustering of signature exposures, followed by Fisher's exact tests (Supplemental Fig. S4A,B). To test its performance, we compared APOBEC enrichment scores (measured by trinucleotideMatrix method) and APOBEC signature weight (measured by signatureEnrichment function). Notably, this comparison revealed a high correlation coefficient (Spearman's correlation coefficient: 0.9438; P-value < 2.2 × 10 −6 ) (Supplemental Fig. S4C), highlighting the concordance between these two methods. We next performed an enrichment analysis to identify genes preferentially mutated along with a particular signature. Applying this method on ESCCs revealed exclusive mutations of the TGFBR2 gene in samples associated with COSMIC Signature 6 (Fisher's exact text; P < 0.001) (  Table S3). Signature 6 has been associated with MMR deficiency and is characterized by C·G→T·A transitions at a NpCpG sequence context and C·G→A·T transversions at CpCpC. Several solid tumors associated with Signature 6-such as breast, colon, and ovarian cancers-are shown to be enriched for small indels, a characteristic feature associated with Microsatellite Instability (MSI) (Helleday et al. 2014). Also, TGFBR2 mutations have been previously shown to occur among MSI tumors in colon cancer (The Cancer Genome Atlas Network 2012a). Specifically, truncating mutations in TGFBR2 affecting BAT-RII domain, cause deficiency in mismatch repair (MMR) pathways (Biswas et al. 2008). Similar to observations in colon cancer, mutations in ecTBetaR2 domain of TGFBR2 were mostly loss of function (frame shift, indels, and splice site), further highlighting the role of TBGBR2associated MMR-deficient pathways in ESCC (Fig. 3E).
Combined observations from signature and enrichment analyses suggest that the use of therapies such as methotrexate or poly(ADP-ribose) polymerase (PARP) inhibitors targeting MMR-deficient pathways may have therapeutic activity in ESCC (Martin et al. 2010).

Cohort comparison and Pfam domain summarization
Different forms of cancers differ in their mutational burden and overall mutational landscape depending on tissue lineage and underlying mutagenic processes (Lawrence et al. 2013b). On the other hand, clinical parameters and histopathology contribute to tumor heterogeneity within a single cancer type. We implemented mafCompare function to identify such differentially mutated genes (DMGs) and pathways between two cohorts, wherein mutation load for each gene is compared by Fisher's exact tests. Comparison of ESCC and EAC cohorts revealed 38 genes to be differentially mutated (P < 0.01) ( Fig. 4A; Supplemental Table S4). Among them, only four genes (NFE2L2, TP53, LRP6, and PDS5B) were significantly enriched in ESCC, and the other 34 genes were significantly enriched in EAC, largely validating our recent reports using different cohorts (Lin et al. 2018a). These 34 EAC-specific genes were enriched in ERBB signaling and pathways associated with sodium channel signaling (Fig. 4B). ERBB signaling pathway has been known to be dysregulated in EAC more frequently than ESCC (Fichter et al. 2014; The Cancer Genome Atlas Research Network 2017; Lin et al. 2018a). Alterations of genes involved in sodium signaling pathways have not been documented in EAC. Mutant genes of this pathway included SCN3A, SCN5A, and SCN9A, which belong to the family of voltage gated sodium channels (VGSCs) involved in action potential initiation and conduction in excitable cells such as cardiac and neuronal cells (Catterall 2012). VGSC factors are aberrantly expressed in several types of cancers and contribute to cell migration and metastasis (Schönherr 2005;Fiske et al. 2006;House et al. 2015). In addition, VGSC mutations in glioblastoma are associated with shorter survival rates of the patients (Joshi et al. 2011). However, their roles in EAC are unknown. Notably, uniformly distributed activating mutations in a mutually exclusive manner ( Fig. 4B; Supplemental Fig. S5) strongly suggest that these may be gain-of-function variants.  Table S5). Interestingly, we observed that a specific domain, pfam00520 (Ion_trans), was affected more frequently in EAC (mutated 95 times across 41 genes) (Fig.  4C) compared to ESCC (mutated 35 times across 25 genes) (Fig.  4D). The Ion_trans domain is a family of six transmembrane helicases responsible for ion transportation, primarily involved in ion (Na + /Ca 2+ ) signaling and L1CAM interaction pathways (Jiang et al. 2003). Importantly, these observations in the Ion_trans domain strongly support our earlier results of mutated pathways involving VGSCs (Fig. 4B) and the potential role of dysregulated sodium channels in EAC.

Somatic interactions, cancer genes, and clinical enrichment analysis
Recent evidence from cancer genomic studies demonstrates that key genes in dysregulated pathways are often mutated in a mutually exclusive manner (Yeang et al. 2008). Identification of such gene sets can reveal de novo pathways and underlying mechanisms of tumorigenesis. Here, we developed a somaticInteractions function, which facilitates the identification of gene sets mutated in either mutually exclusive or a co-occurring manner. To demonstrate its performance, we used the TCGA AML cohort, which has well-documented gene sets mutated in either a mutually exclusive or co-occurring manner (The Cancer Genome Atlas Research Network 2013). Importantly, our method accurately reproduced a number of these observations, such as TP53/FLT3 (mutually exclusive; P-value = 0.012) and NPM1/FLT3 (co-occurring; P-value = 0.00104) ( Fig. 5A; Supplemental Table S6). Extending this method further to gene sets of larger size (N = 3) identified several pathways mutated in an exclusive manner (exact test P < 0.001), including well-known transcription factor fusion genes, NPM1, TP53, and RUNX1 ( Fig. 5B; Supplemental Table S7; The Cancer Genome Atlas Research Network 2013).
Cancer driver genes are defined by their capability of providing selective growth advantage to cancer cells when genetically altered. Several mathematical approaches have been developed to identify such driver genes, based on mutational frequencies, functional impact, or by clonality modeling (Dees et al. 2012; Gonzalez-Perez and Lopez-Bigas 2012; Lawrence et al. 2013b;Zapata et al. 2017). To facilitate such detection, we built a function oncodrive based on OncodriveCLUST algorithm, which leverages the observation that a majority of the activating mutations within oncogenes are clustered around mutational hotspots (Tamborero et al. 2013). Applying oncodrive on the TCGA AML cohort identified 11 well-known oncogenes as significantly mutated (FDR < 0.1) (Fig. 5C). However, as noted in the original study, oncodrive is biased toward oncogenes with mutational hotspots and has less optimum performance in identifying potential tumor suppressors, such as TP53, with randomly distributed mutations across the gene body. Despite this, the majority of oncodrive results overlapped with driver genes identified by the widely used program MutSigCV (9 of 11) (Supplemental Fig. S6A; Lawrence et al. 2013b). The two drivers identified exclusively by oncodrive are known leukemic oncogenes, KIT and KRAS (Supplemental Fig.  S6B,C). These results underscore the performance of oncodrive and suggest that combined approaches of multiple programs can increase the sensitivity of identifying driver genes. Via the pancanComparison function, users can analyze MutSigCV results and compare them to pan cancer driver genes identified by analyzing more than 4000 samples across 21 different cancer types (Lawrence et al. 2014). Such comparison identifies common driver genes as well as tumor-specific ones. Applying pancanComparison on ESCA cohort revealed several common driver genes, including TP53 and CDKN2A. Importantly, tumorspecific driver genes, such as TGFBR2 and ZNF750, were uniquely mutated in ESCC, validating the findings from us and others ( Fig.  5D; Lin et al. 2014).

Variant annotations, format conversions, and subset operations
Maftools also includes functions to perform quick variant annotations and format conversions. Oncotate function takes raw variants stored in a simple tabular format and annotates them using oncotator's REST web API. However, often this process is relatively slower and time-consuming for larger inputs (most time is needed for connecting to API and retrieving annotations). Another widely used tool for variant annotations is Annotate Variation (ANNOVAR), capable of annotating a putative variant with several gene, region, and filter-based annotations. Tabular output files generated from ANNOVAR can be converted to MAF with function annovarToMaf, which parses and maps values from gene-based annotations into MAF-specific values. For converting VCF files to MAF, we recommend vcf2maf utility (https://github.com/mskcc/ vcf2maf), which processes variant annotation (with Variant Effect Predictor) and transcript prioritization.
Similar to TCGA, ICGC is an international consortium providing large-scale cancer genome data. However, somatic variants from ICGC are in "simple somatic mutation format," a standard tab delimited text file introduced by the ICGC consortium to store and distribute somatic variants. We implemented icgcSimpleMuta-tionToMAF function, which converts ICGC mutation format files Significantly altered pathway identified in AML by CoMEt exact test involving NPM1, RUNX1, and TP53 genes mutated in mutually exclusive manner (P < 0.001). (C ) Disease-associated driver genes identified by oncodrive in TCGA AML cohort (FDR < 0.1). Number of closely spaced mutational clusters are highlighted within brackets. (D) Pan cancer comparison of significantly mutated genes in ESCA identified by MutSigCV against Pan cancer driver genes. TGFBR2 and ZNF750 are exclusively mutated in esophageal cancer, whereas other drivers, such as TP53, and CDKN2A, are mutated in the global pan cancer cohort. (E,F) Bar plots displaying the association between genes and clinical features, French-American-British (FAB) classification and age group, respectively (P < 0.05, Fisher's exact test). Bars are annotated with the ratio of mutated samples to total samples. Error bars display 95% CI of binomial ratios. The y-axis denotes the fraction of samples associated with the phenotype.
into MAF, thereby streamlining ICGC data processing. subsetMaf is another function in Maftools that allows subset operations of MAF files based on genes of interest, samples, or user-specified queries.

Discussion
With the rapid increase in genomic analysis, exponential growth has occurred in the availability of software to reanalyze and understand these data . MuSiC is one of the most widely used tools that offers several analysis modules to perform tasks such as identification of significantly mutated genes, pathway analysis, and clinical correlation tests (Dees et al. 2012). However, its application is rather limited by its dependency on alignment files, platform specificity, and large computational requirements. This further limits its application on public data sets, which are often restricted by the availability of raw alignment files. Several tools such as EmU and deconstructSigs offer statistical frameworks to identify de novo signatures and enrichment of known signatures (Fischer et al. 2013;Rosenthal et al. 2016). Other tools such as Dendrix, PathScan, and HotNet2 can be used for pathway and network analysis (Wendl et al. 2011;Vandin et al. 2012). However, these tools have different requirements of data formats and much time is required in preparing data sets and file format conversions, significantly hindering the efficiency of cancer researchers.
Here, we introduce an R Bioconductor package, Maftools, which integrates standard analyses modules that are frequently performed in cancer genomics and provides a plethora of visualization options to generate publication-quality images. Along with the analysis of somatic variants, Maftools allows easy integration of copy number data generated from either GISTIC or segmentation algorithms. Maftools utilizes well-established dimensional reduction and statistical methods, enabling reproducible datadriven research with the requirement of few lines of code. In addition, inclusion of clinical data can identify novel clustering and enrichment patterns. Moreover, Maftools is easy to use, platform independent, and does not rely on large alignment files, greatly facilitating the exploration of genomic data sets such as those from TCGA and ICGC.
Here, we reproduced many known results utilizing only MAF files from the published TCGA data sets. More importantly, we showed that Maftools can also be used to uncover novel findings through integrative analysis. Via implementation of wellestablished computational and statistical methods, Maftools provides a wide range of functions for cancer genomic analyses. In future updates, we will include gene expression as well as DNA methylation data for integrative multiomic analysis.

Data sets
TCGA MAF files (ESCA: http://dx.doi.org/10.7908/C1BV7FZC; LAML: http://dx.doi.org/10.7908/C1D21X2X) along with the clinical data, MutSigCV (v1.41), and GISTIC2 results for ESCA and AML cohort were obtained from Broad Firehose using fireho-se_get utility (analysis stamp: 2016_01_28). Somatic mutations for BRCA WGS samples were obtained from a published study (D'Antonio et al. 2016) and were later converted to MAF using vcf2maf utility. All raw input data and reproducible R code used to generate the results are provided as Supplemental Data S1. Computational time required for each function is provided in Supplemental Table S8.

R/Bioconductor package
Maftools is implemented as an open source R package and available as a part of the Bioconductor project (Gentleman et al. 2004). MAF file and the associated clinical data along with summary statistics are stored in an S4 class container. For faster data processing and summarization of larger data sets, Maftools uses data.table library, which offers performance of several magnitude higher than base R functions. Package workflow is simple, with every function taking MAF object as an input, and comes with a detailed vignette including a case study describing the usage of available modules.

Visualization
Visualization module in Maftools facilitates generation publication-quality images with easy to use and customizable functions. Plots such as oncoplots and oncostrips are generated using ComplexHeatmap Bioconductor package, whereas plots generated to display results from analysis modules, such as forest plots, lollipop plots, somatic interactions, rainfall plots, among others, are generated using either ggplot2 or base R plot functions (Gu et al. 2016).

Signature and enrichment analyses
Mutational signature analysis in Maftools begins with the process of extracting 5 ′ and 3 ′ adjacent bases surrounding the mutated base, thereby constructing a count matrix M of dimension 96 × n, where n is the number of samples available in input MAF file. This process is implemented in the function trinucleotideMatrix which uses BioStrings and GenomicRanges Bioconductor packages for efficiently reading reference genome and extracting adjacent bases while simultaneously classifying mutations into transition and transversion events (Lawrence et al. 2013a).
Once the matrix is generated, the extractSignatures function uses NMF to factorize count matrix M into two smaller matrices -W (96 × r) and H (r × n)-such that product of W and H sufficiently recomposes the original matrix M (Gaujoux and Seoighe 2010).
A key step in factorization is identifying optimal rank, r, used to approximate the target matrix M. Several methods have been described in recent studies such as Bayesian, and Expectation Maximization-based approaches to estimate the ideal value of r (Fischer et al. 2013;Kim et al. 2016). Here, we use the method described by Brunet et al. (2004), wherein NMF is run on a range of incremental values of r, and for each value, the cophenetic correlation coefficient (measure of goodness of fit) is calculated. A final optimal r is chosen such that further increase in r results in decreasing values of the coefficient (Supplemental Fig. S3A,B). Signatures identified following matrix factorization are scaled and compared to known mutagenic processes (COSMIC signatures), and a cosine similarity value is estimated for the best possible match (Alexandrov et al. 2013b). For signature enrichment analysis, we use matrix H, containing signature exposures for every sample in every signature. Using k-means clustering, we group the samples into r clusters, thereby assigning samples to an identified signature.
For APOBEC-based enrichment analysis, we integrated the method described by Roberts et al. (2013) to estimate an enrichment score, which defines the strength of APOBEC-related mutagenic processes for every tumor sample in MAF. Briefly, enrichment of C > T mutations occurring within tCw trinucleotide context over all of the C > T mutations in a given sample is compared to background cytosines and tCw occurring around ±20 bp of mutated bases. We further extended this method to identify genes associated with the APOBEC enrichment by classifying samples as APOBEC-enriched (enrichment score >2) and non-APOBEC-enriched (enrichment score <2), followed by one-way Fisher's exact tests to identify genes overrepresented among APOBEC-enriched samples.

Somatic interactions and pathway analysis
Somatic interactions function in Maftools allows users to identify gene sets mutated in a mutually exclusive or co-occurring manner. For a pair of genes, the pattern of exclusiveness or co-occurrence is estimated by performing a Fisher's exact test on a 2 × 2 contingency table containing frequencies of mutated and nonmutated samples. However, Fisher's exact test is limited by degrees of freedom for larger contingency tables (n > 2; with 2 n -n-1 degrees of freedom). To identify pathways involving gene sets (n > 2) mutated in an exclusive manner, we included Combinations of Mutually Exclusive Alterations (CoMEt) algorithm implemented in the cometExactTest R package (Leiserson et al. 2015). CoMEt exact test takes a contingency table of 2 n and computes the probability for mutual exclusiveness. As an input to the CoMEt exact test, we first identify all pair of genes that are mutated in a mutually exclusive manner using Fisher's exact test, and CoMEt exact test is run on all unique combination of gene sets of size n. For example, given a gene list of six genes and gene set of size three, the CoMEt exact test is performed on 20 combinations (6C3) to identify gene sets mutated in an exclusive manner. Options are included for users to manually specify the desired gene set size.

Identification of cancer genes
To identify disease-associated cancer genes, we reimplemented OncodriveCLUST algorithm as previously described by Tamborero et al. (2013). Briefly, candidate amino acid positions with mutations above a background threshold (based on binomial distribution accounting for protein length and mutation rate) are compiled for every gene. Candidate positions occurring within a distance of five amino acids are grouped together to form a cluster. Clusters are further refined and extended by merging neighboring amino acid positions occurring within five amino acids. Once cluster formation is completed, a cluster score is calculated, which is a ratio of mutations occurring within the identified clusters to the total number of mutations. Finally, a P-value is estimated based on t-statistic and Z-score.
OncodriveCLUST has originally been implemented in Python framework requiring users to prepare input data containing gene symbols and amino acid positions. Its reimplementation as a part of the Maftools package serves two purposes. The package takes care of file parsing and input data preparation allowing novice users to run the program easily. Its implementation in R programming language exposes it to the larger bioinformatic community and also allows visualization of results in an intuitive manner (Fig. 5C).

Cohort comparison and enrichment analysis
mafCompare function in Maftools allows comparison of two independent cohorts to identify differentially mutated genes or to perform association between clinical features. A 2 × 2 contingency table of frequencies of mutation is calculated for every gene from the input cohorts followed by Fisher's exact test to identify genes showing significant differences in their mutation frequencies.
Similarly, for clinical enrichment analysis, once again contingency tables are generated for every categorical variable followed by Fisher's exact test to calculate P-values. Results from cohort comparison and enrichment analysis are visualized as forest plots or frequency bar plots (Figs. 4A, 5E,F).

Pfam domain summarization
Inspired by the Pfam annotation module of MuSiC, we implemented pfamDomains function to identify and group mutations by affected protein domains. However, the approach for grouping mutations is much simpler. MuSiC requires a large database containing protein foci translated to genomic loci, which is later queried by external tools such SAMtools tabix to identify affected protein domains. Here, we use protein change information (in HGVSp format) to parse the protein position and transcript annotations. These positions are later mapped onto the Pfam domain database via rapid data.table foverlaps function for summarization. Plotting options are available to display frequently affected domains by means of bubble plots (Fig. 4C,D).

Change-point detection for identification of hypermutated genomic regions
Hypermutated genomic regions are segments along the chromosome where the mutation rate is significantly higher than the average mutation rate across the genome. To identify such regions, we utilize change-point detection method in which genomic boundaries with the sharp decrease in distance between consecutive mutations are identified. Briefly, mutations are ordered for every chromosome based on the position, and distance between consecutive mutations are calculated followed log transformation. Log transformed inter-event distances are later fed into cpt.mean function implemented in the changepoint R package to identify potential change points (Killick and Eckley 2014). Consecutive change points are merged into genomic segments and annotated as Kataegis if the segment contained six or more consecutive mutations with an average inter-mutation distance of <1000 bp.

Software availability
Maftools is implemented as an R package, released under MIT license. Source code is available on GitHub (https://github.com/ PoisonAlien/Maftools) and can be installed via Bioconductor project (https://bioconductor.org/packages/release/bioc/html/maftools .html) (R Core Team 2018). R package source code is also available in the Supplemental Data S2. In addition, we also provide an R data package containing ready-to-use, precompiled somatic variants from Broad Firehose, and Multi-Center Mutation Calling in Multiple Cancers project for all 34 TCGA cohorts along with the relevant clinical information (https://github.com/PoisonAlien/TCGAmutations) (Ellrott et al. 2018).