Pathway and network analysis of more than 2,500 whole cancer genomes

The catalog of cancer driver mutations in protein-coding genes has greatly expanded in the past decade. However, non-coding cancer driver mutations are less well-characterized and only a handful of recurrent non-coding mutations, most notably TERT promoter mutations, have been reported. Motivated by the success of pathway and network analyses in prioritizing rare mutations in protein-coding genes, we performed multi-faceted pathway and network analyses of non-coding mutations across 2,583 whole cancer genomes from 27 tumor types compiled by the ICGC/TCGA PCAWG project. While few non-coding genomic elements were recurrently mutated in this cohort, we identified 93 genes harboring non-coding mutations that cluster into several modules of interacting proteins. Among these are promoter mutations associated with reduced mRNA expression in TP53, TLE4, and TCF4. We found that biological processes had variable proportions of coding and non-coding mutations, with chromatin remodeling and proliferation pathways altered primarily by coding mutations, while developmental pathways, including Wnt and Notch, altered by both coding and non-coding mutations. RNA splicing was primarily targeted by non-coding mutations in this cohort, with samples containing non-coding mutations exhibiting similar gene expression signatures as coding mutations in well-known RNA splicing factors. These analyses contribute a new repertoire of possible cancer genes and mechanisms that are altered by non-coding mutations and offer insights into additional cancer vulnerabilities that can be investigated for potential therapeutic treatments.


Introduction
Over the past decade, cancer genome sequencing efforts such as The Cancer Genome Atlas (TCGA) have identified millions of somatic genetic aberrations; however, the annotation and interpretation of these aberrations remains a major challenge 1 . Specifically, while some aberrations occur frequently in specific cancer types, there is a "long tail" of rare aberrations that are difficult to distinguish from random passenger aberrations in modestly sized patient cohorts 2,3 . In many cancers, a significant proportion of patients do not have known coding driver mutations 4 , suggesting that additional driver mutations remain undiscovered. To date, the vast majority of known driver mutations affect protein-coding regions; only a few non-coding driver mutations, most notably mutations in the TERT promoter [5][6][7] , have been identified. Recent studies from the Pan-Cancer Analysis of Whole Genomes (PCAWG) project of the International Cancer Genome Consortium (ICGC) reveal few recurrent non-coding drivers in analyses of individual genes and regulatory regions 7 .
Cancer driver mutations unlock oncogenic properties of cells by altering the activity of hallmark pathways 8 . Accordingly, cancer genes are known to cluster in small number of cellular pathways and interacting subnetworks 3,9 . Previously, pathway and network analysis has proven useful for implicating infrequently mutated genes as cancer genes based on their pathway membership and physical/regulatory interactions with recurrently mutated genes [10][11][12][13][14] . However, the interactions between coding and non-coding driver mutations have not been systematically explored.
We performed pathway and network analysis of coding and non-coding somatic mutations from 2,583 tumors from 27 tumor types compiled by the Pan-Cancer Analysis of Whole Genomes (PCAWG) project of the International Cancer Genome Consortium (ICGC) 15 , the largest collection of uniformly processed cancer genomes to date. We derive a consensus set of 93 high-confidence pathway-implicated driver genes with non-coding variants (PID-N) and a consensus set of 87 pathway-implicated driver genes with coding variants (PID-C) using seven pathway and network analysis methods. Both sets of PID genes, particularly the PID-N set, contain rarely mutated genes that were not identified by individual recurrence tests but interact with other well-known cancer genes. In total, 121 novel PID-N and PID-C genes are revealed as promising candidates, expanding the landscape of driver mutations in cancer.
Furthermore, we examined the contribution of coding and non-coding mutations in altering biological processes, finding that while chromatin remodeling and some well-known signaling and proliferation pathways are altered primarily by coding mutations, other important cancer pathways, including developmental pathways such as Wnt and Notch pathways, are altered by both coding and non-coding mutations in PID genes. Intriguingly, we find many non-coding mutations in PID-N genes with roles in RNA splicing, and samples with these non-coding mutations exhibit similar gene expression signatures as samples with well-known coding mutations in RNA splicing factors. Our analysis demonstrates that somatic non-coding mutations in untranslated and cis-regulatory regions constitute a complementary set of genetic perturbations with respect to coding mutations, affect several biological pathways and molecular interaction networks, and should be further investigated for their role in the onset and progression of cancer..

