Cancer-driving mutations are enriched in genic regions intolerant to germline variation

Putative somatic variants in gnomAD

Population-level catalogs of human genetic variation allow for the investigation of selective constraint and mutational patterns in the exome. We used the gnomAD database of 125,748 human exomes to survey both germline and somatic variants (1). Although the gnomAD variant calling pipeline was tuned to detect germline variation, we reasoned that we may also be able to identify somatic variants that reach a sufficiently high variant allele frequency to be detected through their germline variant caller. Inherited heterozygous germline variants are expected to have an allelic ratio close to 50%. We observed that the distribution of median allelic balance (AB_median) values for gnomAD variants followed a bimodal distribution, with one distribution centered around 50% and another smaller distribution centered around 20% (Fig. 1A).

Defining OncMTR

We previously introduced a sliding window–based metric, the MTR, which measures purifying selection on missense variation in genic subregions (6). This score demonstrably detects crucial functional domains of proteins that can cause Mendelian disease when mutated in the germ line. Motivated by the overlap between mutations associated with Mendelian disease and cancer, we set out to create a cancer-relevant version of MTR (see Methods) that captures regions that are depleted of germline variation relative to somatic variation. In this study, we defined another variation of the MTR score, namely, MTR^germline. In its construction, MTR^germline is restricted to only those variants achieving an AB_median > 0.3. Taking the well-known cancer gene TP53 as an example, we can observe those genic subregions where the two MTR formulations diverged (Fig. 1B). We then define OncMTR as the difference between these two MTR formulations for each codon and using a 31-codon sliding window (Fig. 1B). Negative scores correspond to regions with the greatest divergence between germline intolerance and somatic variant enrichment. Overlaying OncMTR scores on the AlphaFold-predicted structure of TP53 (13) illustrated that the strongest negative scores correspond to the DNA binding domain, which is the domain enriched for mutations known to drive hematologic malignancies (Fig. 1C).

Using OncMTR to prioritize driver mutations in hematologic malignancies

Motivated by the positive proof of concept demonstrated for TP53, we next tested whether the MTR and MTR^germline distributions differed across other oncogenes and tumor suppressor genes included in the Catalogue of Somatic Mutations in Cancer (COSMIC) Cancer Gene Census (CGC). The CGC is divided into two tiers, with tier 1 containing bona fide cancer genes (n = 556) and tier 2 containing genes that have strong indications of playing a role in cancer but with less expansive evidence than tier 1 (n = 137). The difference between MTR and MTR^germline distributions per gene, calculated via cross-entropy, was significantly higher for tier 1 genes than a random selection of 556 non-CGC genes (P = 5.7 × 10⁻³¹), the remainder of the exome (P = 2.8 × 10⁻⁶⁷), and tier 2 genes (P = 1.1 × 10⁻⁷) (Fig. 2A). The cross-entropy was also significantly larger for tier 2 genes than the remaining genes in the exome (P = 2.6 × 10⁻⁴) (Fig. 2A). We also compared the cross-entropy between genes annotated in IntOGen as loss of function (LoF) versus “activating” (i.e., gain of function) and explored the signals from the CGC tier 1 and 2 gene sets, excluding genes acting only through LoF mechanism (fig. S10). Overall, there is a significant difference between MTR and MTR^germline distributions for all genes regardless of the mechanism of action of the mutational cancer driver genes, when compared against the rest of the exome. However, we observe that restricting analysis to genes with a gain-of-function mechanism (i.e., excluding those with LoF mechanism), there is a stronger enrichment signal compared to the rest of the exome or between activating and “LoF” gene sets. For example, IntOGen activating cross-entropy was more divergent from the rest of the exome than IntOGen LoF and with a significant difference (P = 0.01, Mann-Whitney U test; fig. S10A). The same pattern was observed for the CGC tier 1 and 2 datasets, although it was significant only for CGC tier 1 (P = 0.037 and P = 0.28, respectively, Mann-Whitney U test; fig. S10, B and C). Together, these results support the hypothesis that mitotically mutable genic subregions that are intolerant to germline variation are broadly relevant to cancer.

