Interaction Landscape of Inherited Polymorphisms with Somatic Events in Cancer

Cancer is a complex genetic disease influenced by both inherited variants in germline DNA and somatic alterations acquired during formation of the tumor (1, 2). Recently, large multicenter efforts such as The Cancer Genome Atlas (TCGA; ref. 3) and the International Cancer Genome Consortium (ICGC; ref. 4) have performed a series of detailed analyses of the somatic alterations affecting tumor genomes. These studies have focused primarily on identifying recurrent somatic alterations, uncovering immense heterogeneity within and between tumors (2).

Prior to tumor genome sequencing, many genes that play a role in cancer were discovered through studies of the germline (5). Linkage studies in families with inherited, typically childhood, cancers identified rare germline mutations in genes related to DNA damage repair, RAS signaling, or PIK3 signaling (2, 6). In contrast to childhood cancers, adult tumors have largely been considered “sporadic”; however, mounting evidence points to a potentially substantial influence from the germline (7). Studies of homogenous populations in Scandinavia identified a genetic component underlying common adulthood cancers (8, 9); more recently, many loci have been implicated through genome-wide association studies (GWAS) of various types of cancer (10–13). Large-scale systematic sequencing of over 4,000 tumors (12 cancer types) from TCGA found rare germline truncations in 114 cancer susceptibility–associated genes, ranging in frequency from 4% (acute myeloid leukemia) to 19% (ovarian cancer), including BRCA1, BRCA2, FANCM, and MSH6, which are associated with increased somatic mutation frequencies (14).

Other evidence has emerged that germline variants and somatic events can be intricately linked. For example, a recent analysis of rare germline variants and somatic mutations in ovarian cancer highlighted novel ovarian cancer genes and pathways (15). Recent reports have also associated specific haplotypes with JAK2^V617F mutations in myoproliferative neoplasms (16) and with EGFR exon 19 microdeletions in non–small cell lung cancer (17). Germline variation was also found to influence gene expression in tumors (18, 19). Together, the variants identified thus far are estimated to explain at most 20% of the likely germline contribution to cancer (11), suggesting the existence of many as-yet-uncharacterized genetic determinants.

Here, we integrate germline genotypes with somatic changes from TCGA to obtain a pan-cancer view of how common inherited variation can prime the later progression of tumors. We seek to identify genetic associations that explain two major classes of somatic events: (i) the tissue site where the tumor develops, and (ii) which specific cancer genes are mutated. This analysis identifies an array of germline polymorphisms that increase or decrease the risk for these somatic events, some of which were previously known, but most of which are new discoveries.

We obtained common germline variants and somatic tumor mutations for 6,908 patients from the TCGA Research Network (20). We restricted our analysis to germline variants present at minor allele frequencies (MAF) of 1% or higher in this patient group, resulting in 706,538 autosomal single-nucleotide polymorphisms (SNP) organized into ∼1.6 million haplotypes (Fig. 1A). Both SNPs and haplotypes were studied, as the two marker types have complementary power to detect common and rare disease-associated variants, respectively (21–23). Examination of patient SNP profiles showed evidence of population substructure (24), even among the majority of patients who self-identified as European ancestry (Supplementary Fig. S1A). Of these patients, we retained 4,165 who also clustered tightly with Europeans from the HapMap III reference population (ref. 25; discovery cohort, Fig. 1B). For later replication of our findings, we used a validation cohort consisting of 1,789 additional TCGA patients for whom full data (tissue type and somatic alterations) became available only after the start of our study (Supplementary Fig. S1B and S1C).

We first sought to identify germline markers associated with the site at which the tumor develops, i.e., incidence of tumors of a particular tissue type. Previous studies have focused on individual tumor types in isolation; thus, it remains unclear to what extent cancer-associated alleles are general versus site specific. All tumors of a particular type were compared with all other types pooled; separate comparisons were performed for each of the 22 TCGA tumor types. SNP or haplotype markers were selected as initial candidates if their P values of association were less than the “suggestive” threshold commonly used in GWAS (26): 1 × 10⁻⁵ by Fisher exact test, corrected for the number of tumor types assessed. Different tumor types were collected at different times by TCGA; to minimize false-positive associations due to such batch effects, we removed markers displaying strong association with batch or plate (Supplementary Fig. S2) and used genomic control to adjust for inflated P values (Supplementary Fig. S3). In cases in which multiple tumor types were each weakly associated with a genotype, these were combined for testing for stronger association as a group. All candidate associations were then tested in the validation cohort, and empirical false discovery rates (FDR) were estimated using matched numbers of random associations. By this procedure, we identified 916 markers of potential interest, 395 of which could be replicated at an empirical FDR < 0.25 (Fig. 1C and 2A; Supplementary Table S1). Additional markers at each locus were imputed to provide a more complete view of the association. Although all TCGA patients have been diagnosed with some type of cancer, these germline markers provide predictive information about where the tumor develops.

As an example, markers at locus 8q24.13 were identified for their specific and significant association with breast cancer (P < 2 × 10⁻¹⁰; Fig. 2B), which was reinforced by a second finding of association with age of breast cancer onset (Fig. 2C). Although this region had not previously been implicated in cancer predisposition, it had been reported to harbor frequent somatic amplifications in breast cancer cell lines (27). In addition, breast tumors overexpressing genes in this region are associated with shorter time to recurrence (28–31). Among the many other loci associated with tumor type, we identified an SNP at 13q14.2 that falls into an intron of the tumor suppressor RB1 and is associated with multiple tumor types, and an SNP at 11q22.3, a region encoding the cancer genes DDX10 and ATM that is associated with breast cancer (Supplementary Table S1). Epigenetic silencing of DDX10 was recently reported in ovarian cancer (32), and rare germline coding variants in ATM have previously been reported in familial breast cancer (33).

To examine the correspondence of this TCGA analysis with previous GWAS of cancer, we analyzed 557 cancer-associated SNPs recorded in the National Human Genome Research Institute (NHGRI) GWAS Catalog (26). Most of these previous SNPs had been identified based on an increase in risk for developing cancer as compared with nondiseased controls, whereas the SNPs identified in our study are based on an increased prevalence of tumors of one particular tissue as compared with other tissues. Despite this difference, the previously published markers of cancer risk had substantial correspondence to the markers associated with tumor site, with the strongest signal seen for approximately 15 markers (Fig. 2D; Supplementary Fig. S4A–S4D; Supplementary Table S2). The loci validated by both types of study can be prioritized as a set of clearly reproducible cancer risk factors that are also specific to tissue of origin.

Even where the new associations recapitulate a previously reported cancer locus, new insights may be gained, as illustrated for the thyroid cancer locus 9q22.23 (refs. 34, 35; Fig. 2E). SNP rs1867277 at this locus, in the promoter region of the FOXE1 transcription factor, was previously reported to be the cause of association; the minor allele was shown to recruit the USF1/USF2 transcription factors, resulting in allele- specific FOXE1 transcription (35). In TCGA, we observed allele-dependent expression of FOXE1, but also of three nearby genes (Fig. 2F): c9orf156 (NAP1), involved in nucleosome assembly with the potential to widely alter gene regulation (36); TRIM14, which may regulate the expression of multiple cancer genes (37); and CORO2A, which derepresses expression of Toll-like receptors to initiate the inflammatory response (38). In TCGA, FOXE1 expression decreased with an increasing number of minor alleles at this locus, whereas Landa and colleagues reported increased FOXE1 expression (ref. 35; Fig. 2F). This discrepancy in FOXE1 expression levels may arise from performing the analysis in patient tumors, whereas the previous report was based on experiments in cancer cell lines, and in part because risk in this region is driven by at least two distinct haplotypes that may affect the expression of different genes (39, 40). Interestingly, targeted overexpression of FOXE1 was not in itself capable of causing thyroid cancer in mice (41), suggesting other genes in this region may play a role. Thus, reanalysis of this locus using genotype and mRNA expression from the TCGA suggests that its role in thyroid cancer may be more complex than previously appreciated.

We next sought to identify associations between inherited germline variants and the occurrence of somatic mutations in particular cancer genes. We hypothesized that germline background could generate a context in which a loss- or gain-of-function event in a particular gene could be advantageous to a tumor. As loss- or gain-of-function is commonly achieved through DNA mutation or amplification/deletion, we considered both somatic mutation and somatic copy-number changes in our analysis. We used a previously published set of 138 cancer genes that had been compiled based on a strong signal of positive selection for subtle somatic mutations (SSM, i.e., point mutations, insertions, or deletions) or copy-number variants (CNV) in tumors (2). Such evidence of positive selection suggests that the alterations affecting these genes are likely to be enriched for functional drivers. Indeed, genes under positive selection in cancer exhibit a bias toward mutations predicted to interfere with protein activity (42, 43). Because these genes play a central role in many cancers, we expected that germline loci associated with their alteration status in tumors would provide insight into specific biological contexts that influence which cancer genes most effectively promote tumor growth and survival. For the most frequently mutated genes in this set, such as TP53, which is altered in approximately 35% of TCGA tumors, we estimated that we had 80% statistical power to detect changes in mutation rate of ≥1.8-fold (Supplementary Fig. S5A). For genes near the median mutation frequency in this set, such as BRCA1, which is altered in ∼3% of tumors, we were powered to detect changes in mutation rate of ≥3 fold.