Results
The long tail of coding and non-coding cancer mutations highlights opportunities for pathway and network analysis We analyzed the genes targeted by single nucleotide variants (SNVs) and short insertions and deletions (indels) identified by whole genome sequencing in the 2,583 ICGC PCAWG tumor samples from 27 tumor types. Our pathway and network analyses focused on a subset of 2,252 tumors that excluded melanomas and lymphomas due to their atypical distributions of mutations in regulatory regions 16 . We analyzed the pan-cancer driver p -values of single protein-coding and non-coding elements predicted by the PCAWG consensus driver analysis 7 including exons, promoters, untranslated regions (5' UTR and 3' UTR), and enhancers. This PCAWG consensus driver analysis integrates p -values from 16 driver discovery methods, resulting in consensus driver p -values for coding and non-coding elements. Among protein-coding driver p -values of the pan-cancer cohort, 75 genes were highly significant (FDR < 0.1; Supplemental Figure S1 ) and an additional 7 genes were observed at near-significant levels (0.1 ≤ FDR < 0.25). These numbers are consistent with previous reports of a "long tail" of driver genes with few highly-mutated genes and many genes with infrequent mutations across cancer types 2,17 . Non-coding mutations exhibit a similar long-tail distribution with even fewer significant genes (8 genes at FDR < 0.1 and 2 genes at 0.1 ≤ FDR < 0.25). No single gene has both significant or near-significant coding and non-coding driver p -values (FDR < 0.25), suggesting that non-coding mutations target a complementary set of genes as coding mutations.
Earlier studies have demonstrated that proteins harboring coding driver mutations interact with each other in molecular pathways and networks significantly more frequently than expected by chance 2,3,[9][10][11]13 . We observed significant numbers of interactions between both coding and/or non-coding elements with more mutations than expected by chance, suggesting 4 that pathway and network methods may be able to identify rare driver events that are not prioritized by single-element analyses ( Supplemental Figure S2; Coding and non-coding mutations cluster on networks in Supplement ).
Consensus pathway and network analysis reveals possible non-coding driver mutations We performed a comprehensive pathway and network analysis of cancer drivers using the results of the single-element driver discovery study of the PCAWG project 7 19 , and SSA-ME 20 ) that utilized information from molecular pathways or protein interaction networks ( Figure 1 , Methods ). Each method nominated genes, and consensus sets of genes with possible coding and non-coding driver mutations were defined as the genes found by at least four of the seven methods ( Supplemental Tables S1-S4 ). All methods were calibrated on randomized data ( Individual pathway and network algorithms in Supplement ) .
When using non-coding mutations alone, a consensus of pathway and network analysis results on non-coding data identified 62 genes. In contrast, and as one might expect, the coding analysis resulted in substantially more genes, producing a set of 87 pathway-implicated driver genes with coding variants (PID-C). To increase the sensitivity for detecting contributions provided by non-coding mutations, we devised a "non-coding value-added" (NCVA) procedure   ; Figures 2A and 2C, Supplemental Figure S7A ). The PID-C genes have significantly higher coding gene scores than non-PID-C genes (rank sum test p = 1.72 ⨉ 10 -58 ; median rank 48 of PID-C genes), and each of the 87 PID-C genes improves the score of its network neighborhood (19.7 genes expected; p < 10 -6 ; Supplemental Table S5 ). This network neighborhood analysis shows that PID-C genes are not implicated solely by their network neighbors 14 but themselves contribute significantly to their discovery by pathway and network methods. The 87 PID-C genes also include 31 genes that are not statistically significant (FDR >   Tables S7, S11 ).
We found significantly more interactions between PID-C genes that expected by chance using a node-degree preserving permutation test (64 interactions observed vs. 40 interactions expected, p < 10 -6 ), a near significant number of interactions between PID-N genes (18 vs. 12 expected, p = 6.8 ⨉ 10 -2 ), and significantly more interactions between both PID-C and PID-N genes (67 vs. 40 expected, p = 6 ⨉ 10 -4 ), demonstrating an interplay between coding and non-coding mutations on physical protein-protein interaction networks ( Network annotation in Methods ).
Overall, we organized the interactions between PID-C and PID-N genes into five biological processes: core drivers, chromatin organization, cell proliferation, development, and RNA splicing ( Figure 4A ). While the high frequency of molecular interactions between PID-C and PID-N genes is expected since such interactions were used as a signal in pathway and network methods, the specific structure of these interactions illustrates the relative contributions of coding and non-coding mutations in individual subnetworks.
We further characterized the molecular pathways enriched among our PID-C and PID-N using the g:Profiler web server 29 ( Figure 4B , Supplemental Figure S9, Supplemental Tables S15-S18, Pathway annotation in Methods ). Since our methods use pathway databases and interaction networks as prior knowledge, enrichment with known pathways is expected.
However, the enrichment results provide clues about the modular organization of the pathways and allow us to assess the relative contributions of coding and non-coding mutations in each pathway. Overall, 63 pathways were enriched for PID-C genes and 13 pathways were enriched for PID-N genes (FDR < 10 -6 ).
We further grouped these pathways into 29 modules using overlaps between annotated pathways recorded in the pathway enrichment map ( Supplemental Figure S11 ).For each enriched module, we examined whether PID-C, PID-N, or both types of genes were responsible for the observed enrichment. This produced a clustering of modules and PID genes into four biological processes: chromatin organization, cell proliferation, development and RNA splicing ( Figure 4B ).
We found that pathways in the chromatin and cell proliferation processes -including chromatin remodeling and organization, histone modification, apoptotic signaling, signal transduction, Ras signaling, and cell growth -were altered primarily by coding mutations in PID-C genes. This is not surprising as these pathways contain many well-known cancer genes, such as TP53 , KRAS , BRAF , cyclin dependent kinase inhibitors, EGFR , PTEN , and RB1 .
Several signaling pathways contain significant numbers of both PID-C and PID-N genes, indicating that non-coding mutations provide additional avenues for disrupting key molecular interactions. These pathways include the Wnt signaling pathway (FDR = 6.8 ⨉ 10 -13 ), which was predominantly targeted by coding mutations but was also targeted by non-coding mutations in and RELN and the PID-C genes ATM and SMAD4 . In these cases, non-coding mutations complement coding mutations that disrupt these pathways, covering significant numbers of additional patients.
Intriguingly, we find that RNA splicing pathways were affected primarily by non-coding mutations (FDR = 7.6 ⨉ 10 -9 ). A total of 17 PID-N genes belonged to splicing-related pathways ( Supplemental Figure S12C ), including several heterogeneous nuclear ribonucleoproteins (hnNRP) and serine and arginine rich splicing factors (SRSFs). None of these PID-N genes were significantly mutated according to single-element tests of the PCAWG driver discovery analysis. We did not find any significant (FDR < 0.1) in cis associations between non-coding mutations and altered expression of these genes. Thus, we explored potential in trans effects on pathway expression changes. We found that non-coding mutations in splicing-related PID-N genes largely recapitulate a recently published association by TCGA 32 between coding mutations in several splicing factors and differential expression of 47 pathways ( Figure 5 ).
Specifically, we identified three clusters of mutations (C1, C2, and C3 in Figure 5A and In addition to the above modules, we also found that transcription factors were well represented among both the PID-C and PID-N genes. In total, 9 PID-C genes are transcription

Discussion
While single-region tests in the PCAWG project identified only a few non-coding driver elements, our integrative pathway and network analysis further expands the list of genes with possible non-coding driver mutations, extending into the "long tail" of rare mutations. In particular, we find that genes with either coding or non-coding mutations are linked in pathways and networks, and that pathway databases and interaction networks can be leveraged as prior knowledge to identify additional possible non-coding drivers that are too infrequently mutated to be detected by single-element tests. In total, our integrative pathway analysis identified 87 pathway-implicated driver genes with coding variants (PID-C) and 93 pathway-implicated driver genes with non-coding variants (PID-N). Importantly, 90 PID-N genes were not statistically significant (FDR > 0.1) by single-element tests on non-coding mutation data, and these genes are key candidates for future experimental characterization. Among them, we find that promoter mutations in TP53 , TLE4 , and TCF4 are associated with reduced expression of these genes.
We find that coding and non-coding driver mutations largely target different genes, and contribute differentially to pathways and networks perturbed in cancer. While some cancer pathways are targeted by both coding and non-coding mutations, such as the Wnt and Notch signaling pathways, other pathways appear to be predominately altered by one class of mutations. In particular, we find non-coding mutations in multiple genes in the RNA splicing pathway, and samples with these mutations exhibit gene expression signatures that are concordant with gene expression changes observed in samples with coding mutations splicing factors SF3B1 , FUBP1 , and RBM10 32 . Together these results demonstrate that rare non-coding mutations may result in similar perturbations to both common and complementary biological processes.
There are several caveats to the results reported in this study. First, there is relatively low power to detect non-coding mutations in the cohort, particularly in cancer types with small numbers of patients. Second, transcriptomic data was available for only a subset of samples, further reducing our ability to validate our predictions using gene expression data. Third, our pathway and network analysis relied on the driver p -values from the PCAWG consensus driver analysis 7 . This analysis accounts for regional variations in the background mutation rate across the genome. However, if these corrections are inadequate and the uncorrected confounding variables are correlated with gene membership in pathways and subnetworks, then the false positive rates in our analysis may be higher than estimated. All of these factors, plus other unknown confounding variables, make it difficult to assess the false discovery rate of our predictions, particularly for PID-N genes. Further experimental validation of these predictions is necessary to determine the true positives from false positives in our PID gene lists.
While pathway and network analysis was successful in revealing potential new cancer-associated genes impacted by non-coding mutations, future investigations that consider the changing landscape of gene regulation and pathway interactions across tissues may offer a new perspective on the data. Specifically, each cell type has a different epigenetic wiring and regulatory machinery, and non-coding mutations may target cell type-specific vulnerabilities.
Approaches that incorporate tissue-specific gene-gene regulatory logic may be successful in revealing new classes of drivers unexplored with our current approaches.
In conclusion, our pathway-and network-driven strategies enable us to interpret the coding and non-coding landscape of tumor genomes to discover driver mechanisms in interconnected systems of genes. This approach has multiple benefits. First, by broadening our mutation analysis from single genomic elements to pathways and networks of multiple genes, we identify new components of known cancer pathways that are recurrently altered by both coding and non-coding mutations, and thus likely to be important in cancer. Second, we identify new pathways and subnetworks that would remain unseen in an analysis focusing on coding sequences. Investigation of the coding and non-coding mutations that perturb these pathways 13 and networks will enable more accurate patient-stratification strategies, pathway-focused biomarkers, and therapeutic approaches.                             Online Methods

Mutation and Pathway Data
We combined several pathways and interaction networks with gene scores derived from the PCAWG drivers analysis 1 for use by pathway and networks methods. Here, we use the term "pathway methods" for those approaches that make use of sets of related genes for their analysis while the term "network methods" are reserved for those that also incorporate the interactions among the genes and/or their products.
In cases where the PCAWG drivers analysis reported multiple p -values for the same genomic element, we used the smallest reported p -value for that element.

Derivation of gene scores
Pathway databases and gene interaction networks typically record information at the level of individual genes. Thus, we formed coding and non-coding gene scores by combining PCAWG driver p -values across coding and/or non-coding (core promoter, 5' UTR, 3' UTR, enhancer) genomic elements as follows. Let p x ( g ) be the driver p -value for element x of gene g from the PCAWG drivers analysis 1 . We combined p -values from multiple elements using Fisher's method, where we selected the minimum p -value min( p promoter ( g ), p 5'UTR ( g )) for overlapping core promoter and 5' UTR elements on gene g and the minimum p -value p enhancer ( g ) of all enhancers targeting gene g . Using this approach, we defined the following gene scores on 1 coding (GS-C), non-coding, (GS-N), and combined coding and non-coding (GS-CN) genomic elements: 1. GS-C: p C ( g ) = p coding ( g ) 2. GS-N: p N ( g ) = fisher(min( p promoter ( g ), p 5'UTR ( g )), p 3'UTR ( g ), p enhancer ( g )) 3. GS-CN: p CN ( g ) = fisher( p coding (g), min( p promoter ( g ), p 5'UTR ( g )), p 3'UTR ( g ), p enhancer ( g )) Here, p = fisher( p 1 , …, p k ) is Fisher's method, i.e., -2 ∑ k i=1 ln( p i )~ 2 2 k , for independently and identically distributed p 1 , …, p k~ U (0, 1), where 2k is the degrees of freedom in the calculation.
Moreover, when the driver p -value for a genomic element was undefined, we omitted that element from the calculation and reduced the number of degrees of freedom.
For the pathway and networks methods that analyze individual mutations, we used mutations from the PCAWG MAF ( syn7118450 ) on the same genomic elements ( syn5259890 ) as the PCAWG drivers analysis, i.e., coding, core promoter, 5' UTR, 3' UTR, and enhancer. We removed melanoma and lymphoma samples as well as 69 hypermutated samples with over 30 mutations/MB ( syn7222520 , syn7814911 ). We also removed mutations in elements that the PCAWG drivers analysis had manually examined and discarded (see above), resulting in lists of mutations used for later assessing biological relevance of our results ( syn8103141 , syn9684700 ).
Network methods used interactions from three interaction networks: the largest connected subnetwork of the ReactomeFI 2015 interaction network 9 ( syn3254781 ) with high-confidence (≥ 0.75 confidence score) interactions, which we treated as undirected; the largest connected subnetwork of the iRefIndex14 interaction network 10 , which we augmented with interactions from the KEGG pathway database 6 ( syn10903761 ); and the largest connected subnetwork of the STRING v10 network 11 ( syn11712027 ) with high-confidence (> 9 confidence score) interactions. The BioGRID interaction network 12   Non-coding value-added (NCVA) procedure The GS-CN results leverage both coding and non-coding mutation data, improving the detection of weaker pathway and network signals. We devised a non-coding value-added (NCVA) procedure to separate the coding and non-coding signals in this combined analysis, resulting in a set of NCVA genes for which the non-coding mutation data makes a statistically significant contribution to their discovery in the GS-CN results. Specifically, we evaluated the statistical significance of genes in the GS-CN results using a permutation test where the driver p -values for coding elements were fixed and the driver p -values for non-coding elements were permuted. This procedure identified the subset of the GS-CN results that were reported infrequently ( p < 0.1) on permuted data and thus more likely to be true positives. Each method's NCVA results were added to that method's overall set of non-coding results (GS-N).

