Novel Driver Strength Index highlights important cancer genes in TCGA PanCanAtlas patients

Introduction

According to the mainstream view, cancer is initiated and promoted by mutations in so-called “driver” genes, whereas mutations in “passenger” genes bear no effect and simply “travel” along, during the course of somatic evolution (Vogelstein et al., 2013). There are multiple approaches to identifying cancer driver genes and to their ranking (Marx, 2014; Raphael et al., 2014). One of the earliest and most intuitive approaches is to use mutation frequency, typically corrected by background mutation rate in a gene, e.g. by the rate of synonymous mutations. Some notable examples of such algorithms are MuSiC (Dees et al., 2012), MutSigCV (Lawrence et al., 2013) and dNdScv (Martincorena et al., 2017). Interestingly, even after background correction, high recurrence of a mutated gene in a cancer cohort does not necessarily translate to the statement that this gene is a strong driver. We can imagine a gene that is mutated in the majority of cancer patients (e.g., because it has multiple suitable sites for a driver mutation, such as catalytic sites, sites of post-translational modifications or protein-protein interfaces) but has a weak contribution to cancer progression in each of these patients (e.g., because this gene is redundant). We can also imagine a gene that is mutated rarely (e.g., because it has only one suitable site for a driver mutation) but if the mutation does occur it immediately leads to cancer (e.g., because this gene is in a key position to control cell growth). The former would be an example of a frequent but weak driver, whereas the latter would be a rare but strong driver. Overall, algorithms based on mutation recurrence cannot reliably determine driver gene strength.

There is a large group of algorithms that aim to predict and rank driver genes according to the impact of mutations on protein structure and activity. Some notable examples include HotMAPS (Tokheim et al., 2016a), 2020plus (Tokheim et al., 2016b), OncodriveFML (Mularoni et al., 2016), OncodriveCLUSTL (Arnedo-Pac et al., 2019) and CHASMplus (Tokheim & Karchin, 2019). Typically, these algorithms are based on mutation clustering, sequence evolutionary conservation, protein domains and three-dimensional protein structure, or combinations thereof, sometimes with the use of machine learning. These methods can reliably identify if a given protein’s structure and function are disturbed by a given mutation and to what degree. However, they are much less suitable for identification of the role of this protein in the cellular and microenvironmental context, which is critical to determine whether it will drive cancer and how strongly. Thus, the strength of driver genes is unlikely to be determined by structural approaches.

Some algorithms aim to predict and rank drivers based on data from protein interaction networks, e.g., NetBox (Cerami et al., 2010), DriverNet (Bashashati et al., 2012), NetSig (Horn et al., 2018) and Hierarchical HotNet (Reyna, Leiserson & Raphael, 2018). The idea is that a gene having multiple connections with other genes, i.e. playing the role of a network hub, will have more dramatic influence on the cell in case of mutation (Bowler, Wang & Ewing, 2015). However, whilst mutations in network hubs are indeed likely to cause severe disturbance in the cell, this would rather lead to cell death than to oncogenic transformation. Thus, network-based approaches are also not likely to uncover the actual strength of driver genes.

We have recently quantified the number of driver events in individual patients in various cancer types and discovered very high variability in this number even within the same cancer type (Vyatkin et al., 2022). We therefore asked a question: what is the reason that for some patients one driver event is sufficient to drive cancer, whereas the others do not develop cancer until dozens of driver events accumulate? We hypothesized that the major reason is driver strength–one strong driver is equivalent in its cancer-promoting potency to several weak drivers. Thus, we define the strength of a cancer driver as its share of the contribution to the initiation and progression of a given cancer relative to other drivers. It immediately follows that a driver that is able to singlehandedly initiate cancer is a strong driver, whereas a driver that is present in cancers only amongst the multitude of other drivers is very likely a weak driver.

We reason that a few strong drivers are sufficient to initiate cancer, and there would be no need to accumulate additional drivers. On the other hand, weak drivers would need to accumulate in much higher quantity, until their combined strength would become sufficient to initiate cancer. Therefore, it should be statistically more likely to find strong drivers in patients that have only few driver mutations in their tumours, and less likely to find them in patients with multiple drivers per tumour. Likewise, it should be statistically less likely to find weak drivers in patients that have only few driver mutations in their tumours, and more likely to find them in patients with multiple drivers per tumour. Hence, we propose the Driver Strength Index (DSI) that takes into account the frequencies of mutation of a given driver gene in groups of patients with different total number of driver mutations, and gives priority weights to groups with fewer mutations. We also propose a modification of this index that is completely independent of the overall frequency of mutation of a given driver gene–the Normalized Driver Strength Index (NDSI).

