Systematic identification of mutations and copy number alterations associated with cancer patient prognosis

1Cancers that arise from the same tissue can exhibit vast differences in clinical behavior. For instance, among individuals diagnosed with early-stage colorectal cancer, about 60% of patients will be cured by surgery alone, while the remaining 40% will experience a recurrence that is frequently fatal. Various pathological and molecular biomarkers are typically analyzed in order to assess patient risk and aid clinical decision-making. In general, these biomarkers are divided into two classes: predictive and prognostic. Predictive biomarkers identify patients who are likely to respond to specific therapies, like the EGFR mutations that sensitize lung tumors to EGFR inhibition. In contrast, prognostic biomarkers provide information on cancer aggressiveness and the likelihood of patient death. Tumor de-differentiation and lymph-node infiltration serve as prototypical prognostic biomarkers due to their strong association with poor outcomes. Yet, these pathology-based biomarkers can suffer from low levels of inter-observer concordance, and even perfect pathological assessment yields incomplete information on a patient’s most likely clinical course. New methods to identify aggressive tumors could lead to improvements in the stratification of patient risk, better clinical management, and a decrease in dangerous and unnecessary over-treatment.

2Advances in high-throughput technologies have yielded unprecedented insight into the diverse array of genomic changes found within every cancer cell. Projects like The Cancer Genome Atlas (TCGA) and the International Cancer Genome Consortium (ICGC) have characterized methylation, mutation, copy number, and gene expression patterns across cancer types. As a result of these studies, many of the genomic differences between normal and transformed cells have been identified and characterized. However, we lack a similar understanding of the genomic differences between indolent tumors and aggressive malignancies. As the cost of DNA sequencing continues to drop, it has become increasingly feasible for hospitals to implement routine targeted and/or genome-wide analyses of patient tumors. But, while several DNA-based, therapy-specific predictive biomarkers have been discovered, the prognostic information contained within tumor genomes is much less clear.

3Previous genome-wide efforts to discover novel prognostic biomarkers have largely focused on the gene expression changes associated with patient mortality. These studies have identified a set of transcripts that encode proteins involved in cell cycle progression that correlate with recurrence and death in several cancer types. Comparatively less is known about how changes at the DNA level affect patient survival. Outcome-associated analyses of genetic mutations have predominantly been conducted on a limited number of known oncogenes from single cancer types and have come to divergent conclusions. Reports in the literature commonly suggest that mutations in driver oncogenes are associated with poor outcomes, including, for instance, KRAS mutations in lung cancer, PIK3CA mutations in breast cancer, and BRAF mutations in colorectal cancer. Other studies of the same genes in the same cancer types have failed to observe any significant associations with outcome. In general, mutation-based biomarker studies may be confounded by small samples sizes, post-hoc hypothesis testing, imprecise clinical endpoints, and the so-called ‘file drawer’ problem, in which negative findings are less likely to be published. The prognostic information captured by sequencing driver oncogenes remains unknown, and a pan-cancer, exome-wide analysis of outcome-associated mutations has not been conducted.

4Previous investigations into the prognostic importance of DNA copy number alterations (CNAs) have indicated that highly-aneuploid tumors tend to have worse outcomes than diploid tumors. However, these analyses have largely focused either on arm-length changes or on alterations that affect single oncogenes or tumor suppressors. The functional importance of copy number changes in most genes at the single-gene level is unknown, and a pan-cancer, gene-by-gene analysis of prognostic copy number alterations has not been conducted. In order to gain a global understanding of the genomic features in a primary tumor that influence cancer prognosis, we collected and analyzed molecular profiles from a ‘discovery’ set of 9442 patients and a ‘validation’ set of 8618 patients with solid tumors. Our comprehensive, gene-centric analysis sheds light on the genomic changes that drive aggressive disease and will provide a useful resource for the development of strategies to improve clinical risk assessment. Additionally, we provide a web portal to facilitate community access to this rich biomarker dataset at http://survival.cshl.edu.

5To determine the differences between benign and fatal tumors, we first analyzed multiple classes of genomic data from 9442 patients with 16 types of cancer from the TCGA (outlined in Figure 1—figure supplement 1A; abbreviations are defined in Figure 1—figure supplement 1B). For every tumor type and every dataset, we generated Cox univariate proportional hazards models linking the presence or expression of a particular feature with clinical outcome (described in Supplemental Text 1). We report the Z score for each model, which encodes both the directionality and significance of a particular association. If the presence of a mutation or copy number amplification is significantly associated with patient death, then a Z score >1.96 corresponds to a P value < 0.05 (Figure 1—figure supplement 2A–C). In contrast, a Z score less than −1.96 indicates that the presence of a mutation is associated with survival or that a gene deletion is significantly associated with patient death.

6We extracted mutation, copy number, gene expression, and clinical information from 16 TCGA cohorts (summarized in Supplementary file 1 and discussed in additional detail in Supplemental Text 2). To assess the validity of our data analysis pipeline, as well as the accuracy of the reported patient annotations, we first examined the overall survival curves for the 16 tumor types that we profiled. As expected, we observed significant differences in clinical outcome according to a cancer’s tissue-of-origin (Figure 1—figure supplement 2D). Prostate cancer had the least aggressive clinical course, with a median survival time that was not reached in this dataset (>4600 days), while pancreatic cancer conferred the worst prognosis (median survival time: 444 days). Overall, the 5 year survival frequencies of patients in the TCGA were highly similar to the national averages reported by NCI-SEER (R = 0.83, p < 0.0001), suggesting that the patients included in this analysis are broadly representative of the general population (Figure 1—figure supplement 2E). Next, we inferred patient sex on the basis of chromosome-specific gene expression patterns. Our analysis exhibited >99% concordance with patients’ self-reported sex, further verifying the overall accuracy of the clinical annotations and our data processing pipeline (Figure 1—figure supplement 2F).

7We first set out to discover whether coding mutations in cancer genomes were associated with patient outcome. We extracted non-silent mutations in each tumor, and then we identified all genes that were mutated in ≥2% of patients in each of the 16 cohorts (discussed in Supplemental Text 1). We next performed Cox proportional hazards analysis to compare the survival times for patients harboring mutant or wild-type copies of each gene. This analysis uncovered very few mutations that were significantly associated with patient outcome (Figure 1 and Supplementary file 2A-B). We first focused on known oncogenes and tumor suppressors, and found that among the 30 most-frequently mutated cancer driver genes, only two (EGFR and TP53) were associated with prognosis in more than two tumor types (Figure 1C). TP53 mutations were linked to outcome in five of 16 cancer types, though the differences in patient survival were generally small (Figure 1—figure supplement 3A–B). In contrast, many other cancer driver genes were not associated with survival time in any tumor type. While mutations in KRAS, PIK3CA, CDKN2A, BRAF, KMT2D, ATM, SMAD4, and many other genes were frequently observed, they were never significantly linked with patient outcome (Figure 1C).

8We next considered the possibility that mutations in specific codons could have prognostic significance not captured when all mutations in a gene are pooled together. To test this, we identified the 30 most-frequently mutated amino acid positions in the TCGA cohorts, and then asked whether patients harboring these alterations had different outcomes than those who did not. IDH1c132 mutations were significantly associated with a favorable prognosis in glioma, but other recurrently-mutated codons (KRASc12, PIK3CAc1047, TP53c273, etc.) were largely uninformative (Figure 1—figure supplement 4A–D). Then, we identified ‘hotspot’ residues that were mutated in at least five different patients across all cohorts. Considering only these ‘hotspot’ mutations in each gene also failed to uncover robust survival associations (Figure 1—figure supplement 4E). Finally, we identified cancer type-specific recurrent mutations, but these alterations (FGFR3c249 in BLCA, CTNNB1c37 in UCEC, etc.) were similarly uninformative (Figure 1—figure supplement 4F).

9Next, we sought to test whether the use of targeted therapies had blunted the deleterious effects of certain driver mutations (e.g., in BRAF or EGFR). However, due to the time-frame of sample collection, very few patients were treated with BRAF or EGFR inhibitors, and removing those patients who had received these therapies failed to significantly affect Z scores (Figure 1—figure supplement 5A). Hyper-mutation within a subset of cancers could increase mutational ‘noise’ and decrease our ability to identify prognostic signatures, but excluding patients with hyper-mutated tumors had minimal effect on the prognostic significance of driver gene mutations (Figure 1—figure supplement 5B and Supplementary file 2C). We then asked whether the presence of mutations in multiple cancer driver genes might cooperate to confer a worse clinical outcome. We found that, in general, patients harboring mutations in two cancer driver genes that were not prognostic alone had the same risk of death as patients with wild-type copies of one or both genes (Figure 1—figure supplement 5C). Lastly, we considered the possibility that the clonality of a mutation might affect its prognostic significance. We calculated the variant allele frequency (VAF) for each cancer mutation and tested whether mutations present at clonal levels in single tumors were more likely to be associated with outcome. We found that restricting our analysis to mutations with high VAFs failed to identify more prognostic genes, indicating that patient stratification is unlikely to be improved by assessing only clonal mutations (Figure 1—figure supplement 6).

10These analyses suggested that, in general, cancer driver gene mutations lacked significant patient stratification power. This led us to investigate whether mutations in genes other than recurrently-mutated oncogenes and tumor suppressors could affect prognosis. We therefore expanded our analysis to include all genes mutated in ≥2% of patients with a particular tumor type. To account for greatly expanding the number of genes tested, we applied a Benjamini-Hochberg correction with a 5% false-discovery rate to the individual Z scores that we obtained. We uncovered several genes that were linked with prognosis in glioma, but found very few genes significantly associated with death or survival in the other 15 cancer types (Figure 1D and Supplementary file 2A). For instance, in breast cancer and lung adenocarcinoma, 128 and 3996 genes were mutated in ≥2% of patients, respectively, but none of these mutations were significantly correlated with patient outcome at a 5% FDR. In total, these results indicate that most mutations in cancer genomes lack significant prognostic power.

11In our above analysis, we noted that the five genes with the strongest survival associations were all observed in the GBMLGG (pan-glioma) cohort. As glioma appeared to be an exception to our overall finding that mutations are seldom prognostic, we investigated this cohort further. Among the top-scoring genes, we found that PTEN and EGFR mutations conferred dismal prognosis, while mutations in IDH1, TP53, and ATRX were associated with favorable prognosis (Figure 1—figure supplement 7A). Mutations in these genes have previously been linked to distinct glioma subtypes, and we verified that mutations in IDH1, TP53, and ATRX were most frequently observed in low-grade gliomas, while mutations in PTEN and EGFR were most frequently observed in high-grade glioblastomas (Figure 1—figure supplement 7B). However, when we analyzed low-grade gliomas and glioblastomas separately, several of these alterations remained prognostic (Supplementary file 2D). For instance, while IDH1 mutations were more common in low-grade gliomas, they were occasionally observed in high-grade tumors as well, and they were independently associated with prolonged survival in both cohorts (Figure 1—figure supplement 7C). In contrast, when EGFR mutations were observed in low-grade gliomas, they were associated poor outcomes, but EGFR mutations were non-prognostic in high-grade glioblastomas (Figure 1—figure supplement 7D). Thus, in gliomas, mutations contain both subtype-dependent and subtype-independent prognostic information. However, outside of this cancer type and the tumor suppressor TP53, mutations in most cancer driver genes are non-prognostic.

12As mutations were largely uninformative, we next set out to determine whether gene copy number conveyed prognostic information. We determined the copy number of each gene at its transcriptional start site and regressed this value against patient outcome in each tumor cohort. We then examined the clinical impact of CNAs affecting the same 30 cancer driver genes that we previously investigated. Surprisingly, we found that the copy number of these oncogenes and tumor suppressors was frequently linked with patient outcome (Figure 2 and Supplementary file 3A-B). Amplification of EGFR, PIK3CA, and BRAF, and deletion of CDKN2A, RB1 and EP300 were strongly associated with shorter patient survival times in four or more cancer types each. Copy number was prognostic even for genes in which mutations were not linked with outcome: for instance, while mutations in PIK3CA were never informative, the copy number of PIK3CA was associated with outcome in breast, colorectal, glioma, lung-squamous, pancreas, and prostate cancers (Figure 2B and D). Overall, among the 30 most frequently-mutated cancer driver genes, we detected 108 significant associations between gene copy number and outcome, compared to 23 associations between mutation and outcome. For 28 out of 30 driver genes, DNA copy number was prognostic in more cancer types than mutational status was. We conclude that determining the copy number of oncogenes and tumor suppressors in a primary tumor can better stratify patient risk than assessing single base-pair mutations.

13In our analysis thus far, we have treated mutations as a binary variable (‘mutant’ vs. ‘not mutant’), while copy number alterations are treated as continuous values. Thus, the greater prognostic significance of tumor CNAs could reflect the fact that individual CNA measurements inherently harbor more information. To test this possibility, we trichotomized CNA values into ‘deletions’ (<−0.3), ‘amplifications’ (>0.3), and ‘copy-neutral’ (≥−0.3 and≤0.3). We then calculated Cox regressions at the same 30 loci using the discretized copy number values. This analysis resulted in 94 significant survival associations, more than four times as many significant features as when mutations were analyzed, and comparable to the number of significant features that resulted using continuous CNA values (Figure 2—figure supplement 1). This analysis suggests that the greater prognostic significance of CNAs is not simply a consequence of the continuous nature of copy number data.

14We next investigated whether these oncogene and tumor suppressor CNAs were likely to drive patient mortality, or whether they were passenger genes that changed in copy number along with other, unknown drivers. To assess this question, we combined Z scores from different cancer types using Stouffer’s method, and then plotted the pan-cancer meta-Z scores along every chromosome (Figure 2C). This analysis revealed multiple sharp peaks and valleys in the data that overlapped with known driver mutations. The most significant survival-associated copy number changes genome-wide were found on chromosome 9p in a valley that precisely included the tumor suppressor CDKN2A. Z score peaks were found at loci that include oncogenes PIK3CA, EGFR, MYC, CCNE1, and others. This overlap suggests that, in many instances, the copy number of these oncogenes and tumor suppressors directly influence the risk of cancer patient death.

15Whole-chromosome aneuploidy has previously been linked to a decreased infiltration of immune cells. We therefore considered the possibility that CNAs are prognostic via an indirect mechanism; namely that they are found in tumors that lack robust immune infiltration, and this deficient immune response was itself driving patient mortality. However, multiple lines of evidence argue against this interpretation. First, we assessed the association between patient survival and three different measures of tumor sample purity: pathologist-assessed tumor cell fraction, sample purity as judged by ABSOLUTE, and leukocyte infiltration, as judged by methylation analysis. We found that sample purity was inconsistently-associated with patient outcome (Figure 2—figure supplement 2). For instance, higher tumor purity determined by either pathological analysis or ABSOLUTE was associated with worse outcome in only one of 16 cohorts, each (Figure 2—figure supplement 2A–B). The lack of a strong correlation between infiltrating cell populations and clinical prognosis suggests that analyte purity is insufficient to explain the relationship between CNAs and patient survival. Secondly, we generated multivariate Cox models that included gene copy number and these three measurements of tumor purity, and we found that driver gene CNAs remain broadly prognostic in these bivariate models (Figure 2—figure supplement 2C and Supplementary file 3C-E). For instance, we discovered that amplification of Cyclin E1 is associated with poor prognosis in ovarian cancer, and this remained true even when our analysis was restricted to high-purity tumor samples and samples that lacked significant leukocyte presence (Figure 2—figure supplement 2D). Thus, while the interrelationship between aneuploidy and immunological tolerance likely plays an important role in tumor development, this analysis suggests that it is not the primary driver of CNA-associated patient mortality.

16Pathological assessment of tumor stage and grade are important sources of prognostic information, though blinded assessments reveal significant inter-observer discordance. We therefore tested whether the CNA biomarkers that we uncovered could affect the stratification conferred by these parameters. We found that Z scores generated from either univariate models or multivariate models that included stage or grade were highly correlated (R = 0.91 and R = 0.96 respectively; Figure 2—figure supplement 3A and Supplementary file 4). Overall, 71% of prognostic CNAs in individual cancer types remained prognostic in these multivariate models (Figure 2—figure supplement 3B). Thus, including gene-level copy-number assessment can significantly improve the stratification of patient risk beyond standard clinical parameters (Figure 2—figure supplement 3C–D). Certain mutations were similarly able to yield prognostic information in a stage- and grade-independent manner. However, due to the lower overall significance of most mutations that we identified, the improvements in patient stratification were generally more modest (Figure 2—figure supplement 3E and Supplementary file 4E).

17Gene-level copy number values also remained prognostic when separating TCGA cohorts by cancer subtype (Figure 2—figure supplement 4 and Supplementary file 4F-G). For instance, CNA Z scores values were highly correlated between the bulk GBMLGG cohort and the individual GBM and LGG subtypes (R = 0.68 and R = 0.86, respectively; Figure 2—figure supplement 4A). While analyzing the GBM cohort separately abolished the prognostic significance of EGFR mutations (Figure 2—figure supplement 4D), EGFR amplifications remained associated with outcome in both the LGG and GBM cohorts (Figure 2—figure supplement 4B–C). Amplifications in MYC and PIK3CA were similarly prognostic in multiple tumor subtypes (Figure 2—figure supplement 4D–E and Supplementary file 4G). At other loci, low patient numbers from certain subtypes may obscure the detection of specific biomarkers. For instance, within the KIPAN cohort, 67% of tumors are clear cell carcinomas, 23% of tumors are papillary cell carcinomas, and 10% of tumors are chromophobe carcinomas. CDKN2A deletion is a strong indicator of poor prognosis in the pan-kidney cohort, in clear cell carcinomas, and in papillary cell carcinomas, but did not reach statistical significance when kidney chromophobe carcinomas were analyzed independently (Figure 2—figure supplement 4F–G). In total, these results underscore the ability of driver gene CNAs to improve patient stratification when controlling for tumor identity, though larger cohort numbers may be needed to identify the strongest biomarkers in rare cancer subtypes.