Consensus results for pathway and network methods
We defined a consensus set of genes for each set of results: GS-C results, GS-N results, GS-CN results, and GS-N combined with NCVA results, across our seven pathway and network methods. Specifically, we defined a gene to be a consensus gene if it was found by a majority (≥ 4/7) of the pathway and network methods. For our analysis, we focused on the consensus GS-C results, which we call the pathway-implicated driver genes with coding variants (PID-C), and the consensus from the GS-N results combined with NCVA results, which we call the pathway-implicated driver genes with non-coding variants (PID-N). We defined PID-C genes as the 87 genes in the consensus of the GS-C results, and we defined PID-N genes as the 93 genes in the consensus of each method's GS-N results combined with its NCVA results. 4 Downstream Interpretation of Pathway-Implicated

Drivers
We performed several analyses to assess the biological relevance of PID-C and PID-N genes.

Identification of mutational signatures of PID genes
We performed a permutation-based enrichment test for mutation signatures from PCAWG mutation signatures analysis 18 . We identified the most likely mutation signature for each non-coding mutation in PID-N genes and compared them to randomly chosen non-coding mutations in non-PID-N genes.
Gene scores improve network neighborhood scores of PID genes To assess the extent to which gene scores on PID genes contribute to their detection by pathway and network methods, we considered the contribution of each PID gene's score to the score of its network neighborhood in the BioGRID interaction network.
For each PID gene g , we used Fisher's method to combine the gene scores of the first-order network neighbors of g both with and without the score of g itself. In particular, for gene g , let p ( g ) be the gene score for g and N ( g ) be the network neighborhood of g . Then is a score for the network neighborhood of g when including gene g and is a score for the network neighborhood of g when excluding gene g .

5
If the network neighborhood of g has a smaller p -value with g than without g , i.e., p N(g) with < p N(g) without , then gene g improves the score of the network neighborhood, suggesting that the gene score of g plays a role in its detection by pathway and network methods. Alternatively, if the network neighborhood of g has a larger p -value with g than without g , i.e., p N(g) with > p N(g) without , then gene g worsens the score of the network neighborhood, suggesting that the gene scores of the network neighbors of g are predominantly responsible for the detection of g by pathway and network methods.
We performed this test for every PID-C gene with GS-C gene scores and every PID-N gene with GS-N gene scores. We also sampled genes uniformly at random from the network (87 for PID-C genes and 93 for PID-N genes; 10 6 trials) to ascertain whether significantly more PID genes that improved the scores of their network neighborhoods than expected by chance.

Expression analysis of PID genes
We evaluated whether mutation status of each PID gene was correlated with RNA expression. We used PCAWG-3 gene expression data ( syn5553991 ), which was averaged from TopHat2 and STAR-based alignments, with FPKM-UQ normalization. Tumor type and copy-number aberrations are known to be covariates for gene expression, so we conditioned on tumor types and annotated copy-number aberrations.
We used the following procedure to assess expression correlations on individual tumor types. We only considered cases with at least 3 mutated samples and 3 non-mutated samples to restrict our analysis to cases with sufficient statistical power. For each PID-C gene or each non-coding element in a PID-N gene, we partitioned the samples with expression data into a set A of samples with mutation(s) in the element and a set B of samples without mutations in the element. We performed the Wilcoxon rank-sum test for the expression of the gene in sets A and B and performed the Benjamini-Hochberg correction on each coding or non-coding element to provide FDRs.
We used the following procedure to assess expression correlations across tumor types.
We only considered cases with at least 1 mutated sample and 1 non-mutated sample to restrict our analysis to cases with sufficient statistical power. For each PID-C gene and each