Calculating these indices requires data on the number of driver mutations in each individual patient. The majority of existing driver prediction algorithms work at the cohort level, i.e. they predict driver genes for large groups of patients, usually having a particular cancer type. This does not allow to look at the composition of driver mutations in individual patients. We wrote specific scripts to convert cohort-level predictions into patient-level events (Vyatkin et al., 2022), which also allowed seamless integration of the results from various third-party algorithms, including 2020plus, CHASMplus, CompositeDriver, dNdScv, DriverNet, HotMAPS, OncodriveCLUSTL and OncodriveFML. This is useful, as each individual driver prediction algorithm is based on the unique combination of theoretical concepts and computational methodology and thus has its own strengths and shortcomings, and combining results from multiple algorithms allows to obtain more complete and balanced picture, ensuring that less driver mutations have been missed. We also used a consensus driver gene list from 26 algorithms (Bailey et al., 2018), applied separately to each cancer cohort of TCGA PanCanAtlas, and a list of COSMIC Cancer Gene Census Tier 1 genes affected by somatic SNAs and CNAs (CGC), as it represents the current gold standard of verified cancer drivers (Sondka et al., 2018). To minimize false positives, we used only driver gene-cohort pairs that were independently predicted by at least two of our sources (Vyatkin et al., 2022).

Methods

Our detailed workflow is described in Methods S1. Most of the methodology, except for the PALDRIC modification, enabling the calculation of driver strength indices, and pathway and network analysis of top-ranked driver genes, has been previously published in (Vyatkin et al., 2022). Here, we describe it briefly.

Results

We calculated the number of various types of driver events in individual genes, chromosome arms or full chromosomes for each cancer type, tumour stage, age group, as well as for patients with each total number of driver events from 1 to 100. We performed the analyses for the total population and for males and females separately, and, for each group, plotted the histograms of top 10 driver events in each class and overall (for data and graphs see Data S1). In Fig. 1 we present the overall ranking of genes for all TCGA PanCanAtlas cohorts combined. It can be seen that PIK3CA is the oncogene with the highest number of SNAs, as well as the highest number of simultaneous occurrences of SNAs and gene amplifications. MYC is the oncogene with the highest number of amplifications. TP53 is the tumour suppressor with the highest number of SNAs, as well as the highest number of instances of simultaneous occurrences of an SNA in one allele and a deletion of the other allele. It is also the top mutated gene when driver events of all classes are counted. CDKN2A and PTEN are tumour suppressors with the highest number of deletions. Losses of chromosomes 13 and 22 are the most frequent cancer-promoting chromosome losses, whereas gains of chromosomes 7 and 20 are the most frequent cancer-promoting chromosome gains. Losses of 8p and 17p arms are the most frequent cancer-promoting chromosome arm losses, whereas gains of 1q and 8q arms are the most frequent cancer-promoting chromosome arm gains. Overall, these results are expected and indicate that our analytic pipelines work as they should.

Next, we calculated the Driver Strength Index (DSI) (see Eq. (2) in Methods). Surprisingly, we do not see much change compared to the simple frequency-of-mutation approach (Fig. 2). The only dramatic difference is that BRAF became the top SNA-based oncogene according to DSI, whereas PIK3CA dropped to the second place, lagging behind by a wide margin. Also, PIK3CA overtook MYC as the top CNA-based oncogene, and PTEN displaced CDKN2A from the top CNA-based tumour suppressor spot. Moreover, members of several gene families appeared in the top 10 lists, such as KRAS, NRAS and HRAS in the SNA-based oncogenic events list, histones HIST2H2BE and HIST1H3B in the CNA-based oncogenic events list, or lysine methyltransferases KMT2C and KMT2D in the SNA-based tumour suppressor events list. This indicates that our approach is indeed meaningfully selecting for some biological attributes, which are not selected by simple frequency sorting. Finally, multiple small changes of ranking positions happened, nevertheless not affecting the overall picture. We think the reason for the limited effect of changes is that DSI is still very much affected by the overall frequency of gene mutation, due to ${\displaystyle \frac{p_{Ai}}{p_i}}$ component. Hence, to uncover the true driver strength unrelated to the mutation frequency, further normalization is required.