Distinguishing between cancer-causing driver mutations and neutral passenger mutations remains a central challenge in cancer genomics. We thus tested whether OncMTR could help prioritize somatic mutations that cause hematologic malignancies. We found that the OncMTR scores of a previously defined list of 546 unique leukemogenic driver mutations (table S1) (14) were significantly lower than a size-matched set of random variants (P = 2.97 × 10⁻⁸⁶, Mann-Whitney U test; fig. S1A). A random forest model using OncMTR achieved an area under the receiving operator curve (AUC) of 0.74 in discriminating between these leukemogenic variants and the random set (Fig. 2B). We also calculated transcript-level percentiles for the MTR scores, in which lower percentiles corresponded to lower OncMTR scores. The AUC or the OncMTR transcript percentiles was 0.76, and a combined model that incorporated both the raw OncMTR scores and transcript percentiles achieved an even higher AUC of 0.78 (Fig. 2B).

We next compared the performance of OncMTR with other population- and evolutionary-based scores, including CADD (15), ncER (16), LINSIGHT (17), and phyloP (18), in predicting these leukemogenic variants. Here, we defined the putatively neutral variant set as a random set that was both size- and transcript-matched with the leukemogenic set. Transcript matching ensures that we can detect the power of each score to distinguish driver from neutral mutations within the same genes, rather than with neutral variants in genes that may have different score profiles than the carrier genes. In comparing the mean AUC across fivefold cross-validation via logistic regression, OncMTR outperforms the other scores, achieving an AUC of 0.71, compared to AUCs of 0.61, 0.58, 0.52, and 0.50, achieved by ncER, LINSIGHT, CADD, and phyloP, respectively. This highlights OncMTR’s ability to capture orthogonal and previously unidentified information compared to existing scores in picking up leukemogenic variants, which were not a focus set of variants in the training and development of many existing variant effect scores.