Network annotation of PID genes
We performed a permutation test to evaluate the statistical significance of the number of interactions in the BioGRID high-confidence interaction network between PID-C genes, the number of interactions between PID-N genes, and the number of interactions between PID-C and PID-N genes, i.e., when a PID-C gene interacts with a PID-N gene. To compute the permutation p -value we sampled random networks uniformly at random from the collection of networks with the same degree sequence as the BioGRID network.
We found connected subnetworks of 46 PID-C genes (31 genes expected, p = 9 ⨉ 10 -4 ) and 16 PID-N genes (10 genes expected, p = 6.1 ⨉ 10 -2 ) in the high-confidence BioGRID 19 protein-protein interaction (PPI) network. The union of the PID-C and PID-N genes formed a larger connected subnetwork of 73 genes ( Figure 4A ). These connected subnetworks were significantly larger than expected by chance according to this permutation test (57 genes expected, p = 2.2 ⨉ 10 -3 ). Further, we observed statistically significant numbers of protein-protein interactions between PID-C and PID-N genes (67 interactions observed vs. 45 expected, p = 6 ⨉ 10 -4 ), suggesting that the associated mutations may target an overlapping set of underlying pathways. The PID-C genes were connected by significantly more interactions than expected (64 vs. 40 expected, p < 10 -4 ) and the PID-N genes were interconnected at a sub-significant level (18 vs 12 expected, p = 6.8 ⨉ 10 -2 ). Thus certain pathways are affected by either coding or non-coding mutations, but some pathways are affected by a complement of both coding and non-coding mutations.

7
Pathway annotation of PID genes Using g:Profiler 20 , we performed a pathway enrichment analysis for PID genes and 12,061 gene sets representing GO biological processes and Reactome pathways. We used the Benjamini-Hochberg correction to control the FDR of the results.
Characterization of PID genes in RNA splicing Separately, we computed a 2D projection of NES scores using t-Distributed Stochastic Neighbor Embedding (tSNE).

Additional Information
See Supplement for more information about data processing and details of individual network and pathway methods. 8