Therefore, we propose the Normalized Driver Strength Index (NDSI) that corrects for the effects of mutation frequencies (see Eq. (3) in Methods). As can be seen in Fig. 3, this time the rankings are completely different from both DSI and frequency-based approaches. GTF2I conquers the top spot amongst SNA-based oncogenes and overall, SPOP becomes number one CNA-based oncogene, and MET occupies the first line of mixed oncogene rating. ATRX, CSDE1 and NF2 become the top SNA-based, CNA-based and mixed tumour suppressors, respectively. NDSI reveals the losses of chromosomes 12 and 3 as the strongest cancer-promoting chromosome losses, whereas the gain of chromosome 17 as the strongest cancer-promoting chromosome gain. NDSI shows that the losses of 19q and 12p arms are the strongest cancer-promoting chromosome arm losses, whereas the gain of 5q arm is the strongest cancer-promoting chromosome arm gain.

Like DSI, NDSI is able to select for specific gene families. Two members of the guanine nucleotide-binding protein family, GNAQ and GNA11, appeared on the top 10 SNA-based oncogenic events and top 10 driver events of all classes lists (Fig. 3). Additionally, one more G protein, GNAS, appeared on the top 50 NDSI-ranked driver list (Table 2). Of note, no members of this family are present on the top 50 DSI-ranked driver list (Table 2). Two members of the isocitrate dehydrogenase family, IDH1 and IDH2, appeared on the top 10 SNA-based oncogenic events list, whereas fibroblast growth factor receptors FGFR2 and FGFR3 appeared on the top 10 mixed oncogenic events list (Fig. 3). The ability of NDSI to prioritize members of specific protein families suggests that this index has an actual biological meaning.

Next, we wanted to analyse top DSI- and NDSI-ranked genes using several common gene list analysis tools. To this aim, we combined the lists of drivers from various classes. If the same gene was affected by more than one kind of alteration, we chose the alteration type with the highest DSI or NDSI, depending on the analysis. Also, we removed the data on chromosome arms and full chromosomes, as external pathway and network analysis tools can work only with genes. Then, we selected top 50 DSI- and NDSI-ranked genes. The resulting lists can be seen in Table 2.

First, we analysed the top 50 DSI- and NDSI-ranked genes for overrepresentation in various Reactome pathways. It can be seen in Fig. 4A that top 50 DSI-ranked genes are significantly overrepresented in such categories as signalling by NOTCH, signalling by PTK6, ESR-mediated signalling, PIP3 activates AKT signalling, signalling by receptor tyrosine kinases, signalling by WNT, signalling by erythropoietin, RAF/MAP kinase cascade, signalling by TGF-beta receptor complex, mitotic cell cycle, meiosis, cell cycle checkpoints, DNA double-stand break repair, generic transcription pathway, epigenetic regulation of gene expression, RNA polymerase I transcription, circadian clock, chromatin modifying enzymes, diseases of signal transduction by growth factor receptors and second messengers, diseases of cellular senescence, diseases of programmed cell death, cellular responses to stress, activation of HOX genes during differentiation, and transcriptional regulation of granulopoiesis. Top 50 NDSI-ranked genes are significantly overrepresented in even fewer categories (Fig. 4B)–signalling by PTK6, extra-nuclear oestrogen signalling, negative regulation of PI3K/AKT signalling, signalling by receptor tyrosine kinases, GPCR downstream signalling, erythropoietin activates RAS, RAF/MAP kinase cascade, cytokine signalling in immune system, adaptive immune system, haemostasis, generic transcription pathway, chromatin modifying enzymes, and diseases of signal transduction by growth factor receptors and second messengers. Surprisingly, several large categories often deemed important for cancer – cell cycle, DNA replication, DNA repair, autophagy, cellular responses to stress, programmed cell death, cell-cell communication and metabolism – are not affected.

Next, we mapped the top 50 DSI- and NDSI-ranked genes to KEGG cancer pathways. Figure 5A and Table 2 together suggest that top DSI-ranked genes comprise the EGFR/ERBB2/FGFR3-KRAS/NRAS/HRAS-BRAF-MYC pathway, PIK3CA-PTEN pathway, CTNNB1-MYC pathway, TP53-CDKN2A-RB1 pathway and MYC-CUL1-RB1 pathway. Figure 5B and Table 2 together suggest that top NDSI-ranked genes comprise the EGFR/FGFR2/GNAQ/GNA11-NRAS/HRAS/KRAS-BRAF pathway, AKT1-MTOR pathway, and TCEB1-VHL-HIF1A pathway.

Finally, we explored protein-protein interactions in the top 50 DSI- and NDSI-ranked genes using the BioGRID network. Figure 6A shows that although CTNNB1 and EGFR are the biggest hubs of the top-DSI-ranked gene network, their DSI values are much lower than those of BRAF and PTEN, which have fewer connections. Notably, TP53 exhibited the highest DSI value and second-highest connectedness. Similarly, Fig. 6B shows that although EGFR, AKT1 and HRAS are the biggest, centrally located hubs of the top-NDSI-ranked gene network, their NDSI values are much lower than those of GTF2I, BRAF and ZMYM3, located on the periphery of the network. Moreover, the top-NDSI-ranked gene network has much fewer edges than the top-DSI-ranked gene network, despite containing presumably stronger drivers. All of this supports our initially proposed notion that network centrality does not equal driver strength.

Discussion

NDSI places GTF2I on the top spot both amongst the strongest SNA-based oncogenes and amongst the strongest drivers averaged across all classes. The GTF2I-encoded protein binds to the initiator element (Inr) and E-box element in promoters and functions as a regulator of transcription. The GTF2I c.74146970 T > A mutation was detected in 82% of type A and 74% of type AB thymomas (Petrini et al., 2014). GTF2I β and δ isoforms are expressed in thymomas, and both mutant isoforms are able to stimulate cell proliferation in vitro (Petrini et al., 2014). Recently, it has been shown that expression of mutant GTF2I alters the transcriptome of normal thymic epithelial cells and upregulates several oncogenic genes (Kim et al., 2020). GTF2I L424H knockin cells exhibit cell transformation, aneuploidy, and increased tumour growth and survival under glucose deprivation or DNA damage (Kim et al., 2020). Our analysis also shows frequent mutations of GTF2I in the TCGA THYM (thymoma) cohort. GTF2I has been recently named gene of the month and its role in cancer reviewed (Nathany, Tripathi & Mehta, 2021).

SPOP is categorized by NDSI as the strongest CNA-based oncogene. SPOP encodes a protein that is a component of a cullin-RING-based E3 ubiquitin-protein ligase complex that mediates the ubiquitination of target proteins, leading most often to their proteasomal degradation. SPOP is the most commonly mutated gene in primary prostate cancer (Barbieri et al., 2012). SPOP mutations in prostate cancer result in impaired homology-directed repair of double strand breaks and are associated with genomic instability (Boysen et al., 2015). As most cancer-associated mutations in SPOP are missense and almost none are frameshift or nonsense, PALDRIC classifies it as an oncogene. However, SPOP is usually viewed as a tumour suppressor (Clark & Burleson, 2020). Recently it has been discussed that SPOP actually has a dual role, and while being a tumour suppressor in prostate cancer it performs as an oncogene in kidney cancer (Wang et al., 2020). Indeed, cytoplasmic accumulation of SPOP leads to the ubiquitination and degradation of multiple regulators of cellular proliferation and apoptosis, including the tumour suppressor PTEN, ERK phosphatases, the proapoptotic molecule DAXX, and the Hedgehog pathway transcription factor GLI2, and is sufficient to induce tumorigenesis in clear cell renal cell carcinoma (Li et al., 2014). Our analysis shows frequent mutations and amplifications of SPOP in TCGA PRAD (prostate adenocarcinoma) and UCEC (uterine corpus endometrial carcinoma) cohorts.

MET is the strongest mixed (SNA+CNA) oncogene and the second-strongest CNA-based oncogene, according to the NDSI rating. MET encodes a receptor tyrosine kinase that transduces signals from the extracellular space into the cytoplasm by binding to the hepatocyte growth factor ligand. MET regulates many physiological processes including proliferation, morphogenesis and survival. Ligand binding at the cell surface induces dimerization and autophosphorylation of MET on its intracellular domain that provides docking sites for downstream signalling molecules. Following activation by its ligand, MET interacts with the PI3-kinase subunit PIK3R1, PLCG1, SRC, GRB2, STAT3 or the adapter GAB1. Recruitment of these downstream effectors by MET leads to the activation of several signalling cascades including the RAS-ERK, PI3K-AKT, or PLCG-PKC. Mutations in MET are associated with papillary renal cell carcinoma, hepatocellular carcinoma, and various head and neck cancers. Amplification and overexpression of this gene are also associated with multiple human cancers (Comoglio, Trusolino & Boccaccio, 2018; Recondo et al., 2020). Our analysis shows frequent mutations and amplifications of MET in TCGA KIRP (kidney renal papillary cell carcinoma) and LUAD (lung adenocarcinoma) cohorts.

ATRX is ranked by NDSI as the strongest SNA-based tumour suppressor, 9^th strongest CNA-based tumour suppressor, 5^th strongest mixed (SNA+CNA) tumour suppressor and 9^th strongest driver averaged across all classes. ATRX (Alpha-Thalassemia/Mental Retardation Syndrome, X-Linked) encodes a protein that contains an ATPase/helicase domain, and thus it belongs to the SWI/SNF family of chromatin remodelling proteins. ATRX together with DAXX encode a complex that deposits the histone variant H3.3 into repetitive heterochromatin, including retrotransposons, pericentric heterochromatin, and telomeres, the latter of which show deregulation in ATRX/DAXX-mutant tumours (Heaphy et al., 2011; Dyer et al., 2017). ATRX loss induces extensive changes in chromatin accessibility in both repetitive DNA regions and non-repetitive regulatory regions (Liang et al., 2020). These changes are highly correlated with changes in transcription, which lead to alterations in cancer-related signalling pathways, such as upregulation of the TGF-β pathway and downregulation of the cadherin family of proteins (Liang et al., 2020). Our analysis shows frequent mutations and deletions of ATRX in TCGA ACC, GBM, LGG and SARC cohorts.

CSDE1 is ranked by NDSI as the strongest CNA-based tumour suppressor and 5^th strongest driver averaged across all classes. CSDE1 encodes for an RNA-binding protein involved in translationally coupled mRNA turnover. CSDE1 not only promotes and represses the translation of RNAs but also increases and decreases the abundance of RNAs (Guo et al., 2020). CSDE1 loss-of-function mutations and deletions define a Wnt-altered subtype of pheochromocytomas and paragangliomas (Fishbein et al., 2017). Our analysis also shows frequent deletions of CSDE1 in the TCGA PCPG (pheochromocytoma and paraganglioma) cohort.

NF2 is ranked by NDSI as the strongest mixed tumour suppressor and 3^rd strongest CNA-based tumour suppressor. NF2 encodes Merlin (Moesin-ezrin-radixin-like protein), also known as Neurofibromin 2 and Schwannomin. Merlin is a tumour suppressor classically known for its ability to induce contact-dependent growth inhibition (Mota & Shevde, 2020). Loss-of-function mutations or deletions in NF2 cause neurofibromatosis type 2, a multiple tumour forming disease of the nervous system, characterized by the development of bilateral schwannomas, as well as meningiomas and ependymomas (Petrilli & Fernández-Valle, 2016). NF2 is also mutated and deleted in mesotheliomas (Sekido et al., 1995), clear cell renal cell carcinomas (Dalgliesh et al., 2010), collecting duct carcinomas of the kidney (Pal et al., 2016), and renal cell carcinomas with sarcomatoid dedifferentiation (Malouf et al., 2016). Our analysis shows frequent mutations and deletions of NF2 in TCGA KIRC (kidney renal clear cell carcinoma), KIRP (kidney renal papillary cell carcinoma) and MESO (mesothelioma) cohorts.

Interestingly, NDSI prioritized three members of the guanine nucleotide-binding protein (G protein) family: GNAQ, GNA11, and GNAS. Guanine nucleotide-binding proteins function as transducers downstream of G protein-coupled receptors (GPCRs) in numerous signalling cascades. The alpha chain contains the guanine nucleotide binding site and alternates between an active, GTP-bound state and an inactive, GDP-bound state. Signalling by an activated GPCR promotes GDP release and GTP binding. The alpha subunit has a low GTPase activity that converts bound GTP to GDP, thereby terminating the signal. The GNAQ-encoded protein, an α subunit in the Gq class, couples a seven-transmembrane-domain receptor to activation of PLC β. Some GNAQ cancer mutants display normal basal activity and GPCR-mediated activation, but deactivate slowly due to GTPase activating protein (GAP) insensitivity (Garcia-Marcos, Maziarz & Leyme, 2018). GNAQ mutations occur in about half of uveal melanomas, representing the most common known oncogenic mutation in this cancer (Onken et al., 2008). The presence of this mutation in tumours at all stages of malignant progression suggests that it is an early event in uveal melanoma (Onken et al., 2008). Mutations affecting Q209 in GNAQ were present in 45% of primary uveal melanomas and 22% of uveal melanoma metastases (Van Raamsdonk et al., 2010). Our analysis also shows frequent mutations of GNAQ in the TCGA UVM (uveal melanoma) cohort. Recently, of the 11,111 patients screened, 117 patients have been found to harbour GNAQ/GNA11 mutations, in melanoma, colorectal, liver, glioma, lung, bile duct and gastric cancers (Yang et al., 2020). GNA11 encodes subunit α-11 in the Gq class and acts as an activator of PLC. Mutations affecting Q209 in GNA11 were present in 32% of primary uveal melanomas and 57% of uveal melanoma metastases (Van Raamsdonk et al., 2010). Our analysis also shows frequent mutations of GNA11 in the TCGA UVM (uveal melanoma) cohort. GNAS encodes subunit α in the Gs class and participates in the activation of adenylyl cyclases, resulting in increased levels of the signalling molecule cAMP. GNAS functions downstream of several GPCRs, including beta-adrenergic receptors. GNAS mutations are found in 67% of intraductal papillary mucinous neoplasms and many associated pancreatic ductal adenocarcinomas (Takano et al., 2014). High GNAS expression in a breast tumour tissue showed a close correlation with the reduced overall survival, frequent distal metastasis, advanced clinical stage, stronger cell proliferation and enhanced cancer cell migration (Jin et al., 2019). Recently, it has been shown that GNAS promotes the development of small cell lung cancer via PKA (Coles et al., 2020). Our analysis shows frequent mutations and amplifications of GNAS in TCGA COAD (colon adenocarcinoma), LIHC (liver hepatocellular carcinoma) and READ (rectum adenocarcinoma) cohorts. The current knowledge on cancer-associated alterations of GPCRs and G proteins has been recently reviewed (Larribère & Utikal, 2020). Strikingly, approximately 36% of all drugs approved by the US Food and Drug Administration during the past three decades target GPCRs (Rask-Andersen, Almén & Schiöth, 2011).

Two members of isocitrate dehydrogenase family, IDH1 and IDH2, appeared on the top 10 SNA-based oncogenic events list as ranked by NDSI. The protein encoded by IDH1 is the NADP⁺-dependent isocitrate dehydrogenase found in the cytoplasm and peroxisomes. The cytoplasmic enzyme serves a significant role in cytoplasmic NADPH production. The protein encoded by IDH2 is the NADP⁺-dependent isocitrate dehydrogenase found in the mitochondria. It plays a role in intermediary metabolism and energy production. The most frequent mutations R132 (IDH1) and R172 (IDH2) involve the active site and result in simultaneous loss of normal catalytic activity, the production of α-ketoglutarate, and gain of a new function, the production of 2-hydroxyglutarate (Yan et al., 2009; Ward et al., 2010; Yang et al., 2012; Bendahou et al., 2020). 2-hydroxyglutarate is structurally similar to α-ketoglutarate, and acts as an α-ketoglutarate antagonist to competitively inhibit multiple α-ketoglutarate–dependent dioxygenases, including both lysine histone demethylases and the 10-11 translocation family of DNA hydroxylases (Yang et al., 2012). Abnormal histone and DNA methylation are emerging as a common feature of tumours with IDH1 and IDH2 mutations and may cause altered stem cell differentiation and eventual tumorigenesis (Yang et al., 2012). In acute myeloid leukaemia, IDH1 and IDH2 mutations have been associated with the worse outcome, shorter overall survival, and normal karyotype (Marcucci et al., 2010). All the 1p19q co-deleted gliomas have mutations in IDH1 or IDH2 (Labussière et al., 2010). Our analysis shows frequent mutations of IDH1 and IDH2 in the TCGA LGG (lower grade glioma) cohort and frequent amplifications of IDH1 in LIHC (liver hepatocellular carcinoma), as well as less frequent mutations and amplifications of IDH1 in CHOL, GBM, PRAD and SKCM.

Two fibroblast growth factor receptors FGFR2 and FGFR3 appeared on the top 10 mixed oncogenic events list as ranked by NDSI. The extracellular region of these proteins, composed of three immunoglobulin-like domains, interacts with fibroblast growth factors, leading to the activation of a cytoplasmic tyrosine kinase domain that phosphorylates PLCG1, FRS2 and other proteins. This sets in motion a cascade of downstream signals, including RAS-MAPK and PI3K-AKT pathways, ultimately influencing cell proliferation, differentiation, migration and apoptosis. FGFR aberrations were found in 7.1% of cancers, with the majority being gene amplification (66% of the aberrations), followed by mutations (26%) and rearrangements (8%) (Helsten et al., 2016). FGFR1 was affected in 3.5% of 4,853 patients; FGFR2 in 1.5%; FGFR3 in 2.0%; and FGFR4 in 0.5% (Helsten et al., 2016). The cancers most commonly affected were urothelial (32% FGFR-aberrant), breast (18%), endometrial (∼13%), squamous lung (∼13%) and ovarian (∼9%) (Helsten et al., 2016). Our analysis also shows frequent mutations and amplifications of FGFR2 in TCGA LUSC (lung squamous cell carcinoma) and UCEC (uterine corpus endometrial carcinoma) cohorts, as well as frequent mutations and amplifications of FGFR3 in BLCA (bladder urothelial carcinoma), HNSC (head and neck squamous cell carcinoma) and LUSC (lung squamous cell carcinoma) cohorts.

Some of the drivers prioritized by NDSI are highly tissue-specific. For example, GTF2I is specific to thymomas, whereas GNAQ and GNA11 are specific to uveal melanomas. As discussed above, in many cases a mutation in one of these genes is able to singlehandedly drive cancer in the corresponding tissue (Amaro et al., 2017). This undoubtedly makes them very strong drivers, despite being able to drive cancer in only one tissue type. A contrasting example would be TP53 which is the most common cancer driver in many tissues, but its individual strength is estimated as relatively low by NDSI, as it requires many other drivers to cooperate for cell transformation. Thus, it is important to differentiate the “weak – strong” axis from the “tissue-specific – universal” axis, as these appear to be orthogonal.

A puzzling question that remains in cancer genomics is why mutations in a given driver gene are typically confined to one or a few cancer types, resulting in each cancer type having its own unique set of driver genes (Iranzo, Martincorena & Koonin, 2018). As mutations are supposed to happen randomly as a result of stochastic mutagenesis processes (Belikov, 2017; Belikov, Vyatkin & Leonov, 2021), it is logical to suggest that mutations in different tissues can affect the same genes. However, the same mutation can be selected for in some tissues and selected against in others (Levine, Jenkins & Copeland, 2019). This selection most likely depends on the tissue-specific epigenetic profiles and microenvironments of the cancer-initiating stem or progenitor cells (Hass, von der Ohe & Ungefroren, 2020; Shlyakhtina, Moran & Portal, 2021). Thus, investigating the interplay between stem cell mutations, epigenetic profiles and microenvironments in various tissues appears to be a promising and exciting avenue for future research.

While both DSI- and NDSI-ranked top 50 genes are significantly overrepresented in such Reactome categories as signalling by PTK6, ESR-mediated signalling, PIP3 activates AKT signalling, signalling by receptor tyrosine kinases, signalling by erythropoietin, RAF/MAP kinase cascade, generic transcription pathway, chromatin modifying enzymes, and diseases of signal transduction by growth factor receptors and second messengers, there are also multiple differences. Top 50 DSI-ranked genes are additionally overrepresented in signalling by NOTCH, signalling by WNT, signalling by TGF-beta receptor complex, mitotic cell cycle, meiosis, cell cycle checkpoints, DNA double-stand break repair, epigenetic regulation of gene expression, RNA polymerase I transcription, circadian clock, diseases of cellular senescence, diseases of programmed cell death, and cellular responses to stress. This suggests that although these pathways are frequently mutated in cancer, none of them possesses strong tumour-promoting activity on its own. On the other hand, top 50 NDSI-ranked genes are additionally overrepresented in GPCR downstream signalling, which suggests that although this pathway is mutated more rarely in cancer, it nevertheless has a very strong tumour-promoting activity. It is also peculiar why neither DSI- nor NDSI-ranked top 50 genes are significantly overrepresented in DNA replication, autophagy, and metabolism categories. This may indicate that the role of these processes in oncogenesis is overestimated.

The major signalling pathway activated by mutations in both top DSI- and top NDSI-ranked driver genes is the RAS-RAF pathway. Although the pathway can be triggered via mutations in EGFR, FGFR, NRAS, HRAS, KRAS and BRAF genes, all of which are in the top 50 of both DSI and NDSI rankings, it can be additionally engaged through mutations in the top DSI-ranked driver ERBB2 and the top NDSI-ranked drivers GNAQ and GNA11. This suggests that ERBB2 driver mutations occur more frequently but are weaker than GNAQ and GNA11 driver mutations. Also, top DSI-ranked driver mutations affect the upper part of the PI3K-AKT-MTOR pathway via constitutive PIK3CA activation or PTEN inactivation, whereas top NDSI-ranked mutations affect the lower part of the pathway by activating AKT1 and MTOR. Similarly, this suggests that PIK3CA and PTEN driver mutations occur more frequently but are weaker than AKT1 and MTOR driver mutations. Moreover, the CTNNB1-MYC pathway, TP53-CDKN2A-RB1 pathway and MYC-CUL1-RB1 pathway are engaged only by top DSI-ranked drivers, indicating their relative weakness in cancer promotion despite high frequency of mutation, whereas the TCEB1-VHL-HIF1A pathway—only by top NDSI-ranked drivers, suggesting that this pathway has very strong tumour-promoting potential whilst being mutated more rarely.

It is interesting to discuss how driver strength should be estimated for genes responsible for the mutator phenotype (Loeb, 2011). On one hand, it appears that mutations in genes responsible for DNA repair should be considered strong drivers as they initiate a cascade of further mutations, many of which occur in other driver genes, thus heavily promoting cancer. On the other hand, by our definition of driver strength, if a driver needs help from many other drivers it cannot be called strong. Thus, it may be more fitting to call them “trigger drivers”, as they trigger the activation of other necessary drivers, rather than call them strong drivers per se. This may explain why top NDSI-ranked genes are not enriched in the Reactome DNA repair pathways. Our algorithm “sees” that there are always too many other driver mutations alongside the DNA repair drivers, thus it classifies the latter as weak. It might be viewed as a pitfall of our method, but we think it just faithfully follows the definition of what a strong driver really is, i.e. a driver able to drive cancer on its own or with a couple more drivers, as opposed to a “trigger driver” which recruits many additional helper drivers.

Another interesting aspect to discuss is the relation of driver strength to the evolutionary history of cancer. i.e. the sequential order of appearance of driver mutations. While our methodology does not allow deciphering the order of driver mutations, this subject was investigated, for example, in a recent article by Gerstung et al. (2020). There appears to be some indication that strong drivers (according to our methodology) overlap with early drivers (according to their methodology). For example, in Fig. 2B they show that, amongst the 50 most recurrent driver lesions, SMARCA4 and NF1 have the highest significant odds ratio (>50) of early versus late clonal driver mutations, and SPOP has the highest significant odds ratio (>50) of clonal vs. subclonal driver mutations. In our results (Table 2), NF1 is the 50th top DSI-ranked gene, whereas NF2 (not shown in Gerstung et al. (2020)), SMARCA4 and SPOP are the 17th, 27th and 35th top NDSI-ranked genes. Many other significantly early/clonal genes—KRAS, VHL, BRAF, PBRM1, SETD2, NRAS, MEN1, IDH1, ATRX—were also found amongst the top 50 NDSI-ranked genes. Thus, despite the bias for recurrence in the Gerstung et al. (2020) figure, many of these genes were recovered in recurrence-independent ranking by NDSI. This suggests that strong drivers might indeed be activated early, which fits with our definition that strong drivers are able to initiate cancer either on their own or with little help.

Finally, it is important to note that our approach is not meant to replace frequency-based ranking, because it serves a different purpose. While the frequency-based approach identifies the most common drivers, our method identifies the strongest ones.

Conclusions

We have introduced a novel concept of cancer driver strength, formulated algebraic equations for Driver Strength Indices, wrote software to calculate these indices and applied it to TCGA PanCanAtlas datasets. As a result, we presented a comprehensive overview on the landscape of cancer driver genes and chromosomes in TCGA PanCanAtlas patients and highlighted particular genes, gene families and pathways deemed strong drivers according to our Driver Strength Indices. These findings should help to direct future research efforts and selection of promising targets for novel therapeutics.

Supplemental Information