18Highly-aneuploid tumors tend to harbor mutations in TP53, and both TP53 mutations and arm-length aneuploidy have previously been associated with poor clinical outcomes. Using an ‘aneuploidy score’ for each tumor based on the total number of arm-length alterations, we verified that TP53-mutant tumors exhibit more aneuploidy than TP53-wild-type tumors (Figure 3—figure supplement 1A), and that total aneuploidy is a poor prognosis factor in several cancer types (Figure 3—figure supplement 1C). To investigate the relationship between gene-level prognostic CNAs, TP53 status, and arm-length aneuploidy, we selected a set of 40 prognostic amplifications and deletions for additional analysis (Figure 3—figure supplement 2A). In multivariate models that included TP53 mutation status, 33 of 40 (83%) gene-level CNAs remained prognostic, demonstrating that these CNAs are not linked with death due to an indirect association with TP53 status (Figure 3—figure supplement 2A–B). Similarly, in multivariate models that included total tumor aneuploidy, 80% of these CNAs were still associated with outcome (Figure 3—figure supplement 2C–D). Finally, as a proxy for the total structural alteration burden, we summed the number of breakpoints (as indicated by discrete copy number values along a chromosome) in each tumor (Figure 3—figure supplement 1B). This metric was associated with outcome in multiple tumor types (Figure 3—figure supplement 1C), but 75% of driver gene CNAs remained prognostic in multivariate models that included this score (Figure 3—figure supplement 2E–F). These results indicate that assessing gene-level tumor CNAs can yield more prognostic information than simply screening for TP53 mutations or measuring bulk levels of tumor aneuploidy (Supplementary file 5).

19We next set out to determine whether focal copy number alterations and broad copy number alterations could have distinct effects on patient outcome. To investigate this possibility, we compared the prognostic power of focal CNAs (defined as an alteration ≤3 Mb in length; ) and broad CNAs (defined as all alterations >3 Mb in length). Among loci at which both broad and focal alterations were observed, we frequently found that broad CNAs were associated with moderately worse outcomes, while focal CNAs were associated with sharp declines in survival (Figure 3A–B). At some loci, broad CNAs had outcomes that were indistinguishable from copy-neutral tumors, while only focal CNAs were associated with death (Figure 3C). We rarely detected instances in which broad CNAs indicated a worse prognosis than a focal alteration (Figure 3A). We interpret these results as a reflection of aneuploidy-induced fitness penalties: large copy number alterations change the dosage of multiple genes at once and can impair tumor growth, while targeted alterations that specifically affect driver gene copy number maximize malignant potential.

20Gene copy number alterations typically result in a proportional change in the expression of the affected loci, though instances of dosage compensation have been reported. To test the effects of prognostic CNAs on gene expression, we compared transcript levels and gene copy number changes at 40 prognostic loci and found a significant correlation between the two at 98% of the analyzed genes (Figure 3—figure supplement 3). Next, we sought to uncover whether these copy number alterations were deadly because they increased or decreased the expression of mutant gene products. That is, we could observe that the amplification of a driver gene is prognostic only in tumors in which that driver gene is also mutated. Interestingly, this is not the case: at 95% of our test loci, gene copy number remained prognostic in multivariate models that also included gene mutation status (Figure 3D). For instance, in colorectal cancer, amplification of EGFR was associated with death even in tumors that lacked EGFR mutations (Figure 3E). In total, these results indicate that even at recurrently-mutated loci, changes in the expression of the wild-type gene can have a profound effect on cancer cell behavior. Together with our observation that focal changes tend to confer a worse prognosis than broad changes, these results support the recently-proposed ‘cancer gene island’ model of tumor genome evolution (discussed in more detail below; ).

21To determine the generality of our findings, we collected independent patient cohorts harboring mutation or copy number data linked to survival outcome (Supplementary file 1). We then performed univariate Cox proportional hazards analysis on these ‘validation’ cohorts and compared the results to the Z scores obtained from our ‘discovery’ set of TCGA data. First, we identified prognostic mutations within a set of 16 patient cohorts from the International Cancer Genome Consortium (ICGC), comprising 3054 patients analyzed by whole-genome or whole-exome sequencing. Overall, the mutation frequencies and the Z scores of recurrent mutations were highly similar between the ICGC and TCGA cohorts (R = 0.67, p < 0.0001, and R = 0.56, p < 0.0001, respectively; Figure 4A–B). Consistent with our TCGA analysis, mutations in TP53 were associated with outcome in more patient cohorts than any other gene (Figure 4C and Figure 1—figure supplement 3C–D). Other mutations, including in known cancer driver genes, were rarely associated with outcome in individual cancer types and harbored minimal pan-cancer significance (Figure 4D–F and Supplementary file 6). Mutations in KRAS, PIK3CA, BRAF, APC, PTEN, CDKN2A, and many others were frequently observed but were never correlated with outcome (Figure 4D). We next analyzed 2431 additional patients with CNA data curated by cBioportal, and found numerous amplifications and deletions associated with patient mortality (Supplementary file 6C). In breast cancer, we found prognostic amplifications that were centered around oncogenes, including ERBB2, MYC, and MDM2, while prognostic deletions encompassed tumor suppressors CDKN2A, PTEN, and TP53 (Figure 4G). Overall, we observed a highly significant correlation between the meta-Z scores obtained from the TCGA and cBioportal datasets (R = 0.42; Figure 4H). Finally, in patient cohorts subjected to both mutation and copy number analysis, we verified that CNAs in driver genes commonly harbored greater prognostic significance than mutations in those same genes (Figure 4I). For instance, in breast cancer, among 25 frequently-mutated genes, mutations in only two genes (TP53 and GATA3) displayed prognostic significance, while CNAs in 12 of those same genes were associated with patient outcome (Figure 4J). In total, these analyses suggest that the survival patterns discovered in the TCGA dataset are conserved across independent cohorts of cancer patients. In particular, while mutations in most cancer driver genes are non-prognostic, copy number alterations in these same genes are tightly linked with patient outcome.

22In order to discover the biomarkers with the greatest potential clinical relevance, we next identified the individual mutations and CNAs that were consistently associated with outcome across independent patient cohorts. To increase our ability to detect these genetic alterations, we performed survival analysis on an additional set of 2701 primary tumors subjected to targeted sequencing and copy number analysis (MSKCC_2017; Supplementary file 7), on 2431 patients from cBioportal cohorts whose tumors had been sequenced (Supplementary file 6D), and on 628 patients from ICGC cohorts subjected to copy number analysis (Supplementary file 6B). Our combined patient dataset therefore included two to six independent cohorts from each of 13 common cancer types, comprising 16,580 total patients. These cohorts were collected at different locations, in different patient populations, using different study designs, and the samples were analyzed using different genomic technologies. We reasoned that alterations that were consistently associated with outcome despite these significant differences would represent highly-penetrant biomarkers of patient prognosis. To identify such alterations, we screened for biomarkers that were associated with outcome (|Z| > 1.96) in ≥2 independent cohorts, and that were highly significant (|meta-Z| > 3.3) across all available cohorts. This approach revealed multiple high-confidence genetic biomarkers of patient outcome that, to our knowledge, were novel, including MDM4 amplifications in prostate cancer, NOTCH2 amplifications in melanoma, and 2q32 deletions in ovarian cancer (Supplementary file 8). These robust biomarkers allowed a striking stratification of patient risk, and top-scoring CNAs remained prognostic in multivariate models that included commonly-measured prognostic criteria (Gleason score in prostate cancer, Hepatitis serology in liver cancer, etc.; Figure 5—figure supplement 1). Consistent with our single-cohort analyses, cross-cohort prognostic CNAs were significantly more common than prognostic mutations, and TP53 was the only gene whose mutation status was associated with outcome in more than one cancer type (Figure 5—figure supplement 2A).

23We hypothesized that some genetic alterations that were sufficient to affect overall patient survival could impact other facets of cancer behavior as well, including, potentially, drug sensitivity. That is, biomarkers harboring significant prognostic information could potentially contain predictive information as well. We therefore sought to discover whether genetic alterations that drove aggressive disease could also sensitize patient tumors to specific therapeutic regimens. By analyzing a cohort of 1000 patient-derived xenografts (PDXs), we identified several instances in which high-confidence biomarkers were associated with vulnerability to particular anti-cancer agents. For instance, we identified Chr9 deletions that encompassed CDKN2A as a robust biomarker for poor prognosis in breast cancer (Supplementary file 8). We found that PDXs harboring CDKN2A deletions were profoundly sensitive to combination therapy with a CDK4/6 inhibitor and an mTOR inhibitor (Figure 5—figure supplement 2B), consistent with the fact that a protein encoded by CDKN2A, p16, functions as a natural inhibitor of CDK4/6, p. 4). In contrast, other biomarkers associated with poor prognosis in breast cancer failed to predict sensitivity to this treatment combination, but instead correlated with sensitivity to other agents (Supplementary file 8). Due to the limited number of drugs tested in PDXs, we expanded our target search to include a recently-described pharmacogenomic profile of cancer cell lines and discovered several additional biomarker vulnerabilities (Figure 5A–B). For instance, we identified mutations in STAG2 as a high-confidence biomarker of poor prognosis in glioma, and we found that STAG2-mutant gliomas were exquisitely sensitive to treatment with the PARP inhibitor olaparib (Figure 5A). In total, we identified highly-significant therapeutic vulnerabilities for 49% of the prognostic biomarkers uncovered by our integrated analysis, providing potential strategies to treat a subset of patients who have the most aggressive cancers.

24Modern medicine has vastly prolonged the survival of individuals diagnosed with cancer. However, increasing evidence suggests that large subsets of patients receive sub-optimal care, and are over-treated or under-treated relative to their level of risk. To date, many of the genetic alterations that differentiate fatal and benign tumors have remained obscure. Our analysis of prognostic biomarkers from 17,879 patients sheds light on these genetic differences, identifies a subset of patients who may benefit the most from aggressive intervention, and suggests therapeutic strategies for tumors harboring certain alterations associated with poor prognosis. A web portal to facilitate access to these results is available at http://survival.cshl.edu/.

25As cancers arise due to the accumulation of mutations in growth-promoting oncogenes and growth-inhibitory tumor suppressors, the presence and diversity of these mutations may be expected to dictate a tumor’s clinical course. However, our data suggest that in many cases, they do not. Substantial disagreements exist in the literature on the value of mutation-based prognostic biomarkers, as the same driver oncogenes have been independently reported to be either adverse or non-significant prognostic features. In this manuscript, we performed an unbiased genome-wide analysis of public datasets with pre-established sample sizes. This approach may therefore bypass certain problems, including post-hoc hypothesis testing, patient-selection bias, and the ‘file-drawer problem’, that can confound targeted biomarker studies. We consider it possible that, with larger sample sizes or more-specific tumor subtypes, additional prognostic mutations could be identified. Importantly, in most patient cohorts that we collected, tumors were analyzed on multiple genomic platforms, and CNAs were commonly prognostic in the same cohorts in which gene mutations were not. These results underscore our ability to successfully detect biomarkers in cohorts of these sizes, and suggest that, in a head-to-head comparison, copy number alterations provide more useful prognostic information than single-gene mutations.

26While we identified very few mutations associated with patient outcome, several lines of evidence underscore the potential benefits of continued clinical sequencing efforts. First, our analysis revealed a subset of mutations with tissue-specific prognostic power, including TP53 mutations in breast cancer, RB1 mutations in bladder cancer, and FBXW7 mutations in colorectal cancer. Secondly, most patients in the TCGA cohorts were treated with standard cytotoxic drugs. As targeted and immuno-therapies are increasingly adopted in the clinic, oncogenic mutations that were non-prognostic in the datasets analyzed here may be able to predict sensitivity to specific therapeutic agents. Thirdly, tumors themselves are composed of sub-clonal populations that harbor distinct sets of mutations, and recent evidence suggests that cancer heterogeneity can influence clinical course. Thus, interrogating the mutational spectrum at the sub-clonal level may identify prognostic mutations not distinguished in bulk analyses.