Association testing identified a total of 62 associations between germline markers and specific cancer genes, of which 17 were found to replicate in the validation cohort at an empirical FDR < 0.25, and of which 35 were deemed of potential interest (empirical FDR < 0.5; Fig. 1C; Supplementary Table S3). These 35 associations covered 28 germline loci and 20 cancer genes. Somatic alteration rates for cancer genes tended to increase with the number of minor alleles at the germline locus, with the largest increases coinciding with the least frequent germline alleles (Fig. 3A–C). The effect sizes were quite large, corresponding to increases in mutation frequency of 1.8- to 14.8-fold (Fig. 3A), and were well correlated between the discovery and validation cohorts (r = 0.43, P < 0.01; Supplementary Fig. S5B). Such large effects linked to common SNPs are particularly striking in comparison with typical GWAS, where effect sizes tend to fall beneath 2-fold (13). For example, a haplotype on 15q22.2 resulted in a 13-fold increased chance of acquiring a CNV affecting GNAQ.

In what follows, we further investigate validated germline–cancer gene associations at 16p13.3 and 19p13.3 as instructive examples. At each of these loci, one or more genes encoded at the germline locus participate in the same biological pathway as the somatically mutated cancer gene (Fig. 3A and C, associations labeled in red).

In one noteworthy connection between germline and somatic gene mutations, we observed a haplotype at 16p13.3 that was associated with a markedly increased risk of somatic mutation in SF3B1, encoding a component of the U2snRNP spliceosome (Fig. 4A, a striking >8-fold increase in mutation rate under the homozygous minor allele). The most strongly associated haplotype encompassed an enhancer in the fourth intron of RBFOX1 (Fig. 4B), the only gene at this locus. Like SF3B1, RBFOX1 also encodes an RNA-binding protein involved in splicing, which influences the inclusion or exclusion of exons in alternatively spliced isoforms (44). This suggested a model whereby the germline configuration of RNA splicing, linked to variation in RBFOX1 expression, modulates the sensitivity of RNA splicing to subsequent somatic mutations in the U2 spliceosome (Fig. 4C).

In support of this hypothesis, we first found that the minor allele at this locus tends to substantially increase RBFOX1 expression, especially in homozygotes (Fig. 4D, normalized across tissue types). Next, we investigated how the germline state of RBFOX1, somatic mutation of SF3B1, or the combination of both of these factors affected RNA splicing patterns in TCGA patients. In patients with the major RBFOX1 allele, SF3B1 mutation had minimal effect on splicing. In contrast, in patients with the minor allele, SF3B1 mutation led to changes in splicing patterns within a number of transcripts (Fig. 4E). Altogether, we identified splice junctions in 12 genes for which there were significant increases or decreases in the fraction of cryptic transcripts, dependent on SF3B1 somatic mutation (Fig. 4E). Such genes include those driving cellular proliferation and metabolism (ribosomal subunits, protein turnover) as well as YBX1, a regulator of alternative splicing. A functional relationship between RBFOX1 germline status and SF3B1 somatic mutation was also supported by statistical modeling, which revealed a significant interaction of these two factors in predicting the overall fraction of cryptic transcripts in a patient (P < 0.04).

A germline haplotype at locus 19p13.3 was associated with a substantial increase in somatic mutation rate of the PTEN tumor suppressor gene, from approximately 5% for the homozygous major allele to 22% for heterozygotes (Fig. 5A). We noted that two genes at this locus, GNA11 and STK11, function in the PIK3CA/mTOR signaling pathway in which PTEN plays a major repressive role (Fig. 5B). In particular, GNA11 can act as an oncoprotein by activating mTOR signaling (45), whereas STK11 inhibits mTOR activity downstream of PTEN (46, 47). The convergence of these three proteins on the mTOR pathway suggested a model in which the minor allele of 19p13.3 affects mTOR signaling, conferring sensitivity of this pathway to later somatic mutation of PTEN.