As constructed, OncMTR compares the ratio of observed missense variants to synonymous variants over an expected ratio. Because the expected number of mutations could be influenced by local sequence context, we also constructed a score that accounts for trimer mutability rate in its construction (“OncMTR-mutrate”; see Methods). OncMTR-mutrate was considerably correlated with the original (nonmutability rate–informed) version of OncMTR, achieving a Pearson’s correlation of 0.75 across all transcripts. We also assessed the predictive power of OncMTR-mutrate in classifying leukemogenic variants against a random size- and transcript-matched set of putative neutral variants. OncMTR-mutrate achieved an AUC of 0.69 compared to the original OncMTR that achieved an AUC 0.71 on the same dataset (Fig. 2C). We also make the OncMTR-mutrate version of scores across all transcripts publicly available through the OncMTR Viewer web app (http://oncmtr.public.cgr.astrazeneca.com).

To further assess the capacity of OncMTR to prioritize driver mutations, we trained random forest models with raw OncMTR scores using fivefold cross-validation. The mean AUC for predicting leukemogenic variants was 0.74 (fig. S2). We next compared the performance of OncMTR in distinguishing between a set of random variants and 200 established driver mutations implicated in acute lymphocytic leukemia (ALL), acute myeloid leukemia (AML), chronic lymphocytic leukemia (CLL), diffuse large B cell lymphoma (DLBCL), or multiple myeloma (MM), achieving an AUC of 0.65 (fig. S2) and having significantly disparate OncMTR distributions from each other (P = 4.89 × 10⁻⁵, Mann-Whitney U test; fig. 1B and table S2). Logistic regression–based classifiers achieved similar, albeit marginally lower, AUCs than the random forest models (with AUCs of 0.73 and 0.62 for the two variant sets, respectively), likely due to a small degree of nonlinear distribution of OncMTR scores (fig. S6). To determine which OncMTR cutoff could maximize sensitivity and specificity, we calculated the Youden index from each of the above learning tasks (see Methods). The optimal cutoff was identified as OncMTR < −0.05. Together, these results demonstrate the utility of our population genetics–based approach in identifying genic subregions relevant to hematologic malignancies.

Because the somatic mutations used to calculate OncMTR arose in the blood, we expected that OncMTR would more reliably prioritize driver mutations in hematologic malignancies than in solid tumors. As expected, the OncMTR scores of driver mutations implicated in heme malignancies were significantly lower (P = 2.53 × 10⁻⁹, Mann-Whitney U test; Fig. 2D). More generally, we find that 7.2% of bases in hematological malignancy genes and 24.7% of the driver missense variants among them correspond to low OncMTR regions (<−0.05). This is in comparison to 3.5% of bases in other genes in the genome (P = 2.50 × 10⁻¹⁸² and P = 5.10 × 10⁻⁴¹, respectively, Fisher’s exact test). To determine whether OncMTR performs better for certain subtypes of heme malignancies, we compared OncMTR distributions of putative driver and passenger mutations identified in a recent comprehensive in silico saturation mutagenesis experiment (19). This dataset includes simulated variants across three genes for CLL, nine genes for AML, two genes for non-Hodgkin lymphoma, five genes for lymphoma, six genes for MM, and two genes for ALL (table S10). The OncMTR scores of predicted driver mutations were significantly lower than those of passenger mutations for each cancer subtype, although we observed the strongest separation in CLL (P < 2 × 10⁻³⁰⁸, Wilcoxon test) and AML (P = 1.4 × 10⁻¹⁵⁵, Wilcoxon test) (fig. S3).

We next assessed whether OncMTR can successfully distinguish between ClinVar pathogenic and benign somatic variants. Logistic regression classification between pathogenic and benign or random variants across all protein-coding genes reached an AUC of 0.60 and 0.58, respectively [P = 815 unique pathogenic versus B = 58 unique benign variants; a set R (random) of equal size to P was sampled to compile the random variants; see also Methods] (fig. S4). We next restricted the set of pathogenic somatic variants to those occurring in genes associated with hematologic malignancies and compared to benign or random variants. The AUC was 0.62 in distinguishing between pathogenic and benign variants in hematologic malignancy genes (P = 64 versus B = 20) and 0.67 when comparing to benign variants across the entire exome (P = 64 versus B = 58). The AUCs for pathogenic hematologic malignancy variants versus random variants were 0.61 for random variants restricted to heme genes (P = 64 versus R = 64) and 0.64 for random variants pulled from all protein-coding genes (P = 815 versus R = 815) (fig. S4). These results provide support to this blood-based sequencing version of OncMTR being more powerful in identifying pathogenic mutations implicated with heme malignancies.

Last, to further explore OncMTR’s power to agnostically detect putative oncogenic regions, we scanned all protein-coding genes in ClinVar in search of transcripts that are preferentially enriched for ClinVar pathogenic somatic variants in regions with OncMTR scores at the bottom 20 percentile of the full OncMTR distribution (see Methods). We identified 101 such transcripts from 24 unique genes (P < 0.05, Fisher’s exact test; table S11), with several known cancer driver genes captured, such as TP53, IDH1, ALK, and HNRNPA1 (20). Many of the top-ranked genes are implicated in hematologic malignancies, including MYC, MSH2, and FBXW7 (fig. S5) (21–23).

Genes carrying mutations implicated in both human Mendelian disease and cancer

There have been an increasing number of observations, which germline mutations in certain genic regions cause severe Mendelian disease, whereas identical somatic mutations—occurring later in life and localized to specific tissue(s)—in these regions can cause cancer. Here, we examined OncMTR distributions for three genes implicated in both neurodevelopmental disease and adult leukemias: GNB1, NRAS, and DNMT3A (Fig. 3, A to C, and table S4).

Germline de novo mutations in GNB1 cause a severe developmental syndrome characterized by intellectual disability and other variable features, including hypotonia, seizures, and poor growth (8). Somatic mutations in this gene have been associated with ALL, CLL, and myelodysplastic syndrome (24). Three of four somatic driver mutations in this gene overlap with de novo mutations implicated in developmental delay (p.Asp76Gly, p.Ile80Thr, and p.Ile80Asn) (Fig. 3A) (8). All four mutated residues reside in a low OncMTR region (OncMTR < −0.05) of the gene, which corresponds to the Gβ protein surface that interacts with Gα subunits and downstream effectors (Fig. 3A).

NRAS encodes a RAS protein with intrinsic guanosine triphosphatase activity that has been implicated in multiple hematologic and solid malignancies (25). There are 28 somatic missense variants in this gene at four distinct amino acid positions associated with juvenile myelomonocytic leukemia and AML, and all are residing in low OncMTR regions (Fig. 3B) (14). Two of these mutations have also been reported as causal germline de novo mutations for Noonan syndrome, a developmental delay syndrome that includes congenital heart defects, short stature, and other features (p.Gly13Asp and p.Gly60Glu) (Fig. 3B) (26, 27).

DNMT3A encodes a DNA methyltransferase essential for DNA methylation during human embryogenesis and, when mutated somatically, increases risk of AML (28). In a large study on clonal hematopoiesis of indeterminate potential (CHIP), DNMT3A was found to harbor the largest proportion of CHIP variants of all CHIP-associated genes (29), suggesting that it is highly mitotically mutable. In line with this, the OncMTR distribution of this gene is highly enriched for negative values, even compared to GNB1 and NRAS (Fig. 3C). The R882 amino acid residue of DNMT3A corresponds to a DNA binding residue that is a major somatic mutation hotspot in CHIP and AML (28). De novo germline mutations at this residue are associated with an overgrowth syndrome called Tatton-Brown-Rahman syndrome characterized by tall stature and impaired intellectual development (30). Mutations at the R882 residue are thought to interfere with DNA binding, resulting in functional impairment of the protein and aberrant DNA methylation patterns (31). As expected, we identify that the leukemogenic variants in this gene are enriched in low OncMTR regions (Fig. 3C). Together, these results support the notion that some critically important genic subregions are exceptionally mitotically mutable, and mutations in these regions result in different phenotypic outcomes depending on timing and cellular context (9).

Enrichment of low OncMTR scores in protein domains

One strength of the sliding window approach implemented in OncMTR is that its estimates are independent of biological boundaries, such as annotated protein domains, which are not always well annotated. However, it is known that cancer-causing missense mutations tend to cluster in certain functional domains. We thus tested whether Pfam domains and domain superfamilies were enriched for low OncMTR regions (defined as OncMTR < −0.05). Across human protein-coding genes, low OncMTR regions were significantly enriched for several protein domains previously implicated in cancer, such as homeodomains (Bonferroni adjusted P = 4.9 × 10⁻⁴⁶, Fisher’s exact test), protein kinase domains (Bonferroni adjusted P = 5.25 × 10⁻¹¹⁰, Fisher’s exact test), RING domains (Bonferroni adjusted P = 3.22 × 10⁻⁴⁸, Fisher’s exact test), and others (Fig. 4, A and B, and tables S5 and S6). Furthermore, we found that proteins that had functional domains enriched for low OncMTR scores are significantly enriched in genes with TOPMed leukemogenic variants and known cancer hotspots (Fig. 4C and tables S1 to S3 and S7 to S9) (32). Among these two lists of genes, zinc finger motifs were found to be the most strongly enriched for low OncMTR scores (most significant adjusted P = 2.3 × 10⁻⁵² from the union list, based on Fisher’s exact test; Fig. 4, D to F), in line with their well-established role in cancer development (33). Although the calculation of OncMTR is agnostic to domain annotations, it independently identifies cancer-relevant functional genic subregions.

As low OncMTR scores can indicate functionally important sites within proteins, we sought to understand whether low OncMTR regions residing outside of annotated protein domains were enriched for any functional motifs. We found that low OncMTR regions share characteristics with structurally disordered regions and were not enriched for currently known sequence motifs (fig. S9).

Because of the high enrichment of protein domains within low OncMTR regions, we next explored whether OncMTR captures additional information orthogonal to the presence of a protein domain. To test this, we used a logistic regression classifier based on the presence or absence of a protein domain in predicting the leukemogenic variants. This classifier achieved an AUC of 0.62 (with fivefold cross-validation) as compared to the superior 0.71 achieved by OncMTR on the same dataset.

Informing rare variant collapsing analysis with OncMTR

With increasing adoption of next-generation sequencing to generate case-control cohorts, rare variant collapsing analysis has emerged as a powerful approach to detecting disease-associated genes for both rare and complex disorders. In this approach, the proportion of cases with a qualifying variant (QV) is compared to the proportion of controls with a QV in the same gene. We have previously shown that incorporating an MTR filter in defining QVs markedly improves rare variant collapsing analyses (2). In that phenome-wide association study (PheWAS) on approximately 300,000 exomes in the UKB, the collapsing analyses detected seven genes associated with hematologic malignancies (2). Here, we sought to test whether OncMTR would further improve collapsing analysis signals for hematologic malignancy associations by performing a collapsing analysis on 394,694 European exomes contained in the UKB focused on 1394 chapter IX (neoplasm) phenotypes. We defined a total of eight collapsing models with and without OncMTR filters (table S12). Imposing an OncMTR filter of −0.05 (i.e., only considering missense QVs that fall below this threshold) significantly increased the effect sizes of gene-phenotype associations (P < 0.0001) for each model (Fig. 5A and table S13). We observed genome-wide significant (P < 1 × 10⁻⁸) associations between several heme malignancies and DNMT3A, FBXW7, IDH2, IGLL5, JAK2, SF3B1, SRSF2, TET2, and TP53; in certain cases, the effect sizes were 10-fold greater than without adopting the OncMTR filter (Fig. 5B). We also found that the association between TP53 and CLL only reached significance in models including our OncMTR filter; for example, in the “raredmg” model, this association had a P value of 1.2 × 10⁻⁷ [odds ratio (OR) = 8.8; 95% confidence interval (CI): 4.8 to 16.0], whereas in the “raredmgoncmtr” model, the same association reached a P value of 3.4 × 10⁻¹⁰ (OR = 33.2; 95% CI: 16.1 to 68.7). Thus, applying the OncMTR filter effectively reduces background variation in the setting of gene-level collapsing analysis for hematological malignancy phenotypes, and we advise future large-scale hematological malignancy discovery studies to consider adopting OncMTR filter for improved signal detection.

Reconstructing the MTR with 125,000 samples from gnomAD

We first reconstruct the MTR using a cohort of 125,748 exomes from the gnomAD Consortium (v2, GRCh38 liftover). The formula for deriving the window-based MTR scores has been introduced in the original paper (6)

where the numerator represents the observed proportion of missense variants among the total observed protein-coding variation. The numerator is then scaled by the same proportion computed from the collection of all possible protein-coding variants in the corresponding protein-coding window. A window size of 31 codons has been used for calculating MTR based on the gnomAD cohort, in agreement with the previously published score (6).

The expected proportion of missense variants in a given protein-coding window was calculated by annotating all possible variants of a protein-coding transcript with SnpEff 4.3t using GRCh38.92 as the reference annotation and assuming that all events were equally likely to occur. Annotation with SnpEff focused on single-nucleotide variants (SNVs) that were flagged as PASS variants in the original gnomAD release (v2). Variants annotated as “missense_variant” or “missense_variant&splice_region_variant” by SnpEff represent the set of “missense” variants in the MTR formula. Variants annotated as “synonymous_variant,” “stop_retained_variant,” “splice_region_variant&stop_retained_variant,” or “splice_region_variant&synonymous_variant” by SnpEff were considered as the “synonymous” variants in the same formula.

OncMTR score construction

Using MTR as our basis, we construct the OncMTR score (i.e., Oncology MTR score) to capture protein-coding subregions that are depleted of germline variation relative to somatic variation. We observe that the total distribution of AB_median values across all gnomAD variants (Fig. 1A) is bimodal, with the main peak centered close to 0.5 and a second one emerging for values approximately around 0.2. The AB_median metric represents the allelic ratio between the alleles for each variant, with values close to 0.5, reflecting an equal number of copies being inherited from each parent in heterozygous settings, while truly biological variants that approach zero increasingly reflect variants that more likely arose somatically.

We leverage this observation to construct an alternative version of the original MTR score: excluding any putative somatic variants and using only germline variants from the gnomAD dataset. We achieve that by selecting only variants with AB_median > 0.3, constructing the MTR^germline version of the score. OncMTR is then defined as the difference of the original MTR score from the MTR^germline version

Negative OncMTR values (i.e., MTR^germline < MTR) represent regions that are depleted of germline variation relative to somatic variation, thus allowing to highlight putative oncogenic subregions in protein-coding genes.

Sensitivity analysis on the definition of putative somatic variants

We base the selection of the AB_median cutoff for defining putative somatic variants on the empirical distribution of AB_median values, derived from the gnomAD dataset (Fig. 1A). Specifically, we observe a bimodal distribution, with the two peaks centered at p1: +0.186 and p2: +0.477 on the x axis. The middle point of the two peaks is at m: +0.332 (fig. S8A). To be more conservative with regards to characterizing variants as putative somatic, we initially suggest a cutoff slightly lower than the middle point of the two peaks, picking 0.30 as our starting cutoff value.

Furthermore, we calculated additional versions of OncMTR scores resulting from different values of AB_median thresholds, specifically 0.20, 0.25, 0.28, 0.30, 0.32, and 0.35. Let us define OncMTR_k as the set of OncMTR scores calculated using k as the AB_median threshold across all canonical transcripts. We calculated the Pearson’s correlations for all available {OncMTR_i, OncMTR_j} pairs and observed that OncMTR_0.30 has a correlation of 0.86 and 0.88 with OncMTR_0.28 and OncMTR_0.32, respectively, showing that the score remains very robust around the 0.30 threshold (fig. S8B). Correlation remains considerably high also with OncMTR_0.25 (r = 0.74), while it falls more rapidly when the threshold is shifted to the other direction (r = 0.68 with OncMTR_0.35). That may indicate that increase of AB_median threshold above 0.32 leads to noticeable mixture of putative somatic and nonsomatic variants. On the basis of the above observations, we select 0.30 as the default AB_median cutoff to calculate OncMTR. We also report Pearson’s correlations between all {OncMTR_i, OncMTR_j} pairs for reference (fig. S8B).

Compilation of variant sets

We used a precompiled set of variants known to be drivers of hematologic malignancies in a total of 160 genes (41). This list was generated from recurrent hematologic somatic mutations in the literature and COSMIC, excluding genes with a relatively high proportion of LoF germline mutations. A second, smaller precompiled list focused on genes that were recurrent drivers specifically for myeloid malignancies (14). A third validation set included a list of annotated driver mutations provided through the IntOGen database (42). We restricted this set to “tier 1” (highest confidence) driver mutations observed in hematologic malignancies, which included ALL, AML, CLL, DLBCL, and MM. To identify variants overlapping between hematologic malignancy and neurodevelopmental disease, we identified genes from the leukemogenic set that had entries in the Online Mendelian Inheritance in Man for both adult leukemias and neurodevelopmental disease, which included GNB1, DNMT3A, and NRAS (43).

Classification of oncogenic variant sets with OncMTR

We have performed classification of different oncogenic variant sets (TOPMed leukemogenic and IntOGen drivers) against random variant sets of equal size. We use two supervised models for the binary classification task, logistic regression with “max_iter” = 1000 and a random forest classifier with “max_depth” = 2, to avoid overfitting on the training set. Each classification was performed as a fivefold cross-validation task, and the mean AUC from all folds is reported to reflect the total average performance of each learning task. The implementations of logistic regression and random forest were derived from the sklearn Python package (v0.22.1).

We also estimated the optimal OncMTR cut point for each classification by calculating the Youden’s index from each learning task. The average Youden index from all classification tasks performed with logistic regression was Y_LR = −0.0409 (SD: 0.00126), while for random forest, it was Y_RF = −0.0614 (SD: 0.00057). The mean of the two averages of Youden indexes is −0.05115 or −0.05, after rounding it up to one decimal point for simplicity. We thus consider OncMTR values of less than −0.05 to have the most distinctive power.

Identifying OncMTR regions significantly enriched for ClinVar somatic variants

For this analysis, we use all ClinVar somatic variants (ORIGIN = 2) from the GRCh38 release (last accessed on 9 June 2019), focusing on those annotated as missense or synonymous. We consider as pathogenic variants those annotated as “Pathogenic” or “Likely_pathogenic” and as benign those annotated as “Benign” or “Likely_benign” (based on ClinVar). Classification between pathogenic and benign (or random) variant sets was performed with a logistic regression classifier with fivefold cross validation (sklearn, Python package v0.22.1). When restricting the classification to heme-implicated genes, we derived those gene sets based on the IntOGen annotation (table S10).

To identify genes/transcripts across the exome that are preferentially enriched for ClinVar somatic pathogenic variants in regions with low OncMTR scores, we use Fisher’s exact test. Specifically, we scan across each transcript and identify what percentage of the codons in each transcript achieves an OncMTR score at the bottom 20 percentile of the full OncMTR distribution (across the entire transcript). Then, we check whether known pathogenic or likely pathogenic ClinVar missense variants preferentially land in these codons (i.e., corresponding to low OncMTR scores) compared to the rest of the transcript. We apply Fisher’s exact test to evaluate the enrichment of each set of regions, i.e., those with low OncMTR scores versus the rest of the transcript. Eventually, we rank each transcript based on the OR and significance of the Fisher’s exact test enrichments (table S11).

Enrichment of low OncMTR scores in protein domains

To describe the functional context of OncMTR, we calculated enrichment of constrained regions in protein domain families. Residues within each canonical transcript (as defined by UniProtKB) were divided into two classes on the basis of their corresponding OncMTR scores: less than −0.05 (constrained; as defined by Youden’s index) and greater or equal to −0.05 (relaxed). Domain and clan annotations for the human proteome were taken from the Pfam 34.0 database. DNA binding domains were pulled from a previous compendium (44). The final set of the canonical human proteome consisted of 18,313 annotated proteins. Enrichments of the constrained regions in protein domains were tested with the Fisher’s exact test, followed by Bonferroni correction and with significance level of adjusted P value of 0.05.

Analysis of low OncMTR regions outside of protein domains

Low OncMTR scores can indicate functionally important protein regions. We analyzed whether any genes had low OncMTR regions (i.e., spanning at least 30 residues with OncMTR lower than −0.05) residing outside of protein domains, but we did not identify any unannotated regions that could constitute previously unknown domains. According to UniProt annotations, 38% of them (372 of 978 in total) overlap with predicted disordered regions.

To find common sequence motifs that would be indicative of functionally important sites, we first obtained clusters of low OncMTR regions of similar sequences (using CD-HIT, at 60% sequence identity). We found only three clusters consisting of multiple (i.e., two) short sequences. To find sequence motifs shared within those clusters, we ran MEME searches, which did not return any significant motifs.

Last, to show that those low OncMTR regions do not significantly differ from other regions outside of Pfam domains in terms of structural disorder, we analyzed their per-residue confidence scores [pLDDT (predicted local distance difference test)] derived from AlphaFold2 models (13). In general, regions with lower pLDDT scores have a higher chance of being unstructured. As shown in fig. S9, pLDDTs for Pfam domains are significantly higher than for the regions outside of annotated domains.

UKB collapsing analysis

Collapsing analyses were performed using the 394,694 exomes available in the UKB cohort (45). The UKB is a prospective study of approximately 500,000 participants aged 40 to 69 years old at time of recruitment. Participants were recruited in the United Kingdom between 2006 and 2010 and are continuously followed. The average age at recruitment for sequenced individuals was 56.5 years, and 54% of the sequenced cohort is of female genetic sex. Participant data include health records that are periodically updated by the UKB, self-reported survey information, linkage to death and cancer registries, collection of urine and blood biomarkers, imaging data, accelerometer data, and various other phenotypic end points. All study participants provided informed consent, and the UKB has approval from the North-West Multi-centre Research Ethics Committee (11/NW/0382).

We performed a gene-based collapsing analysis on 1394 chapter IX (neoplasm) phenotypes adopting our previously described approach (2). We implemented a total of eight dominant collapsing models with and without OncMTR filters (table S12). Using SnpEff (46), we defined protein truncating variants (PTVs) as variants annotated as exon_loss_variant, frameshift_variant, start_lost, stop_gained, stop_lost, splice_acceptor_variant, splice_donor_variant, gene_fusion, bidirectional_gene_fusion, rare_amino_acid_variant, and transcript_ablation. We defined missense as follows: missense_variant_splice_region_variant and missense_variant. Nonsynonymous variants included the following: exon_loss_variant, frameshift_variant, start_lost, stop_gained, stop_lost, splice_acceptor_variant, splice_donor_variant, gene_fusion, bidirectional_gene_fusion, rare_amino_acid_variant, transcript_ablation, conservative_inframe_deletion, conservative_inframe_insertion, disruptive_inframe_insertion, disruptive_inframe_deletion, missense_variant_splice_region_variant, missense_variant, and protein_altering_variant. We derived allele frequencies from gnomAD (1). The raredmg, raredmg_OncMTR, flexdmg, and flexdmg_oncMTR models incorporated a REVEL cutoff of REVEL ≥ 0.5 to restrict to putatively damaging missense variants (47).

To compute P values, the carriers of at least one QV in a gene were compared to the noncarriers. The difference in the proportion of cases and controls carrying QVs in a gene was tested using a Fisher’s exact two-sided test. Variants were required to pass the following quality control criteria: minimum coverage of 10×; annotation in consensus coding sequence (CCDS) transcripts (release 22; approximately 34 Mb); at most 80% alternate reads in homozygous genotypes; percent of alternate reads in heterozygous variants ≤0.25 and ≥0.8; binomial test of alternate allele proportion departure from 50% in heterozygous state P > 1 × 10⁻⁶; genotype quality score (GQ) ≥ 20; Fisher’s strand bias score (FS) ≤ 200 (indels) ≤ 60 (SNVs); mapping quality score (MQ) ≥ 40; quality score (QUAL) ≥ 30; read position rank sum score (RPRS) ≥ −2; mapping quality rank sum score (MQRS) ≥ −8; DRAGEN variant status = PASS; the variant site achieved 10-fold coverage in ≤25% of gnomAD exomes; and if the variant was observed in gnomAD exomes, then the variant achieved exome z score ≤ −2.0 and exome MQ ≤ 30. We excluded 46 genes that we previously found associated with batch effects (2).

For all models, we applied the following quality control filters: minimum coverage of 10×; annotation in CCDS transcripts (release 22; approximately 34 Mb); at most 80% alternate reads in homozygous genotypes; percent of alternate reads in heterozygous variants ≤0.25 and ≥0.8; binomial test of alternate allele proportion departure from 50% in heterozygous state P < 1 × 10⁻⁶; GQ ≤ 20; FS ≥ 200 (indels) ≥ 60 (SNVs); MQ ≤ 40; QUAL ≤30; read position rank sum score ≤ −2; MQRS ≤ −8; DRAGEN variant status = PASS; the variant site achieved 10-fold coverage in ≤25% of gnomAD exomes; and if the variant was observed in gnomAD exomes, then the variant achieved exome z score ≤ −2.0 and exome MQ ≤ 30. We excluded 46 genes that we previously found associated with batch effects 10.

We defined the study-wide significance threshold as P < 1 × 10⁻⁸. We have previously shown using an n-of-1 permutation approach and the empirical null synonymous model that this threshold corresponds to a false positive rate of 9 and 2, respectively, of ∼346.5 million tests for binary traits in the setting of a collapsing analysis PheWAS (2).

This PDF file includes:

Other Supplementary Material for this manuscript includes the following: