Exome-Scale Discovery of Hotspot Mutation Regions in Human Cancer Using 3D Protein Structure

Missense mutations are perhaps the most difficult mutation type to interpret in human cancers. Truncating loss-of-function mutations and structural rearrangements generate major changes in the protein product of a gene, but a single missense mutation yields only a small change in protein chemistry. The impact of missense mutation on protein function, cellular behavior, cancer etiology, and progression may be negligible or profound, for reasons that are not yet well understood. Missense mutations are frequent in most cancer types, accounting for approximately 85% of the somatic mutations observed in solid human tumors (1), and the cancer genomics community has prioritized the task of identifying important missense mutations discovered in sequencing studies. Whole-exome sequencing (WES) studies of cancer have created new opportunities to better understand the importance of missense mutations. This enormous collection of data now allows detection of patterns with power that was unheard of a few years ago.

The first approaches to identify cancer drivers from WES mutations looked for significantly mutated genes (SMG), harboring a larger number of somatic mutations than expected by chance (2–5). Metrics to call SMGs now appear to be underpowered given the size of current cohorts in The Cancer Genome Atlas (TCGA) and International Cancer Genome Consortium (ICGC). A recent study suggested that approximately 1,500 cases of endometrial cancer would need to be sequenced to attain 90% power to detect mutations in 90% of genes with a mutation frequency of 2% with the SMG approach (5). The recognition of the limitations of the SMG paradigm has motivated interest in orthogonal analysis techniques to detect mutational patterns associated with drivers (1, 6–9).

Recurrence of somatic missense mutations in cancers at the same amino acid residue position is well known to be a characteristic feature of both oncogenes (OG) and tumor suppressor genes (TSG; ref. 10). The observation that somatic mutations also frequently occur in positions proximal in protein sequence to the most highly recurrent positions has suggested that positional clustering of somatic missense mutations might be used to identify drivers (7). These clusters, known as "hotspots," are regions where somatic missense mutations occur closer together in protein sequence than would be expected by chance. Hotspot regions can be rationalized as areas in a protein under positive selection in the cancer environment; missense mutations occurring in these regions are selected for because they alter protein function in a manner advantageous to the cancer cell. Several algorithms have been developed to identify protein functional domains and genes in which these regions are enriched (8, 11), and to identify specific missense mutations in hotspots (7, 12). These algorithms consider mutations in a coordinate system based on one-dimensional (1D) protein sequence. WES studies of cancer cohorts are increasingly incorporating missense "hotspot" detection as a routine analysis step in the search for new drivers.

Finding missense mutation hotspot regions in 1D is limited by the fact that functional proteins tend to fold into three-dimensional (3D) structures (with the exception of intrinsically disordered regions). Thus, positional clustering done in 1D will likely miss many hotspots that are present in 3D after folding. Gene- and protein domain-level testing may indicate the possibility of a 3D hotspot but cannot identify the specific positions in the hotspot. An algorithm that leverages 3D protein structure information, but still performs clustering in 1D through a dimensionality reduction step, has shown utility in detecting OGs (9). A recent study of an aggregated collection of TCGA cancer mutations from 21 tumor types presented an algorithm to identify cancer genes based on 3D clustering of somatic missense mutations, yielding ten such genes. They reported low correlation between 3D and 1D hotspot regions (13).

Here, we present HotMAPS (Hotspot Missense mutation Areas in Protein Structure), a new, sensitive algorithm and a web-based community resource for high-throughput analysis of cancer missense mutation 3D hotspot regions. HotMAPS finds clusters of amino acid residues with significantly increased local mutation density in 3D protein space, compared with an empirical null distribution. The statistical model is designed to handle higher-order protein complexes and can capture regions that span protein–protein interfaces. We apply HotMAPS to missense mutations from 23 tumor types sequenced by TCGA. By careful use of both experimentally derived protein biologic assemblies in the Protein Data Bank (PDB) and theoretical protein structure models, we substantially increase the number of amino acids that can be mapped into 3D protein space and the number of detectable hotspot regions (13).

HotMAPS systematically delineates 3D hotspot regions on the level of amino acid positions, and we provide a detailed catalog of 216 tumor-type–specific regions. We show how the catalog can be used as a discovery tool so that the links between 3D protein structure and the biologic functions of missense mutations in cancer can be better utilized by the community. The catalog provides comprehensive identification of hotspot regions that overlap with many key biologic features of proteins available in the literature (e.g., residue positions at active sites, small-molecule and metal-binding sites, protein interfaces, positions with published experimental mutagenesis results). This information can potentially provide a researcher with more fine-grained mechanistic understanding of missense mutation cancer relevance than is possible by 1D clustering or domain and gene enrichment approaches. Using the catalog, we were able for the first time to systematically analyze characteristic properties of 3D hotspot regions and differences between 3D hotspot regions in OGs and TSGs.

TCGA mutation annotation format (MAF) file data for 23 tumor types was downloaded from Xena data store (https://genome-cancer.soe.ucsc.edu/proj/site/xena/hub/) using their API.

PDB structures were obtained from the Worldwide Protein Data Bank (10/17/2015). Only structures solved by x-ray crystallography and containing at least one human protein chain were used. Single-domain, theoretical protein structure models constructed on the basis of homology to nonhuman proteins were included to increase coverage over a greater proportion of genes. Theoretical models were obtained from the ModPipe human 2013 dataset (ftp://salilab.org/databases/modbase/projects/genomes/H_sapiens/2013/), built with Modeller 9.11 (14). In addition to criteria required by ModPipe, we filtered the theoretical models to increase the quality of structures used in our assessment based on minimum length, target-template sequence identity, loop content, and radius of gyration (Supplementary Materials and Methods).

Models were assessed by comparing 3D hotspot regions identified by HotMAPS in experimental structures with those identified in theoretical models of the same protein. First, we found all pairs of experimental structures and theoretical models of the same protein, in which there was overlap of the same amino acid residues. The agreement of a structure/model pair was the overlap of their hotspot region–mutated residues. A false-positive error was called when a model had a mutated residue in a hotspot region that was not in a hotspot region for any protein structure that it had been paired with. A false-negative error was called when a structure had a mutated residue in a hotspot region that was not in a hotspot region for any of the models it had been paired with.

HotMAPS identifies residue positions with higher local mutation density in each protein structure or model than expected from an empirical null distribution, based on simulations of a discrete uniform distribution. Residues are considered significant for increased local mutation density at FDR threshold of 0.01, after correction for multiple testing (Benjamini Hochberg). Three-dimensional missense mutation hotspots are identified as groupings of significant residues according to the principle of maximum parsimony, based on connected components in a neighbor graph. Construction of the neighbor graph and connected components are illustrated in Supplementary Fig. S1. HotMAPS is designed to run on both single-chain protein structures and biologic assemblies with multiple chains originating from the same gene. Mathematical details are provided in Supplementary Materials and Methods.

Applying HotMAPS to 19,368 PDB protein structures (PDB bioassemblies in which in vivo protein structure is represented) and 46,004 theoretical models, we identified 107 unique 3D mutation hotspot regions (aggregated across tumor types), of which 30 were only detectable by clustering in 3D (Supplementary Table S1). When stratified by tumor type, 216 3D missense mutation hotspot regions were found in 19 of the 23 TCGA tumor types, with none in adrenocortical carcinoma (ACC), kidney renal papillary cell carcinoma (KIRP), liver hepatocellular carcinoma (LIHC), or kidney chromophobe (KICH; Supplementary Table S2). KICH is known to be driven by alterations other than point mutation, such as structural breakpoints in the TERT promoter. Among all 23 tumor types, sample and mutation count is lower for these four tumor types (P = 0.02 for sample count, P = 0.04 for mutation count; Wilcoxon rank-sum test), suggesting that at least for tumor types driven by missense mutations, larger sample size might increase our power to find more 3D regions. Our approach enabled us to consider the 3D protein environment of a much higher fraction of TCGA mutations than has been described previously. We were able to map and analyze approximately 53% of the missense mutations in 23 TCGA tumor types (Supplementary Table S3). Of these missense mutations, approximately 10% could be mapped to PDB protein structures and an additional 42% mapped to theoretical models, in the absence of PDB structure. Using hotspot regions identified in the PDB structures as a control, we estimate that the hotspots called in the models have a false-positive rate of 0.058 and a false-negative rate of 0.138. Therefore, very few hotspot regions found in the models are the result of modeling errors, justifying the increase in mutation coverage obtained.

When hotspot regions were stratified by tumor type, 91 genes contained at least one hotspot region in at least one tumor type (q ≤ 0.01) and 40 of these genes had regions that were only discoverable by consideration in 3D. Of the 91, 19 genes were previously annotated as OGs and 11 genes as TSGs with the 20/20 rule, a ratiometric method based on the proportions of different mutation consequence types observed in a gene (1). Twenty-five of the genes are listed in the Cancer Gene Census (CGC; Table 1; ref. 15). Of the remaining 58 genes, five (KLF5, SMARCA2, RASA1, TGBFR2, KEAP1) have been identified as candidate TSGs in the literature (16–24), and six (deacetylated KLF5, MAPK1, FSIP2, RANBP2, MTOR, CHEK2) as candidate OGs (16, 25–30) . Four of the genes are current or potential drug targets (SMARCA2, HDAC4, PARG, HLA-A; refs. 18, 31–33). Two genes (ERCC2, CHEK2) are involved in hereditary cancer susceptibility when mutated in the germline (1). GTF2I is a prognostic biomarker in thymic epithelial tumors (genes with literature support in Table 2, other genes in Supplementary Table S4; ref. 34).

Among the set of genes with available protein structure or models (n = 15,697), the genes harboring a 3D hotspot region are enriched for OGs and TSGs (P = 6.1E−30 for OGs and P = 2.4E−13 for TSGs; one-tailed Fisher exact test). They are also enriched for genes in the CGC list (P = 1.4E−30; one-tailed Fisher exact test). The subset of these genes harboring only a 3D hotspot region not detectable in 1D is also significantly enriched (P = 4.3E−09 for OGs, P = 7.9E−12 for TSGs, P = 8.0E−11 for CGC genes; one-tailed Fisher exact test). An additional 23 genes that are proposed OGs, TSGs, and/or drug targets or hereditary cancer genes contained at least one 3D hotspot region. This enrichment of known and candidate driver genes supports our claim that many of the regions are biologically relevant and not simply artifacts. While regions were detected in only approximately 18% of established cancer genes, we expect that many of these genes harbor drivers other than missense mutations, some are drivers in tumor types not represented in our study and many lack structural coverage.

We examined whether the amino acid residue positions and the missense mutations in the 3D hotspot regions had distinctive features suggestive of a special biologic importance, when compared with the remaining mutations in our study. Four candidate distinguishing features were tested: (i) vertebrate evolutionary conservation; (ii) occurrence at a protein–protein interface, which increases the potential for a missense mutation to disrupt protein–protein interactions; (iii) in silico cancer driver scores generated with the CHASM algorithm (6); and (iv) in silico pathogenicity scores generated with the VEST algorithm (35), which are predictors of increased missense mutation impact (Fig. 1). In comparison with mutated residues not in 3D hotspot regions, vertebrate evolutionary conservation was higher and protein–protein interface occurrence was higher in the 3D hotspot regions (conservation P = 2.9E−29, Mann–Whitney U test; protein interface P = 5.2E−13, one-tailed Fisher exact test). In silico driver scores and pathogenicity scores were higher for missense mutations in 3D hotspot regions (driver score P = 3.0E−47, pathogenicity score P = 3.0E−16; Mann–Whitney U-test) than for the remaining mutations (Fig. 1).

The catalog contains 37 regions stratified by tumor type in bonafide TSGs and 77 in bonafide OGs (114 regions in 30 genes), using as a benchmark the classifications of Vogelstein and colleagues (landscapes benchmark; ref. 1). We used these data to explore possible differences between TSG and OG regions at amino acid resolution. We found that in TSGs, 3D hotspot regions were larger than in OGs (region size P = 9.6E−06; Mann–Whitney U test). They were also more mutationally diverse (mutational diversity P = 2.1E−07; Mann–Whitney U test). In addition, OG 3D hotspot regions were more conserved in vertebrate evolution than TSGs and more solvent accessible in protein structure, meaning that they tend to occur at the protein surface (evolution P = 4.7E−07, solvent accessible P = 1.5E−06; Mann–Whitney U test). TSG hotspot regions harbored increased mutation net change in hydrophobicity (P = 3.3E−07; Mann–Whitney U test) and mutation net change in volume (P = 2.2E−07; Mann–Whitney U test), suggesting that their impact on protein function could be due to decreased stability. The in silico missense mutation cancer driver scores were higher for OG regions (P = 0.003; Mann–Whitney U test). We also tested differences between in silico pathogenicity scores and occurrence at protein–protein interfaces between OG and TSG regions, but these were not significant (pathogenicity scores P = 0.37, protein interface P = 0.34; Mann–Whitney U test).

The fact that these differences between OG and TSG regions were statistically significant suggested that they might have predictive value. Principal components analysis (PCA) of the six significant features indicated some separation (Fig. 2A). Next, we trained a Naive Bayes machine learning classifier to discriminate between OG and TSG hotspot regions, using region size, mutational diversity, vertebrate conservation, residue solvent accessibility, mutation net hydrophobicity change, and residue volume change as features. A rigorous gene-level holdout protocol was used to avoid overfitting (Supplementary Materials and Methods). A Naive Bayes score closer to 1.0 indicates that the hotspot region is likely in an OG while a score closer to 0.0 indicates that it is in a TSG. Area under receiver operating characteristic (ROC) curve or AUC, a standard measure of classifier performance, was 0.84 out of 1.0, a result that supports our claim that 3D hotspot regions in OGs and TSGs have distinctive characteristics (Fig. 2B). AUC of a classifier with random performance is 0.5. Performance did not improve when the other features were included in the classifier. Table 3 lists the 30 genes and the median Naive Bayes score, across all regions in each gene. The median values for each predictive feature are also shown.

The ROC performance and PCA plot support our claim that characteristic differences between OG and TSG hotspots can be quantified. However, some hotspot regions remain misclassified, according to their labels in the landscapes benchmark (Discussion).

The larger size and mutational diversity of hotspot regions in TSGs versus OGs suggest that they could be more difficult to detect and perhaps they have been underreported by 1D approaches. OG hotspot regions consisting of recurrent missense mutations at one or two residues can be seen by eye with lollipop plots and are straightforward to detect computationally based on 1D primary sequence. We hypothesized that detection of many TSG hotspot regions might require a 3D algorithm. To maximize the interpretability of this analysis, regions that occurred in multiple tumor types were merged so that each region was represented only once in each gene (Materials and Methods).

For a well-controlled comparison of 3D and 1D hotspot region detection, we applied a 1D version of our method to the protein chain sequences of the same set of PDB protein bioassemblies and theoretical protein structure models to detect nonuniform clustering patterns on primary protein sequence (Supplementary Materials and Methods, Supplementary Tables S5 and S6). Seventy-two percent of hotspot regions identified in 3D were identifiable in 1D.

Next, we compared the number of hotspot regions identified in OGs and TSGs. We considered regions identified in 3D only, in both 3D and 1D, and in 1D only. Using the bona fide OGs and TSGs (Table 1), there were significantly more OG regions that TSG regions identified by the 1D algorithm (P = 0.03; one-sided Fisher exact test). The 1D-only version of the algorithm detected 5 OG and 2 TSG regions; 1D further detected an additional 25 OG and 7 TSG regions that were also identified by the 3D algorithm. The 3D algorithm identified an additional 4 OG and 6 TSG regions. To increase our power, we repeated this test again using the bona fide OGs and TSGs plus additional regions in five candidate OGs and TSGs reported in the literature (OGs were FSIP2, MTOR, RANBP2, CHEK2, and MAPK1; TSGs were RASA1, SMARCA2, KEAP1, CUL1, TGFBR2; all are listed and cited in Table 2), yielding increased statistical significance (P = 0.009, one-sided Fisher exact test). The results suggest that 1D detection methods may be better suited to detecting regions in OGs rather than TSGs.

A further problem with sequence-based 1D hotspot region detection is that larger regions detectable in 3D may be only partially characterized and/or split into multiple pieces. Figure 3 shows an example of a TSG hotspot region in FBXW7 found in 3D by HotMAPS that has been split into two pieces by the 1D algorithm. In 1D protein sequence, residue 465 is not close enough to residues 502 and 505 to be identified in one hotspot region. On the 3D protein structure of FBXW7 (PDB code 2OVQ), the three residues are spatially close and a single hotspot region is detected.

Three-dimensional consideration of hotspot regions in protein structure can potentially provide researchers with a rich source of hypothesis generation about driver mechanisms. While gene- or domain-level mutation enrichment analysis can point to potential protein functions, interactions, biologic processes, and pathways important for cancer etiology and progression, more detailed information may be available once a specific set of mutated amino acid residues has been identified as significant.

For many of the 3D hotspot regions found by HotMAPS, the literature contains evidence that they are in direct contact with or proximal to amino acid residues of known functional importance. Figure 4 shows six cancer-associated proteins in which the hotspot region is either overlapping or proximal to important functional sites.

RAC1 is a Rho GTPase important in signaling systems that regulate the organization of actin cytoskeleton and cell motility. The hotspot overlaps the GTP/GDP-binding site and could impact regulation of normal RAC1 cycling between GTP- and GDP-bound states (Fig. 4A). It contains a previously identified recurrent mutation in melanoma (P29S), which dysregulates RAC1 by a fast cycling mechanism (36).

SPOP is the substrate recognition component of a cullin3-based E3 ubiquitin–protein ligase complex, which targets multiple substrates for proteasomal degradation. The hotspot overlaps with a binding groove harboring five residue positions (pink) where mutagenesis has strongly reduced affinity for the substrate (annotated in the UniProtKB).

ERCC2 is an ATP-dependent helicase that is part of the protein complex TFIIH involved in RNA polymerase II transcription and nucleotide excision repair (NER). We identified a hotspot region, proximal to the DEAH box, a highly conserved motif containing residues that interact with Mg²⁺ and are critical for ATP binding and helicase activity (Fig. 4C). This proximity suggests that the hotspot mutations could disrupt ATPase activity and yield defective NER (37).

PTEN is a phosphatase for both proteins and phosphoinositides, and it removes a phosphate from PIP3, critical for signaling to AKT. The hotspot region identified in endometrial cancer (UCEC) spans two functionally important loops in the protein (P and WPD loops) at the boundary of the active site pocket. Residues in these loops are critical for catalysis (blue dot) and are important for the P-loop's conformation. Mutagenesis of residues in the WPD loop reduces phosphatase activity and increases colony formation in cell culture (38). Pink dots show residues that impact phosphatase activity.

RHOA is a small GTPase oncogene, and like RAC1 is a member of the Ras superfamily (39). We identified hotspot regions in bladder cancer (BLCA), head and neck squamous cell cancer (HNSCC), and stomach adenocarcinoma (STAD). The hotspot regions overlap with the RHOA effector region, a highly conserved motif that is involved in Ras superfamily signaling with downstream effector proteins. The regions are immediately proximal to a magnesium ion, which has been implicated in regulating the kinetics of Rho family GTPases (40).

VHL is a component of an E3 ubiquitin protein ligase complex, and it ubiquitinates the OG transcription factor HIF1A, targeting it for proteasomal degradation (41). One impact of VHL loss of function with failure to ubiquitinate HIF1A is increased protein expression of HIF1A. The hotspot region is proximal to its interaction site with HIF1A and could potentially have an impact on this interaction (Fig. 4F). The TCGA kidney cancer (KIRC) samples were stratified on the basis of their missense mutation status: VHL hotspot, non-hotspot, or no missense (WT). HIF1A protein expression was not significantly different between VHL non-hotspot and VHL WT groups (P = 0.5; Mann–Whitney U test), but was significantly higher between VHL hotspot and VHL WT groups (P = 0.03; Mann–Whitney U test). This result is consistent with a special role for VHL hotspot missense mutations in regulating HIF1A protein expression. However, increased HIF1A expression in these KIRC samples is likely impacted by additional genetic and other factors. We might see a substantially lower P value if VHL hotspot mutations were the only cause of the observed increase. Also, there are many VHL missense mutations outside of the hotspot region, and it is likely that several of these also have a functional impact. In particular, several of them are at the interface of VHL and the TCEB1 and TCEB2 in the complex and could impact VHL/TCEB binding.

The large-scale WES and mutation calling efforts of the TCGA have identified hundreds of thousands of somatic missense mutations in human cancers. Although some of these mutations are private, many are shared across multiple patients and multiple tumor types. The biologic and therapeutic relevance of these shared mutations is of great interest to the cancer research community. For example, patients can be stratified for clinical trials and treatment protocols selected on the basis of missense mutation status in several key driver genes, including BRAF, KRAS, and EGFR. A special type of shared missense mutations are those that occur recurrently not only at the same genomic codon, but at neighboring codons in translated protein sequence and more generally, neighboring amino acid residues in protein 3D structure. These clusters of neighboring missense mutations are known as missense mutation hotspot regions. They have been proposed to have particular relevance to oncogenic processes (12), as the increased frequency of missense mutation at a highly localized region in a protein may be a signature of positive selection (42). Missense hotspot regions may be informative in detecting driver genes (7). A number of groups have developed algorithms to detect enrichment of these regions on the gene- and domain-level (7–9, 11–13), but until now, there have been not been systematically characterized on a large number of protein structures and models at the resolution of individual amino acid residues.

We systematically identify 3D missense hotspot regions using TCGA somatic mutation data from 6,594 samples in 23 tumor types. HotMAPS identified 107 unique, tumor-type–aggregated gene-level regions, and 216 unique tumor-type–specific gene-level regions (Materials and Methods; Supplementary Tables S1 and S2). This catalog enables assessment of how the specific missense mutations in a hotspot contribute to cancer-associated molecular mechanisms.

We compared features of residues in 3D hotspot regions to other missense mutated residues in the TCGA data. The hotspot regions have characteristic features that support their putative functional importance: high evolutionary conservation, high in silico-predicted missense mutation impacts, and increased frequency of occurrence at protein interfaces. Genes containing the 3D hotspot regions appear to be particularly relevant to cancer. Landscapes benchmark TSG and OGs are overrepresented and the list includes many candidate TSGs, OGs, drug targets, and hereditary cancer genes (Tables 1 and 2). For several TSGs and OGs, the regions coincide with enzymatic active sites, positions that have been shown to alter protein function in experimental mutagenesis assays and sites of interaction with protein and nucleotide interaction partners (Fig. 4; Supplementary Table S7).

Although recurrent missense mutations have long been known to occur in both OGs and TSGs (10), they have been observed more frequently in OGs. We show that there are systematic differences in hotspot regions found in OGs and TSGs. OG regions are smaller, less mutationally diverse, more evolutionarily conserved, and more solvent accessible than TSG regions. TSG regions are more likely to harbor mutations that may impact protein stability through changes in hydrophobicity or volume. Potential explanations for these differences are that there are more ways to lose the function of a protein than to gain function (43). Loss-of-function tumor suppressor mutations can occur at many residue positions and involve many types of amino acid residue substitutions, while oncogene mutations will occur at a few functionally important positions and involve fewer substitution types.

A consequence of these differences is that TSG regions are harder to detect visually or by 1D clustering approaches than OG regions. Thus, they have been missed by 1D analysis methods. A major contribution of 3D analysis is enabling detection of hotspot regions in TSGs in addition to those in OGs. We suspect that novel hotspot regions in TSGs will continue to be discovered as more samples are sequenced in more tumor types.

We are able to leverage the characteristic differences to distinguish between hotspot regions in TSGs and OGs, with a simple machine learning method, achieving an area under the receiver operating characteristic curve (ROC AUC) of approximately 0.80. However, not all regions are correctly classified by this method. Interestingly, we find that some of these "undistinguishable" genes may act as both TSGs and OGs, depending on context or be atypical of their class. PIK3R1 has been described as an OG (44) and SPOP as a TSG (45), in agreement with our Naive Bayes scores, but not with the landscapes benchmark. The OGs IDH1 and IDH2 both have high net hydrophobicity changes, which are protein destabilizing and characteristic of TSGs. IDH1/IDH2 hotspot mutation may cause a (TSG-like) partial loss of enzymatic function, yielding accumulation of 2 hydoxyglutarate (2HG), a carcinogenic catalytic intermediate (46). EGFR has two regions in GBM and LGG, which are scored as TSG-like. In these tumor types, EGFR mutation patterns are atypical, because EGFR amplification is an early event. This amplification has been linked to increased mutation load in EGFR itself, including in aberrant extrachromosomal copies of EGFR (43). Supplementary Figure S2 indicates the locations of these misclassified regions on the plot.

A disadvantage of working with experimentally derived protein structures is that they are available for a limited number of human proteins (39%). For many of these genes, the structure data are incomplete, so that only a single protein domain or small fragment is represented in PDB. In this work, by careful use of biologic assemblies of PDB structures and also theoretical protein models, we mapped approximately 53% of unique residue positions harboring a TCGA missense mutation into 3D protein space. In a recent study of 21 TCGA tumor types that used a different algorithm and PDB structures only, 11.2% positions were mapped (13). We note that theoretical protein models are well suited for this kind of analysis. HotMAPS considers the center of geometry for each amino acid residue, a metric that is not highly sensitive to atomic-resolution errors common in theoretical protein models (47). The increased sensitivity and coverage of HotMAPS are supported by the number of tumor types in which 3D hotspot regions were detected (19 of 23), the total number of regions detected (107 unique, tumor-type–aggregated gene-level regions and 216 unique tumor-type–specific gene-level regions), and the number of genes in which regions were detected (91). The only previous systematic attempt to find 3D hotspots in TCGA data detected statistically significant regions in 10 genes, based on 21 TCGA tumor types (13), and 9 of these were also detected by HotMAPS.

We hope that some HotMAPS regions found by our algorithm point to novel driver genes; however, functional studies are warranted to find out if they are discoveries or false positives.

An interactive 3D protein viewer where users can submit their own mutations and compare with the HotMAPS catalog (48) is available at http://mupit.icm.jhu.edu/.

HotMAPS software is open source at https://github.com/karchinlab/HotMAPS.

Additional material is available in Supplementary Materials and Methods: detailing mapping from genomic coordinates to protein structures (Supplementary Fig. S3), overview flow chart of HotMAPS (Supplementary Fig. S4), an example of single-residue hotspot region discovery in 1D versus 3D (Supplementary Fig. S5), a stratified analysis of HotMAPS properties by solvent accessibility (Supplementary Fig. S6), and in Supplementary Tables, a list of blacklisted residues (Supplementary Table S9).

No potential conflicts of interest were disclosed.

Conception and design: C. Tokheim, D.L. Masica, R. Karchin

Development of methodology: C. Tokheim, R. Bhattacharya, M.C. Ryan

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): C. Tokheim

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): C. Tokheim, N. Niknafs, R. Kim, M.C. Ryan, D.L. Masica, R. Karchin

Writing, review, and/or revision of the manuscript: C. Tokheim, R. Bhattacharya, R. Karchin

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): D.M. Gygax, M.C. Ryan

Study supervision: R. Karchin

The authors thank Drs. Bert Vogelstein and Jing Zhu of the UCSC Xena team for valuable discussion on the manuscript.

This work was supported by NIH, National Cancer Institute fellowship F31CA200266 (C. Tokheim) and grant U01CA180956 (R. Karchin).

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.