In support of this hypothesis, we observed that mRNA expression of GNA11 was higher in the presence of the minor 19p13.3 allele in a majority of tumor types examined (P < 0.05, one-sided Mann–Whitney test with multiple testing correction; Fig. 5C). To investigate the impact of changing GNA11 expression, we next placed transcription of GNA11 under exogenous control in HEK293T cells and measured the relationship between GNA11 expression level and mTOR signaling activity. Increases in GNA11 mRNA led to corresponding increases in mTOR signaling in wild-type cells; however, this effect was greatly magnified by PTEN knockdown (Fig. 5D). Interestingly, we observed that STK11 loss of function (mutation or deletion) was significantly more likely in the presence of a GNA11 gain-of-function event (mutation or amplification; OR: 2.87, P < 1 × 10⁻¹², Fisher exact test). Based on this result, we anticipated that STK11 might also interfere with GNA11 activation of mTOR signaling, because like PTEN it serves as a repressor of mTOR signaling. Indeed, we observed that the increase in mTOR signaling due to GNA11 was greatly magnified under STK11 knockout (Fig. 5E), just as we had observed for PTEN knockdown. Although further work is needed, these findings lend support to a model by which alterations at locus 19p13.3 increase the activity of GNA11, which, in turn, increases the selective advantage provided by PTEN-inactivating mutations during tumorigenesis.

Given this particular sensitivity of GNA11-driven mTOR signaling to STK11 knockout, we hypothesized that GNA11 may also act upstream of STK11. Indeed, GNA11 expression was sufficient to increase the level of STK11 protein phosphorylation with a concomitant increase in phosphorylation of AMPK (Fig. 5F), a direct target of STK11 (Fig. 5B; ref. 48). These findings suggest that GNA11 increases mTOR function while indirectly stimulating AMPK-based inhibition of mTOR via STK11 (Fig. 5B, dotted line).

Because we had identified 28 germline loci that increase the occurrence of somatic alterations in well-known cancer genes, we hypothesized that these same loci might influence mutation rates of other genes not previously linked to cancer. To explore this idea, we used MutSigCV (49) to identify all genes mutated more frequently than expected when grouping TCGA patients in the discovery cohort not by tissue, as has been performed successfully many times in the past (49), but according to the state of their germline at each of the 28 loci. This analysis identified 20 additional genes that had a significantly higher somatic mutation rate than expected when analyzed in a specific germline context (Fig. 6). Some genes had an elevated mutation rate in the presence of the minor allele (Fig. 6, purple blocks), whereas others were mutated only with the homozygous major allele (Fig. 6, green blocks). An elevated mutation rate relative to the background expectation is a signature of positive selection; thus, this approach has the potential to identify genes under positive selection in cancer on a particular genetic background.

Of the 20 genes identified by MutSigCV, 15 had not previously been identified as frequently mutated in any TCGA or ICGC cancer genome study. Evidence suggests that several of these genes could be relevant to carcinogenesis; for instance, CD86 plays a role in the early T-cell response (50), TPTE encodes a phosphatase with very high sequence similarity to the known tumor suppressor PTEN (51), and DEFB115 is often comutated with the protein kinase C isozymes PRKCG, PRKCH, and PRKCQ, which have been implicated as tumor suppressors (52). Thus, mutation analysis based on cohorts defined by genotype, rather than tissue, can provide a powerful strategy to identify novel genes in cancer.

Thus far, most studies of the cancer genome have been concerned with understanding the somatic mutations or transcriptional changes that arise during tumor progression. Recently, however, there has been a growing focus on the role of inherited variation in adult cancer, with a view toward next-generation risk assessment and prevention (53). Here, we have described a genome-wide analysis of germline–somatic interactions, based on the availability of germline genotypes and somatic phenotypes for most TCGA patients. We found evidence that genetic background can influence the somatic evolution of a tumor in at least two ways: in determining the site of tumorigenesis and by modifying the likelihood of acquiring mutations in specific cancer genes. The ability to analyze matched genotypes, somatic genome alterations, and mRNA expression profiles enabled us not only to identify germline–somatic genome linkages but also to investigate which genes at a given locus are affected transcriptionally and thus may mediate the effects of germline variants. Moreover, by grouping tumor samples according to the state of inherited genetic variants, it was possible to discover recurrent somatic mutations that were not identifiable by any previous cancer analysis, highlighting novel cancer gene candidates.

Collectively, this resource of germline–somatic interactions in cancer will generate many testable hypotheses about the molecular mechanisms underlying adulthood cancer risk. Interactions that are most readily interpreted are those for which genes at the germline locus function in the same biological pathway as the gene that is somatically altered. For example, a locus associated with PTEN mutation encoded two cancer genes in the same pathway, GNA11 and STK11; all three of these genes regulate growth signaling via mTOR (Fig. 5). We were able to show experimentally that both STK11 and PTEN interfere with GNA11 activation of mTOR, suggesting that germline variants increasing the activity of GNA11 could increase the selective advantage of PTEN mutations during tumor progression. Another locus associated with SF3B1 mutation encoded RBFOX1; both genes regulate alternative RNA splicing (Fig. 4). Further investigations showed that SF3B1 mutation is associated with significant differences in cryptic RNA splicing, exclusively in individuals with the minor allele. We expect that further analysis of the many additional germline–somatic interactions reported in this resource will provide clues about the underlying molecular relationships that promote cancer.

Some frequently altered cancer genes, such as TP53 and PIK3CA, were not found to be influenced by common germline variants. This result is somewhat puzzling because, due to their frequent mutations in TCGA, our analysis was highly powered to find germline associations with these genes (Supplementary Fig. S5A). One explanation is that the mutation of these genes is so critical to cancer that it provides a selective advantage to tumor cells regardless of germline background. Another possibility is that germline interactions take place with specific mutation sites within the gene, but not with the gene as a whole, as was previously reported for germline control of JAK2 mutation (16).

An important question is to what extent germline variants, including the ones identified here, can affect precision medicine. Thus far, translating cancer GWAS to the clinic has been challenged by the generally small effect sizes of risk variants identified (13). Here, among individuals with cancer, we identified germline alleles that had very large effects on the progression of later somatic events (e.g., a 15q22.2 allele that increases somatic alteration of GNAQ by >10-fold, Fig. 3A). Moreover, a number of the reported associations were with minor germline alleles that are quite common in the human population (e.g., markers at 10p15.1 or 18q21.2 at frequencies >30%, Fig. 3B), or involve cancer genes with very high somatic mutation frequencies, such as PTEN, which is mutated in >60% of some cancer types (54). Such interactions may be of particular interest for their ability to stratify patients with cancer, although the true utility of these large effect sizes in a multigenic disease setting remains to be determined. Future studies of germline associations with other clinically important somatic phenotypes, such as site of metastasis, measures of aggressiveness, or therapeutic response, may also reveal large networks of informative interactions. Ultimately, the influence of many germline variants on specific somatic changes in tumors suggests it might be possible to anticipate key events during tumor development, enabling a preventative rather than reactionary approach to therapy.

Data were obtained from the TCGA and divided into a discovery and a validation cohort based on availability of all data types before or after May 2014. Distribution of tumor types in each cohort is shown in Supplementary Fig. S1C. Assignment of each sample to discovery or validation cohorts is provided in Supplementary Table S4. We acquired genotype, clinical, copy-number, and somatic mutation data for all available samples. Human genome assembly GRCh37/UCSC hg19 coordinates were used for all genomic data.

Normal (nontumor) level 2 genotype calls generated from Affymetrix SNP6.0 array intensities using the BirdSuite software (55) were retrieved from the TCGA data matrix. In these files, each SNP was annotated with an allele count (0 = AA, 1 = AB, 2 = BB, −1 = missing) and a confidence score between 0 and 1. Genotypes with a score larger than 0.1 (corresponding to an error rate of >10%) were set to missing and the data were reformatted with PLINK v1.9 (56).

Information on tumor type was available for all samples. Exome-wide profiles of somatic DNA mutations and genome-wide CNVs were accessed from the Broad Firehose Analysis Pipeline (ref. 57; April 16, 2014, release; Supplementary Fig. S1B). Somatic mutations for 1,099 validation samples were obtained using MuTect (58) and Somatic Indel Detector from GATK release 2.2-2 (59).

In addition to TCGA CNV calls, we used the PennCNV (v1.0.3) Affymetrix pipeline to call CNVs using default parameters (60). The Affymetrix Power Tools software package was used to generate signal intensity data from raw CEL files. PennCNV was used to split the signal intensity files by individual, generate CNV calls, merge adjacent CNVs, and annotate CNV calls with UCSC hg19 knownGene annotation.

For each patient, a gene was considered mutated if a DNA mutation mapped to a protein sequence change with SNPEffect4.0 (ref. 61; GRCh37.31, canonical transcripts only). A gene was considered CNV altered if it received a GISTIC2.0 (62) thresholded score of −2 or 2 (0 or 4+ with PennAffyCV) corresponding to homozygous deletion or high level amplification similar to previous studies (63). SSMs, CNVs, and their union were each recoded as a binary vector for every tumor, denoting the presence or absence of an event for each of 138 cancer genes.

CNV and somatic mutation rates were used as covariates for statistical modeling. Somatic mutation rate was modeled by the nonsynonymous mutations per MB reported in the Firehose MutSig analysis (49), and CNV rate was modeled as the number of CNVs per sample, approximated by the number of times each sample appeared in the focal_input.seg.txt file from the Firehose GISTIC2.0 analysis (62). In each case, values could be assigned to most tumor samples (Supplementary Fig. S1B).

We discarded 322 SNPs with probe names that did not match the hg19 UCSC genome browser Affymetrix track (track: SNP/CNV Arrays, table:snpArrayAffy6). Allele counts were converted to alleles using the definitions in metadata distributed with Birdsuite, and negative strand genotypes were flipped to the positive strand using PLINK.

European ancestry samples were identified using the TCGA metadata (Supplementary Fig. S1A). Genotypes were filtered with PLINK to remove SNPs with call rate < 95%, SNPs with MAF < 1% and individuals with genotype coverage < 95%. Additional samples were dropped due to ambiguous or conflicting gender assignment and unexpectedly high or low rates of heterozygosity. SNPs not in Hardy–Weinberg equilibrium (P < 10⁻⁵) and nonautosomal SNPs were discarded. Batch effects associated with processing groups of samples together (plate effects) can lead to bias in estimates of allele frequencies between groups of samples (64). Thus, we discarded 10,436 additional SNPs demonstrating strong plate associations (P < 10⁻⁸). SNPs associated with > 4 plates were retained. After filtering, 706,538 out of 906,600 SNPs remained. Post association testing, any tumor type–associated markers with plate or batch associations (P < 10⁻⁴) were excluded from further analysis (Supplementary Fig. S2).

Affymetrix SNP6.0 genotypes for 1,789 validation samples were filtered for call rate, coverage, and MAF. SNPs not present in the discovery set were removed, and SNPs missing from the discovery set were imputed using PLINK. Validation set genotypes were phased using Beagle v3.3.2, and haplotype markers were assigned based on agreement of phased SNP sequences with best associated SNP sequence for each discovery set haplotype as determined by the cluster2hap utility and as previously described (65).

To further control for population substructure, we performed principal component analysis with the combined TCGA and HapMap Phase III populations. We discarded 107 samples that did not cluster closely with HapMap III European populations. HapMap Phase III genotypes (25) were obtained from the NCBI HapMap ftp site and lifted to hg19 using the liftOver utility (66). Genotypes were merged and reduced to a set of independent SNPs by linkage-based filtering using PLINK. The reduced set of 33,724 independent SNPs was used to calculate pairwise identity-by-state (IBS) between all individuals. After performing PCA on the IBS matrix, we removed 10 TCGA individuals clustering with Masai and Yoruban samples and 97 individuals clustering more closely with individuals from the Mexican and Gujarati populations.

Haplotypes were inferred using Beagle software v3.3.2 (67) and encoded as binary markers using the psuedomarkers.jar java utility (parameters: edgecount = 60 and othercount = 10), resulting in 1,598,830 haplotype markers.

Power to detect association with gene alteration status was estimated using the Genetic Power Calculator (68) via the ldDesign R package (69). Empirical estimates of case–control ratio were modeled by the number of TCGA samples harboring alterations in each gene. Allele frequency and prevalence were set to 0.01. D′ was set to 1 and the same frequency was used for the marker and quantitative trait locus, simulating the best-case scenario where the causal SNP is genotyped.

When testing for marker association with tumor type, samples were partitioned using a one-versus-rest strategy. For marker association with somatic alteration status at each of the 138 cancer genes, samples were partitioned into two groups based on the presence or absence of a somatic alteration (Fig. 1A).

Association testing between all marker–gene and marker–tumor type pairs was performed using a two-sided Fisher exact test with PLINK. This test is not affected by class imbalance, which occurs for most phenotypes in this study. Candidate associations were selected based on the “suggestive” GWAS significance threshold, adjusted by the number of phenotypes tested (1 × 10⁻⁵/22 for primary tumor type and 1 × 10⁻⁵/138 for mutation status of cancer genes). P values were adjusted for inflation using genomic control (Supplementary Figs. S3 and S6). Because allelic tests are biased for SNPs violating Hardy–Weinberg equilibrium, candidate associations were subjected to permutation to obtain empirical P values; somatic alterations were permuted across samples 1 × 10⁸ times to generate an empirical null distribution capable of providing a P value with a resolution of 1 × 10⁻⁸.