27Though large-scale changes in tumor ploidy have previously been recognized as an indicator of poor outcome, the contributions of copy number alterations in most single genes have remained unexplored. Despite the limited stratification value of mutations in cancer driver genes, we found that copy number alterations of many of these same genes are broadly prognostic. Focal CNAs tended to confer a worse prognosis than broad CNAs, consistent with a model in which large-scale gene dosage imbalances trigger proteotoxic stress and impose a fitness penalty on cancer cells. Moreover, while prognostic CNAs commonly caused proportional changes in target gene expression, most CNAs remained prognostic whether or not they affected the expression of a mutated gene. These results support a ‘cancer gene island’ or ‘cumulative aneuploidy’ model of tumorigenesis, in which cancers accumulate a series of limited copy number changes affecting haplo-sensitive and triplo-sensitive regions. Identifying the functional consequences of these prognostic CNAs on tumor physiology is a key future goal.

28Patients whose tumors harbor genetic alterations that drive mortality are in urgent need of improved treatment options. We discovered many instances in which high-confidence biomarkers of aggressive disease also sensitized tumors to specific anti-cancer therapies. By taking advantage of these vulnerabilities, a precision-medicine approach could be applied to both stratify patient risk and identify drug combinations most likely to provide a clinical benefit. Several predicted sensitivities from our work have clinical or mechanistic support, including the use of CDK4/6 inhibitors to treat CDKN2A-deleted tumors, the use of PARP inhibitors to treat STAG2-mutant tumors, and the use of SYK inhibitors to treat RB1-mutant tumors. Treatment with targeted agents significantly alters the cellular epigenetic and genetic landscape, often culminating in the development of resistance to the applied therapies. We speculate that secondary alterations that tumors evolve to tolerate these drugs could also alter or blunt the aggressive phenotype caused by the original driver alteration. In this way, targeting a biomarker that confers poor prognosis could both directly lead to improved patient outcomes by triggering a robust clinical response, and indirectly help patients by forcing tumor evolution away from dependence on a driver of aggressive disease.

29Patient cohorts analyzed in this study are listed in Supplementary file 1. For the TCGA analysis, pre-processed files from the Broad Institute TCGA Firehose were used (https://gdac.broadinstitute.org/). For the TCGA genomic copy number analysis, we used the HG19 segmented SCNAs, corrected for germline SCNAs. Overall survival time was used as a clinical endpoint for all cancer types except PRAD. Overall survival was chosen because it reflects an objective and unambiguous event, it is the gold-standard for oncology clinical trials, and it is widely-available across different studies. However, as fewer than 2% of the patients in the PRAD cohort died during the follow-up period, ‘days to biochemical recurrence’ was used as a surrogate endpoint. For all cancers, survival or follow-up time from diagnosis were corrected for the days to sample procurement. Primary tumors (indicated with a ‘01’ in the patient barcode) were used for every cancer type except SKCM; for this cancer, few primary samples were available, so metastatic samples (indicated with a ‘06’ in the barcode) were included for patients in which no primary tumor was available. For additional discussion of the TCGA samples, see Supplemental Text 2. Pathology-assessed tumor cell fraction was obtained from the TCGA clinical files under ‘Percent_tumor_cells’. Tumor stage and grade were similarly obtained from the appropriate TCGA clinical files.

30Mutation, copy number, and clinical data from Release 25 of the International Genome Consortium were downloaded from the ICGC Data Portal. Overall survival was used as a clinical endpoint for all cohorts except EOPC-DE; due to the few deaths in this cohort, recurrence-free survival was used as an endpoint. Cohorts were chosen based on the availability of WGS or WES data, and were included if they came from a cancer type comparable to the types that were studied in our TCGA analysis.

31Copy number, mutation, and clinical data from cBioportal were downloaded as pre-processed files from www.cbioportal.org. For the patients described in (the cBioportal/MSKCC_2017 cohorts), only primary tumors were included for all cancer types except melanoma.

32All processing and analysis was performed using Python. Cox proportional hazard analysis used the R survival package (https://cran.r-project.org/web/packages/survival/index.html) to compute Z scores and p values. Justification and further explanation for the use of Cox proportional hazards modeling can be found in Supplemental Text 1. The rpy2 project was used to control R from python, allowing seamless integration of Z score calculations with data processing and pan-cancer analysis. Pandas DataFrames were used as the primary structure for storing and manipulating data. Additionally, native numpy methods and arrays were for used occasionally for efficiently storing strictly numerical data, for example, as input to Cox proportional hazards models. The statsmodels package (www.statsmodels.org) was used for false discovery correction using the Benjamini-Hochberg procedure. Microsoft Excel was occasionally used for final data processing and examination, so a single apostrophe was added before gene names in intermediate data processing steps to protect genes from auto-formatting.

33Code was structured to allow ease of internal reuse and reproducibility of results. Cox univariate proportional hazards, Cox multivariate proportional hazards, Kaplan-Meier, and Stouffers analysis methods were factored into an analysis library, taking as input the data required to perform the computation as numpy arrays or pandas DataFrames.

34In addition to the code for statistical analyses, code for processing TCGA clinical files was factored into a common library. This approach allowed the same TCGA clinical file processing code to be executed across a variety of platform analyses, ensuring identical behavior for each platform. The TCGA clinical processing code selected the relevant clinical endpoints and sample procurement data. The processing translated the available clinical data into the required format for Cox proportional hazard models: an endpoint/survival time value and a censor value for each patient. Code to select tumor samples based on cancer type was also included in this library.

35Raw input data for the mutation analysis needed additional preprocessing before Cox proportional hazard models could be constructed. This preprocessing included removing per-patient headers throughout the data and some data transposition. For all analyses using TCGA mutation data, mutations annotated as silent were excluded. Genes were only included in downstream analyses if they were mutated in 2% or more of the patients in a cancer type cohort.

36Raw input data for copy number analysis also required substantial preprocessing. Copy number input data consists of per-patient, per-chromosome location maps of copy numbers (hg19 downloaded from the UCSC Genome Browser; ). These maps were converted to a single copy number value for each gene. We created an interval tree (using the intervaltree python package, https://pypi.python.org/pypi/intervaltree) of the location maps for each chromosome and used the appropriate HGNC to convert chromosome locations to genes for each patient. We used the gene’s transcriptional start site position to look up in the interval tree the copy number value for a gene. This analysis produced an intermediate file of a similar form to the other TGCA platforms, which allowed for straightforward Cox analysis. Note that Cox proportional hazards models are a threshold-independent method of performing survival analysis, and so no minimum or maximum threshold for a copy number alteration was specified.

37A tumor was defined as having a focal amplification or deletion if its copy number was greater than 0.3 or less than −0.3, and the chromosomal interval with a copy number greater than 80% of the copy number at the gene of interest was less than or equal to 3 Mb.

38To calculate the number of structural alterations per tumor, the number of distinct copy number values per chromosome in the DNA segmentation file was summed for each patient.

39For each platform and analysis type, we performed a pan-cancer analysis. This analysis created a single Z score for each gene by combining the per gene Z scores from each cancer type using Stouffer’s method. To perform Stouffer’s method, we took the sum of the Z scores for a single gene and divided that sum by the square root of the number of cancer types with Z scores for the gene. This meta-Z score was then compared against meta-Z scores obtained similarly from other platform analyses.

40We performed several additional analyses on mutation data, including double mutation combination Z scores, hotspot codon Z scores, and Z scores corrected for VAFs. For double mutation Z scores, we took the top 30 most common cancer driver genes and performed pairwise combinations. We then calculated Cox proportional hazards for each pair of genes, where a patient was considered to have a pairwise mutation if and only if both genes were non-silently mutated for that patient. Z scores were only calculated for a pair if (1) neither gene in the pair was statistically significant alone in the univariate analysis and (2) if both genes were mutated together in at least 10 patients.

41Per-codon Z scores were calculated for a selected set of hotspot codons. Most cancer types were available in HG37, so HG37 mutation positions were used to locate codons. Mutations for OV and COADREAD were only available in HG36, so gene positions were converted to HG37 before codon processing. Per codon Z scores were calculated by first identifying patients with mutations in the relevant gene, then selecting from that set of patients those whose mutations were in the codon of interest. If 2% of patients or more had mutations in the selected codon, a Z score was calculated.

42VAFs were calculated for 10 of the TCGA cancer types. We analyzed VAF data in two ways. First, we calculated Z scores, only counting a gene as mutated if its VAF was greater than or equal to 0.4. Secondly, we identified the median VAF score per gene, and calculated Z scores only counting a gene as mutated if its VAF was equal to or greater than the median VAF for that gene.

43CBioPortal was structured similarly to the TCGA analyses, though data processing was not factored into an independent library since each of these datasets was only used in one analysis. Copy number data from one CBioPortal cancer type, blca_mskcc, required initial preprocessing in the manner described above for TCGA copy numbers. Mutations were included if they were annotated as one of these types: In_Frame_Ins, Nonstop_Mutation, Translation_Start_Site, In_Frame_Del, Splice_Region, Frame_Shift_Ins, Frame_Shift_Del, Splice_Site, Nonsense_Mutation, or Missense_Mutation.

44ICGC analysis was structured similarly to CBioPortal analysis. Mutations were only included in downstream analyses if they were annotated as one of these types: disruptive inframe deletion, disruptive inframe insertion, frameshift variant, inframe deletion, missense variant, splice acceptor variant, splice donor variant, stop gained, or stop lost. Z scores were calculated if a gene was mutated in 2% or more of the patients in a particular cohort.

45Across independent datasets, cohorts of patients from related cancer types were identified. Mutations or CNAs significantly associated with patient prognosis (Z > 1.96 or Z < −1.96) in two or more independent cohorts from each cancer type were determined. Then, the subset of these alterations that remained highly-significant (Z > 3.3 or Z < −3.3) across all cohorts from the same cancer type were classified as high-confidence biomarkers. In some instances, amplifications that spanned continuous chromosomal regions were found to correlate with patient prognosis. These segments were identified manually. For the determinations of therapeutic sensitivity described below, the gene with the minimum meta-Z score (for deletions) or maximum meta-Z score (for amplifications) within a segment was chosen to represent the segment as a whole.

46Therapeutic sensitivity data for PDXs was acquired from. To identify mutations that correlated with therapy sensitivity, for each drug or drug combination, a comparison was performed if five or more PDXs had a mutation in a gene of interest, and if five or more PDXs were wild-type for a gene of interest. For genes and therapies fitting these criteria, we next identified instances in which the therapy resulted in a clinical response in the mutant population, defined as an average ‘Best Average Response’<15% tumor growth among PDXs with a mutation in the gene of interest. Finally, for genes and therapies fitting these criteria, we performed a t-test for the ‘Best Average Response’ between PDXs with mutant and wild-type copies of a gene of interest. We reported therapies in which these criteria were met and tumors with mutation were more sensitive to the therapy than tumors with wild-type copies of the gene of interest (p < 0.01).

47To identify CNAs that correlated with therapy sensitivity in the PDX cohort, amplifications and deletions (CNA >|.3|) were called, and then considered separately. As above, CNAs were included if five or more PDXs exhibited an alteration, and if five or more PDXs did not exhibit that alteration. For genes and therapies fitting these criteria, we next identified instances in which the therapy resulted in a clinical response in the altered population, defined as an average ‘Best Average Response’<15% tumor growth among PDXs with an amplification or deletion in the gene of interest. Finally, for genes and therapies fitting these criteria, we performed a t-test for the ‘Best Average Response’ between PDXs with mutant and wild-type copies of a gene of interest. We reported therapies in which these criteria were met and tumors with a mutation were more sensitive to the therapy than tumors with wild-type copies of the gene of interest (p < 0.01).

48Therapeutic sensitivity data from cancer cell lines was acquired from. For this data, two different comparisons were used. First, the calculations described below were performed for cell lines from the specific cancer type that the high-confidence biomarker was identified in. If this analysis yielded no significant vulnerabilities, then the calculations were repeated across all cancer types (pan-cancer).

49High-confidence mutations were assessed if five or more cell lines in the set of interest had a non-synonymous mutation in that gene, and if five or more cell lines had wild-type copies of that gene. CNAs were assessed if five or more cell lines had an alteration (deletion or amplification) of that gene, and if five or more cell lines lacked that alteration. For each comparison, T-tests were performed between the log(IC50) value of every tested compound. For single-cancer type analyses, a threshold of p < 0.01 was used to identify significance, while for pan-cancer analyses, a threshold of p < 0.0001 was used to identify significance.

50Code is available on GitHub at https://github.com/joan-smith/genomic-features-survival (; copy archived at https://github.com/elifesciences-publications/genomic-features-survival).

51Kaplan-Meier plots were generated using Graphpad Prism. Deletions and amplifications in Kaplan-Meier plots correspond to CNAs > |0.3|; deep deletions and high-copy gains correspond to CNAs > |1|. P values reported in KM plots were generated by the log-rank test in Prism. Note that Kaplan-Meier plots are displayed in this manuscript primarily for the ease of visualizing patient outcomes. Z scores were always generated with Cox proportional hazards modeling, which does not require the selection of artificial cut-offs or thresholds for continuous data.

52The 30 frequently-mutated cancer driver genes were acquired from. NCI-SEER statistics were downloaded from https://seer.cancer.gov. Total tumor aneuploidy scores, ABSOLUTE-determined purity values, and leukocyte infiltration was obtained from. Hyper-mutated samples were obtained from. Lollipop plots were generated using Lollipops software. Density plots were generated with Python scripts using matplotlib (https://matplotlib.org/). Single base-pair mutations were mapped to codons using PolyPhen-2.

53Multiple statistical techniques have been developed to perform survival or ‘time-to-failure’ analysis (reviewed in ). These include Kaplan-Meier analysis, Cox proportional hazards regression, accelerated failure time modeling, and many others. In this paper, we chose to apply Cox proportional hazards regression to analyze cancer survival data. The Cox model is represented by the following function:

54Where t is the survival time, h(t, X) is the hazard function, h0(t) is the baseline hazard, Xi is a potential prognostic variable, and βi indicates the strength of the association between a prognostic variable and survival. In this model, patients have a baseline, time-dependent risk of death [h0(t)], modified by time-independent prognostic features that either increase (βi>0) or decrease (βi<0) risk of death. In this paper, we report Z scores, which are calculated by dividing the regression coefficient (βi) by its standard error.

55Cox proportional hazards modeling was chosen for several reasons. First, unlike Kaplan-Meier analysis, Cox models do not require the selection of a threshold or cut-off, so continuous data like gene expression values do not need to be dichotomized. (Note that in this manuscript, Kaplan-Meier plots are provided for visualization purposes, but the reported Z scores are always from Cox models). Secondly, Cox models can accept both continuous and discrete input data, allowing this approach to be used to analyze both binary (e.g., mutant vs. non-mutant) and continuous (e.g., gene copy number) genomic features. Thirdly, Cox models are amenable to both univariate (i = 1) and multivariate (i > 1) analyses. Fourthly, Cox regression allows us to calculate Z scores and a p value for each association, as Z scores represent the number of standard deviations from the mean of a normal distribution. Fifthly, Z scores encode the directionality of an association: poor prognostic factors will exhibit βi values greater than 0, while favorable prognostic factors will exhibit βi values less than 0. This allows ‘positive’ and ‘negative’ survival features to be directly compared. Sixthly, Z scores are useful for meta-analyses, as they can be combined using Stouffer’s Method:

56Seventhly, Cox proportional hazards modeling is commonly used in both previous genome-wide survival analyses and in numerous clinical biomarkers studies, facilitating comparison with other biomarker discovery efforts.

57To verify the underlying normality of the Z scores, we generated qq plots for gene copy number values (Figure 1—figure supplement 2C). The resulting distributions for CNAs were generally linear, as expected, with occasional shoulders at low and high Z scores. We similarly calculated Z scores for all genes harboring coding-sequence mutations; however, we discovered that this resulted in plateaus around the origin in multiple cancer types. These aberrations were caused by the occurrence of rare, random mutations in multiple genes that lacked any prognostic power. To eliminate these plateaus, we experimented with different thresholds for mutational analysis. Considering only mutations that occurred in a certain percentage of cancer patients diminished the appearance of the plateaus, but high thresholds also eliminated from consideration mutations in a number of known cancer drivers. We selected a 2% threshold to balance between maintaining the normality of the Z score distribution while also retaining infrequent but significant mutations in driver genes.

58Note that in many survival analysis papers, a ‘feature selection’ step is included to identify a minimal number of features that can accurately identify at-risk patients. We performed an unbiased, whole-genome analysis without feature selection, to generate a Z score for every gene and for every feature type in the genome. No feature selection step is applied in this work.

59Patient cohorts that were assembled for the TCGA were collected in order to allow a molecular analysis of the major cancer subtypes found within the United States. Though clinical information was collected for nearly all patients, these cohorts were not specifically chosen in order to conduct survival studies. We posit that our survival analysis is appropriate for several reasons. First, we verified that the overall survival times of patients within the TCGA is highly consistent with national epidemiological data collected by the NCI (Figure 1—figure supplement 2D–E). Secondly, we found that many well-established biomarkers hold prognostic significance in TCGA cohorts, including IDH1 mutations in glioma (Figure 1—figure supplement 7), TP53 mutations in breast cancer (Figure 1—figure supplement 3), tumor stage and grade in multiple cancer types (Figure 2—figure supplement 3), and more. Thirdly, we validated the survival patterns that we describe in the TCGA in several independent patient cohorts, indicating that these are not TCGA-specific phenomena (Figure 4). Fourthly, in an independent analysis of the quality of clinical annotations in the TCGA, none of the cohort/endpoint combinations chosen for this study were classified as ‘not recommended for use.’ Fifthly, our efforts build upon a robust body of work that has also performed survival analyses on TCGA cohorts, and, in some cases, similarly validated findings from the TCGA in independent patient populations. Finally, we note that the TCGA has several benefits over standard investigator-initiated survival studies. Patient samples were collected and analyzed in an unbiased manner, precluding the possibility of the ‘file-drawer problem’ (failing to publish negative results) or post-hoc sample size adjustment (ending patient enrollment when a significant result is found). Significantly more molecular data is available from TCGA tumors than in any other comparably-sized dataset, which allows for multivariate and correlational analyses of different facets of tumor genomes. All data from the TCGA and all code from this manuscript are publicly-available, allowing easy replication and extension upon this analysis.