To gain power to identify SNPs associated with multiple tumor types, tumors with correlated ORs were tested for association as a group. Specifically, loci receiving at least one weak but insignificant association (P < 0.05) were reevaluated by grouping correlated tumor types. All tumors at the locus with an OR > 1.2 were grouped and a Fisher exact test was performed to determine their combined OR and P value. This joint P value was compared with an empirical distribution obtained by permuting tumor type label 1,000 times, grouping tumor types with OR > 1.2 and repeating the Fisher exact test. Markers were selected for validation if they were among the top 5 most significant tests after permutation and had a combined OR > 1 and a P value < 1 × 10⁻⁵.

Firth's penalized logistic regression (R logistf package; ref. 70) was used to control for covariates, including population substructure (first two principal components PC1 and PC2), individual-specific somatic mutation and/or copy-number alteration rates, and gender. For associations with cancer genes, we also controlled for primary tumor type. The Wald test was used to reject the hypothesis that the slope contribution of the genotype in the fitted model was 0 (i.e., that the genotype does not provide information about the phenotype) in the presence of covariates.

For both gene alteration status and tumor type associations, validation testing was performed using a one-sided Fisher exact test consistent with the odds ratio observed in the discovery screen. Candidate associations achieving an empirical FDR < 0.5 were considered to be of potential interest and are listed in Supplementary Tables S1 and S3. We considered a marker to have “validated” at an empirical FDR < 0.25. These associations are visualized as a custom track in the UCSC Genome Browser (http://ideker.ucsd.edu/apps/germline/).

For each candidate marker–phenotype pair carried forward for validation testing, phenotype labels were permuted 10,000 times. Each of these 10,000 sets of permutated marker–phenotype pairs was then evaluated in the validation cohort. The expected number of false discoveries at a particular significance threshold was estimated as the mean number of marker–phenotype pairs detected at that threshold across the 10,000 permutations. The empirical FDR was estimated as the ratio of false discoveries to the total number of candidates detected at a particular significance threshold. Empirical FDRs were estimated at significance threshold intervals of 0.05.

MutSigCV version 1.4 (49) with default parameters and configuration files was used to identify genes with an elevated mutation rate on each germline background. Discovery-phase samples were divided into two groups: those with one or more copies of the minor allele, and those with none. MutSigCV was run separately on both groups to identify genes that were mutated at higher frequency than expected according to their gene-specific background mutation model.

Expression data (level 3 normalized data) for all TCGA samples were obtained from the April 16, 2015, Firehose run (71). To determine a pan-cancer relationship between gene expression and genotype, log₂ expression was regressed on minor allele count and tumor type. Significance of effect size was determined by first regressing expression on primary tumor type alone. Variance due to tumor type was removed by summing the expected value over all samples and the residuals after regressing on tumor type. Effect sizes were estimated using a one-sided t test to evaluate the difference in distributions for samples with one or more copies of the minor allele versus samples homozygous for the major allele. T tests were performed when at least 5 samples were present in both groups.

TCGA RNA-sequencing data for 95 samples with SF3B1 mutations and 105 controls (Supplementary Table S5) were processed as previously described (72). The percentage of cryptic reads at each splice site across the genome was calculated for each patient. Only splice sites in the top 85th quartile of coverage were considered for the analyses. Percentages were log transformed and averaged across patients with matching genotypes in order to determine the variation in splicing due to the minor allele and SF3B1 mutation status. Specifically, the number of sites that were differentially spliced between WT and SF3B1 mutant tumors was determined for individuals with one or more copies of the minor allele, and compared with the same statistic for individuals homozygous for the major allele. A site with 2-fold change in cryptic splicing was considered differentially spliced. In order to control for the effect of sample size differences when partitioning tumors on germline genotype on identifying differentially spliced junctions, samples with the major allele were downsampled 1,000 times. A two-sample t test was used to determine whether cryptic splicing occurred at the same frequency for samples with the minor allele and without an SF3B1 mutation as compared with all other samples. To correct for any biases, the observed genotypes were randomly reassigned 10,000 times and the two-sample t test reperformed.

A gene was considered to be in cis with a germline marker if it was encoded within 1 MB of the marker. Coordinates of promoters, exons, and introns were determined from Gencode V19 basic obtained from the UCSC Genome Browser (73). Promoters were assumed to fall within 2 KB 5′ of a gene's transcription start site. Enhancer coordinates were obtained from Fantom 5 (ref. 74; robust enhancer set).

Post association testing, markers were assumed to occur at the same locus if pairwise linkage disequilibrium (LD) was sufficiently large (r² > 0.2; ref. 75). To determine overlap of new markers with existing cancer GWAS SNPs, GWAS SNPs within 1 MB of candidate loci were pulled from the NHGRI GWAS Catalog (downloaded April 2015; ref. 26). Pairwise LD between SNPs at each candidate locus and NHGRI cancer GWAS SNPs was assessed using the SNAP server (76). The NCBI Genome Decoration Page was used for ideogram construction. Published cancer-associated SNPs from the NHGRI GWAS Catalog and additional SNP markers at loci implicated in this study were imputed using the Michigan Imputation Server (ref. 77; 1,000 genomes Phase 1v3 reference, Shapeit2). Imputed markers were used for visualizing loci using LocusZoom Software (78). For LocusZoom plot construction, we used default recombination rates, gene locations, and NHGRI catalog data. LD values were generated from the 1,000 Genomes 2012 European population, and promoter and enhancer regions were determined as described above.

HEK293 cell lines were obtained from the American Type Culture Collection in 2015 and stored in liquid nitrogen vapor. Cells were authenticated using short tandem repeat fingerprinting within 6 months of freezing. HEK293 cells for the current study were thawed from this stock.

HEK293 control and HEK293 STK11 knockout cell lines were grown and maintained in DMEM supplemented with 10% FBS, 100 units/mL penicillin, 100 μg/mL streptomycin, and 250 ng/mL amphotericin B (Sigma Aldrich) at 37°C in humidified air with 5% CO₂. Cells were grown to 70% to 80% confluence prior to replating for transfection experiments. STK11 knockout cells were engineered by the CRISPR/Cas9 gene editing system, using the pSpCas9 (BB)-2A-Puro (PX459) V2.0 vector (purchased from Addgene #62988). The corresponding single-guide RNAs were designed using the CRISPR design website following published protocols (79). Cells were transfected with the CRISPR construct and selected for 3 days using puromycin then replaced by fresh medium to maintain the growth of the cell. Mass culture of cells was confirmed to lack STK11 protein expression by Western blotting (Supplementary Fig. S7A).

Twenty-four hours before transfection, HEK 293 cells were plated in 6-well plates at 40% confluence. When indicated, cells were transfected with siRNA targeting PTEN and control siRNA (purchased from Dharmacon, ON-TARGET plus SMART pool) at a final concentration of 10 nmol/L, using Turbofect transfection reagent (Thermo Fisher) according to the manufacturer's instructions. The following day, cells were transfected with increasing amounts of pcDNAlll–GNA11 where indicated, to control GNA11 protein expression levels at 0, 0.25, 0.5, 1, 2, and 4 μg/well. Cells lysates were harvested after 24 hours, serum starving for the last 4 hours. Reproducibility of PTEN siRNA knockdown was confirmed by Western blot (Supplementary Fig. S7B).

Immunodetection was carried out using antibodies from Cell Signaling Technology against PTEN (1:1,000), STK11 (1:1,000), phospho-STK11 (Ser428; 1:1,000), ribosomal protein S6 (1:20,000), phospho-S6 (Ser240/244; 1:5,000), AMPK (1:5,000), and phospho-AMPKα (Thr172; 1:5,000). Antibodies targeting GAPDH (1:20,000) and α-tubulin (1:3,000) were included as a loading control. GNA11 (D-17; 1:400) was detected using a primary antibody from Santa Cruz Biotechnology (sc-394). Secondary horseradish peroxidase–linked goat anti-rabbit IgG antibodies (1:20,000) were obtained from Southern Biotech.

A.M. Gross is a scientist at Illumina. No potential conflicts of interest were disclosed by the other authors.

Conception and design: H. Carter, M. Hofree, T. Ideker

Development of methodology: H. Carter, T. Ideker

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): X. Wu, E.L. Van Nostrand, J.S. Gutkind

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): H. Carter, R. Marty, M. Hofree, A.M. Gross, J. Jensen, K.M. Fisch, C. DeBoever, E.L. Van Nostrand, E. Wheeler, G.W. Yeo, J.S. Gutkind, T. Ideker

Writing, review, and/or revision of the manuscript: H. Carter, R. Marty, K.M. Fisch, J.F. Kreisberg, S.M. Lippman, J.S. Gutkind, T. Ideker

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): Y. Song, T. Ideker

Study supervision: T. Ideker

We thank all members of the Ideker Lab, as well as Cherie Ng, Kelly Frazer, Erin Smith, and Richard Kolodner for scientific feedback and discussion. The results presented here are in whole or part based upon data generated by the TCGA Research Network (http://cancergenome.nih.gov/).

This work was supported by NIH grants DP5 OD017937 to H. Carter, T32 GM008666 to C. DeBoever, and U54 CA209891 and U24 CA184427 to T. Ideker. Additional support was provided by UCSD CTRI grant UL1TR001442 to K.M. Fisch and the Fred Luddy Family Foundation.