Benchmark enables identification of optimal settings and best practices for systematic pathway discovery in experimental datasetsWe developed Benchmark as an evaluation platform to explicitly evaluate performance of pathway/geneset analysis tools (Fig. 1a, b; please note that the terms ‘pathway’ and ‘geneset’ are used interchangeably throughout the manuscript). Briefly, Benchmark comprises three components: 1) “input genesets” (IGS), representing genesets to be investigated – these are most commonly user lists of genes derived from high throughput assays, such as RNA-seq or Chromatin Immunoprecipitation sequencing (ChIP-seq) performed during projects of interest; 2) “target genesets” (TGS), representing curated biological pathways (e.g. in the KEGG and Gene Ontology databases) to be examined for their relationship to the IGS; 3) pathway enrichment analysis tools which test statistical relationships between IGS and TGS (Fig. 1a). The genesets used in Benchmark were extracted from large numbers (~1000) of high-throughput sequencing experimental assays from ENCODE10 (Supplementary Data 1a). Every geneset in Benchmark was composed of genes bound in ChIP-seq by a defined transcription factor (TF), in eCLIP-seq by a defined RNA binding protein (RBP), or in RNA-seq as genes differentially expressed after small interfering or hairpin RNA-mediated knockdown of a specified target (gKD). For each TF, RBP or gKD, at least two genesets from distinct cell lines (e.g. K562 and HepG2) or species (e.g. human and mouse) were represented. Importantly, Benchmark was designed to simulate a common scenario where the investigated genesets (i.e. IGSs) are derived from different sources (e.g. cell lines, species, conditions) than established biological pathways (i.e. TGSs). To achieve this, Benchmark was created such that genesets from the same TF, RBP or gKD but from distinct cell lines are better matched to each other than any other genesets. As an example, the IGS that contains STAT1-bound genes in K562 cells is more closely related to the TGS that contains STAT1-bound genes in HepG2 cells than any other TGS.Fig. 1: Benchmark for evaluating and optimizing existing pathway analysis methods.a Schematic representation of Benchmark for each pathway analysis method. For every input geneset (IGSi), the outputs, of pathway enrichment analysis against all target genesets (TGSi) are recorded in the score matrix M. Benchmark is created such that TGSi is the most related geneset to IGSi compared to all other TGSs. Thus, the rank of diagonal entries in M represent how well the correct genesets are scored for each IGS. The best score is 1. b, c Cumulative distribution function (CDF) plots of correct pathway ranks comparing indicated pathway analysis methods on shRNA-seq genesets. Median ranks for correct pathways, the fraction of correct pathways among top 10 reported pathways (precision@10/P@10, and average precision@10/AP@10) are shown. d Plots showing the correct pathway ranks of the indicated methods in shRNA-seq genesets. Boxes span from the first to the third quartile, with the median highlighted, and whiskers covering the minimum to maximum values. The precision statistics of each method at finding the correct pathway amongst reported pathways are reported below each bar. Target genesets in b–d were defined as the top 200 upregulated genes by DESeq2 in K562 cells; input genesets in b–d for default and optimized settings were extracted using signal2noise (default setting of GSEA) and DESeq2 in HepG2 cells, respectively. e Plots showing the performance of multiple methods as increasing percentage of unrelated random genes (“noise”) are introduced to target genesets. This evaluates the robustness of methods to inaccuracies (or noise) in target genesets. Lines and error bars show median +95% confidence interval at each noise level. *p < 0.05; **p < 0.01; ****p < 0.0001 by two-sided Wilcoxon rank sum test corrected for multiple hypothesis. ora! and ora are essentially identical, with the only difference being that they are implemented using different programming languages. Specifically, ora! is implemented in R by egsea, while ora is implemented in Python by PET. Source data are provided as a Source Data file.Consequently, for each IGS, an ideal pathway analysis method should identify the most closely matched (i.e. “correct”) TGS and to rank it as number one compared to all other TGSs. Thus, in Benchmark, the enrichment of TGS1 to n is calculated for each IGS using the user’s pathway analysis method of choice to generate a score matrix, from which the rank of correct pathways can be determined for that method (Fig. 1a). Thus, by evaluating how accurately pathway analysis methods rank matched IGSs and TGSs compared to all other TGSs, Benchmark can assess the performance of different methods in identifying and ranking correct biological pathways in experimental settings.Using Benchmark, we first evaluated the performance of different pathway analysis tools, namely, decoupler11, piano12, egsea6, ssgsea4, camera13, ora6, enrichment browser14 (e.browser), zscore15, gsva16, padog17, roast18, plage19, gage20 and safe21. We compared three statistics for all pathway analysis methods. First, we determined the median rank of the correct pathway returned by each method. As experimental biologists typically focus on the top 10 reported pathways, we also determined the frequency by which the correct pathway is among the top 10 reported pathways, termed “Precision@10” (P@10). We also adopted another reliable performance metric commonly used in information retrieval, known as Average Precision at 10 (AP@10). AP@10 is the average (mean) of the precision scores at each of the first 10 positions, i.e. P@1, P@2,…, P@10. We tested four different settings for how the genes can be ranked and genesets can be generated in these methods (i.e. ratio of classes or signal2noise or edgeR- or DESeq2-based ranking (Fig. 1b, Supplementary Fig. S1a). The genes generated from matched IGSs and TGSs using these settings exhibited a significantly higher degree of overlap compared to randomly generated control IGSs and TGSs (Supplementary Fig. S1b). The top three methods, decoupler11, piano12, egsea6, notably all ensemble approaches, achieved a median correct pathway rank of 1–8, P@10 of 52–76%, and AP@10 of 44–69%, under the best setting (Fig. 1b). This creates a scenario where pathways with significant biological importance are regularly not among the top 10 reported pathways, hindering the utility of these methods for unbiased discovery.We next focused on evaluating the performance of the three most used individual pathway analysis methods ora6, GSEA4 and Enrichr5 in more detail. We tested the four different settings described above (Supplementary Fig. S1c) as well as 3 different options for sizes of genesets (i.e. 200, 500 and 1000 genes) (Supplementary Fig. S1d). Specifically, under the default settings (i.e. signal2noise ranking), the median ranks of the correct pathway (7–14) and P@10 (45–54%) values were similar to the top three performing methods above (Fig. 1c, Supplementary Fig. S1c), indicating that ora, GSEA and Enrichr are not ideal for unbiased pathway discovery under default settings. We found that optimization of input parameters did significantly improve pathway discovery – a geneset size of 200 genes was near optimal for both GSEA and Enrichr (Supplementary Fig. S1d) and using DESeq2-ranked lists of differentially expressed genes returned the correct pathway at a median rank of 3 and 63–68% of the time among the top 10 reported pathways (Fig. 1c, Supplementary Fig. S1c). Similar improvement could be achieved when selecting the optimal binding intensity statistics in ChIP-seq or eCLIP-seq. For example, using the signal value of binding in ChIP-seq data to rank genes and create genesets usually leads to an improved pathway discovery performance compared to utilizing adjusted p-values (or q-values) of binding signals (Supplementary Fig. S1e). It is worth noting that the overlap size of matched IGSs and TGSs generated from ChIP-seq and eCLIP-seq data was significantly larger compared to that of shRNA-seq data (Supplementary Fig. S1f), which could explain the overall trend of better performance observed in pathway analysis methods on these genesets.Collectively, these data indicate that Benchmark can be used for evaluating existing and future pathway analysis methods for discovery. Notably, although the most widely used individual methods can be optimized, they remain sub-optimal for unbiased pathway discovery, and ensemble approaches typically outperform these individual methods.Pathway Ensemble Tool (PET) enables unbiased pathway discoveryWe next sought to develop a high precision method to enable unbiased systematic pathway discovery. Leveraging the accuracy and robustness of ensemble approaches, we created Pathway Ensemble Tool (PET). PET works by calculating the unweighted mean of the rank statistics from GSEA, ora, and Enrichr, all operating under the optimal settings determined above (see Methods). We tested the accuracy of PET using Benchmark. PET significantly improved the median rank of the correct genesets to 1, indicating that PET tends to rank the correct pathway as number 1 (Fig. 1d). PET also significantly enhanced P@10 and AP@10 statistics, with 73–82% of the time the correct pathways being among the top 10 pathways reported by PET (Fig. 1d).Pathway analysis tools are expected to conform to the false discovery rate (FDR) estimates, meaning that at 5% FDR levels, 5% of reported pathways are expected to be false. To determine the actual false discovery rates of pathway analysis methods, we randomly sampled and divided replicated RNA-seq datasets into two groups and ran pathway analysis methods comparing the two groups for enrichment of biological pathways (see Methods). Since the samples were randomly selected biological replicates, any reported pathways significantly different between the groups represent false discoveries. Surprisingly, we found that some of the pathway analysis methods had much higher than the expected false discovery rate, in some cases as high as 25% (Supplementary Fig. S1g, h). In fact, specific pathways, such as the KEGG Ribosome pathway, were common false discoveries (Supplementary Fig. S1i and Supplementary Data 1b), suggesting that caution should be exercised when these pathways are identified as significant during unbiased discovery approaches. However, we found that PET operated at the expected 5% false discovery rate (Supplementary Fig. S1g, h).To further validate PET reported results and ensure its ability to discern pathways associated with alterations that drive differences in transcriptional programs within groups, we applied PET to compare fibroblast growth factor receptor 3 (FGFR3) mutant vs. wild-type (WT) bladder cancer cell lines, estrogen receptor (ER)+ vs. ER– breast cancer cell lines, and epidermal growth factor receptor (EGFR) mutant vs. WT lung cancer cell lines, available in CCLE22 (Supplementary Data 1c), assessing all canonical pathways. In all three comparisons, the top reported pathways were supported by the literature (Supplementary Fig. S2). For instance, the top two enriched pathways in ER+ breast cancer cell lines were “Estrogen dependent gene expression” and “Estrogen receptors (ESR)-mediated signaling”, which corresponds with the established understanding of estrogen receptor signaling in ER+ breast cancer23,24. In EGFR mutant lung cancer cell lines, the top two enriched pathways were “Surfactant metabolism” and “Integrin5 pathway”, which are both recognized for their cross-talk with EGF signaling in lung cancer25,26,27,28,29. Moreover, the “MYC activity pathway” and its regulated “Ribosomal RNA processing” pathway were among the top enriched pathways in FGFR3 mutant bladder cancer cell lines, aligning with known interactions between FGFR3 and MYC in bladder cancer, as supported by existing literature30,31,32,33,34 (Supplementary Fig. S2). These data demonstrate that PET is able to discern subtle yet significant variations in pathway activities across different genetic backgrounds.Because most pathways are defined under real experimental conditions, there are multiple variables, such as cell type, species and stimulation, that contribute a degree of differences or inaccuracies in genesets. This “noise” manifests as missing or unrelated genes among matched IGSs and TGSs. To test the robustness of common pathway analysis tools, we evaluated their performance under conditions of simulated noise by replacing genes from the target genesets with increasing proportions of unrelated genes selected randomly from the genome. We found most available pathway analysis methods to be highly sensitive to such differences as their performance significantly dropped after >30% of unrelated random genes were introduced into existing genesets (Fig. 1e). Importantly, however, PET was highly robust to even as high as 60% added noise (Fig. 1e). Collectively, these results indicate that PET is robust, has controlled false discovery rate and high precision, thus can enable reliable and systematic identification of key dysregulated pathways in omics datasets.PET identifies prognostic pathways in 12 distinct cancer typesWe next sought to deploy PET for unbiased pathway discovery in cancer. We selected 12 distinct cancer types from the cancer genome atlas (TCGA) that had cancer tissue RNA-seq data from the earliest stages of disease (based on the pathological stage or primary tumor size/extent) in at least 20 subjects that had unfavorable prognosis (i.e. succumbed to death) during the follow-up period (Fig. 2a; Supplementary Data 2a, b). We used PET to identify the pathways enriched in samples from patients with favorable (survivors) or unfavorable (deceased) prognosis for each cancer type, since these pathways could act as potential biomarkers of outcomes. To ensure that age, sex and cancer stage were not confounders, patients with favorable prognosis were randomly matched with those with unfavorable prognosis (see Methods). We confirmed that the groups were matched in age and sex and that age and sex were not determinants of outcome (Supplementary Fig. S3). Since tumor infiltrating lymphocytes, classically depicted by CD8-expression, are known determinants of outcome, we also noted that CD8 expression did not differ between subjects with favorable and unfavorable outcomes and that CD8 expression could not distinguish survivors from deceased subjects (Supplementary Fig. S4a). We next examined >1000 canonical biological pathways curated in MSigDB4 that could potentially differentiate cancer patients that survived from those that did not using PET. Out of these pathways, PET identified tens of biological pathways significantly associated with favorable or unfavorable prognosis for each cancer type (Fig. 2b, Supplementary Data 2c). Examples of the top favorable and unfavorable prognostic pathways in each cancer types are shown in Supplementary Fig. S4b and Fig. 2c, respectively. As proof-of-principle, we validated the reported pathways associated with unfavorable outcome in bladder cancer in three independent cohorts35,36,37 (Supplementary Fig. S4c). Some pathways were significantly dysregulated in multiple different cancers (e.g. cell cycle in bladder carcinoma (BLCA), cervical squamous cell carcinoma (CESC), kidney renal papillary cell carcinoma (KIRP) and pancreatic adenocarcinoma (PAAD)) whereas others were cancer-specific (e.g. glycerolipid metabolism in breast invasive carcinoma (BRCA) and arginine and proline metabolism in KIRP. Some of these pathways have been shown to have biological significance in previous studies but have not yet been formally recognized as prognostic markers. For example, we identified triacylglyceride synthesis as a significant poor prognostic pathway in breast carcinoma (Fig. 2c). This is consistent with fatty acid synthesis as a requisite for breast cancer metastasis to brain38 and association of upregulated cellular triacylglycerol with enhanced survival of human breast cancer cells39. Another identified poor prognostic pathway that was shared between multiple cancer types including bladder and pancreatic cancers was extracellular matrix (ECM) organization (Fig. 2c). This is in line with prior knowledge about the role of ECM components, such as laminin, collagen and fibronectins in invasion, progression, and metastasis of bladder and pancreatic cancers and restriction of antitumor drug delivery40,41. Similarly, we found that Vitamin D (VitD) receptor pathway to be a strong favorable prognostic pathway in stomach cancer (Supplementary Fig. S4b). This is congruent with the fact that in gastric cancer patients, VitD levels serve as an independent prognostic factor and VitD supplementation can reduce the risk of developing gastric cancer42,43.Fig. 2: PET identifies pathways and gene combinations associated with unfavorable prognosis.a Cancers used in this study, along with tissue type, pathological stage (roman numerals) or pathological primary tumor TNM T score (e.g. T1, T2, T3) and number of samples. BLCA bladder carcinoma, BRCA Breast invasive carcinoma, CESC Cervical squamous cell carcinoma, COAD Colon adenocarcinoma, HNSC Head and Neck squamous cell carcinoma, KIRC Kidney renal clear-cell carcinoma, KIRP Kidney renal papillary cell carcinoma, LIHC Liver hepatocellular carcinoma, LAUD Lung adenocarcinoma, LUSC Lung squamous cell carcinoma, PAAD Pancreatic adenocarcinoma, STAD Stomach adenocarcinoma. To adjust for confounders, samples from later stages of diseases have also been included in cancers denoted by #. b Plot showing the number of significantly enriched pathways identified using PET among >1000 MSigDB canonical pathways in patients with favorable (blue) or unfavorable (red) prognosis in the indicated cancers. c Top pathways from each cancer associated with unfavorable survival. FDR values are from PET. d Area under the curve (AUC) distribution of indicated sets used to predict deceased status in each cancer type. Leading genes were selected from top 20 significant pathways associated with unfavorable prognosis. DEGs, differentially expressed genes; error bars show mean ± s.d. AUCs of leading gene combinations are higher than all genes and the AUCs of 5-gene combinations are higher than individual leading genes (Wilcoxon rank sum test p-value < 1e−9) (Supplementary Data 2d). e Elbow plots showing significant AUCs (low to high) from all 1–5 leading gene combinations. Statistics were based on FDR-adjusted logrank test of Kaplan–Meier (KM) survival curves. BRCA and LIHC did not yield any significant unfavorable predictors. KM plots showing of the top combination biomarker (darker lines in f), the top DEG (lighter lines in f) of unfavorable survival, PAI-1 expression (orange lines in g) and uPA expression (blue lines in g) to stratify overall survival. Samples in f and g were evenly split into two groups according to z-score/expression values of the indicated biomarker. Hazard ratio (HR) and logrank test p-values are reported. *p < 0.01; **p < 0.001; ***p < 0.0001. See Supplementary Fig. S4 for information on favorable biomarkers. Source data are provided as a Source Data file.PET-identified biological pathways guide combination biomarker discoveryWe next explored the genes in significant prognostic pathways to identify a subset of genes with high prognostic power that might have clinical utility. Single gene expression values have provided some level of prognostic power for different cancers44,45. However, expression of multiple biomarkers are consistently more accurate diagnostic/prognostic tools to stratify patients and improve cancer management46,47. This has fueled a growth in studies identifying gene signatures associated with response to therapy and/or survival. However, these signatures are typically derived from lists of differentially expressed genes and characteristically consist of hundreds of genes, making clinical application of these large signatures challenging48. The alternative brute force strategy, which involves iteratively uncovering combinations of genes with prognostic power among ~20,000 protein-coding genes in the genome, is computationally challenging (or impossible) and statistically unfeasible. For example, an exhaustive search space for a 3-gene combination has 8 trillion permutations, which is statistically unfeasible given the number of cancer samples available for study. We postulated that limiting the search space of combinatorial biomarker discovery to leading genes from prognostic pathways might alleviate this problem. We therefore selected up to 20 of the most significant prognostic pathways reported by PET and extracted their leading genes from the corresponding GSEA results. We then generated all possible 1–5 gene combinations from the leading genes and quantified prognostic power by calculating the area under the curve (AUC) statistic of the receiver operating characteristic (ROC) curve for the ability of a gene or gene combinations to distinguish survivors from deceased subjects (see Methods). Combinations of 3–5 leading genes from prognostic pathways were found to be the best predictor of prognosis. Notably, none of the individual gene biomarkers displayed a higher AUC than the combinations identified. These are shown in Fig. 2d, e, Supplementary Fig. S4d, e and Supplementary Data 2d, e, together with their corresponding overall survival plots (Fig. 2f, Supplementary Fig. S4f). As an example, high expression of two collagen genes (COL5A1, COL3A1) and vitronectin (VTN) provides a significantly unfavorable overall survival outcome (<20% 10-year survival; HR 9.5; p < 0.0001) in renal papillary cancers (KIRP; Fig. 2e, f). The ability of this 3-gene combination to distinguish survivors from non-survivors over a prolonged period of follow-up is superior to any other biomarkers of renal cancer in current use49. Previous studies have linked the expression of individual collagen genes such as COL5A1 and COL23A1 to renal cancer50,51. Collagen expression has also been associated with tumor grade52 and plays a key role in promoting tumor cell invasion53 in renal cancers. Furthermore, vitronectin is linked to inducing the differentiation of cancer stem cells and the formation of tumors54. Thus, it is biologically conceivable that high expression of this 3-gene combination biomarker provides the observed dramatically unfavorable overall survival outcome, surpassing the prognostic power of individual collagen genes alone.Moreover, the combination biomarkers outperformed any individual gene including differentially expressed genes (DEGs; Supplementary Data 2f) in all cancers, and no individual gene was found to be the most effective on its own (Fig. 2d, e, Supplementary Fig. S4d, e). Illustrative examples are indicated in Fig. 2f and Supplementary Fig. S4f, showing Kaplan–Meier overall survival curves of the top combination biomarker compared to the top DEG in each cancer type. As further illustration, we compared the performance of the PET-derived combination biomarkers (Fig. 2f) to that of PAI-1 (Plasminogen activator inhibitor 1) and uPA (Urokinase plasminogen activator) (Fig. 2g), because these two genes are frequently investigated as potential cancer prognostic markers55. PAI-1 is associated with poor prognosis in some cancer types, such as breast, colon, and lung cancer. uPA is a protease involved in the degradation of the extracellular matrix, and linked to tumor invasiveness, angiogenesis, and poor prognosis in some cancers, including breast, prostate, and ovarian cancers. We found that although elevated expression of PAI-1 in CESC, LUSC, PAAD, and STAD, and uPA in PAAD may have prognostic value (Fig. 2g), they are both inferior to the identified top combination biomarkers (Fig. 2f). Moreover, the identified combination biomarkers accurately predicted outcomes across independent cohorts of patients. For example, >70% of biomarkers in TCGA cohorts exhibited consistent favorable or unfavorable overall survival outcomes in several independent cohorts of bladder35,36,37, breast56, pancreatic57, and stomach58,59 cancers (Supplementary Fig. S4g; Supplementary Data 2g). We also investigated the percentage of combination biomarkers that consistently exhibited predictive power in later-stage samples. We found that in certain cancer types (e.g. HNSC, PAAD), the early-stage biomarkers consistently predict outcomes even in later stages of cancer. However, in other cancer types (e.g. BLCA, KIRP), only a small fraction of the early-stage biomarkers demonstrate consistency in their predictive ability (Supplementary Data 2h, i).Mutations play an important role as biomarkers in clinical settings, therefore we sought to compare the common cancer mutations versus gene expression-based biomarkers for their power in determining overall survival outcome. To this end, we focused on the three most frequently mutated genes in cancer samples, BRAF, TP53 and KRAS, that often confer poor prognosis60,61,62. For each TCGA cancer type with sufficient number of samples exhibiting mutations (n > 15), we split the samples into two groups based on their mutational status: wild-type samples without the mutation and samples containing mutant oncogenes. To assess whether these mutations could serve as indicators of survival, we conducted log-rank tests between the wild-type (WT) samples and mutated (mut) samples for each cancer type. Among all the cancer types examined, pancreatic cancer (PAAD) was the only one in which the presence of KRAS mutation served as a significant indicator of poorer overall survival. (Hazard ratio = 1.8; logrank p-value = 0.02) (Supplementary Fig. S4h). Nevertheless, the KRAS mutation exhibited less discriminatory power compared to gene expression-based biomarkers in distinguishing between favorable (Supplementary Fig. S4f) and unfavorable (Fig. 2f) survival outcomes.We also assessed the susceptibility of the identified biomarkers or conventional DEGs and their combinations to overfitting, a commonly encountered pitfall in the search for trustworthy biomarkers that leads to poor performance when tested on different samples. We found that combination biomarkers were significantly less prone to overfitting compared to those derived from conventional DEGs (Supplementary Fig. S4g). This indicates that combination biomarkers have a better balance between complexity and robustness, which makes them more suitable for real-world applications. Collectively these data indicate that combination biomarkers identified using PET have great potential to stratify patient outcomes and perform better than existing known biomarkers.Identifying resilient and vulnerable cancer patients is crucial as it allows healthcare providers to personalize treatment plans, improve prognosis predictions, prioritize resources, and match patients with appropriate clinical trials. This leads to more effective and efficient care, better outcomes, and improved quality of life for patients. We thus considered whether combining the top biomarkers of favorable with the top biomarkers of unfavorable outcome would provide added prognostication power. To do this, we divided the samples into two groups based on expression of the top favorable combination biomarker and two groups based on expression of the top unfavorable combination biomarker. By overlapping these groups, patients can be further divided into subgroups according to simultaneous low and high expressions of the top biomarkers of favorable and unfavorable survival. The resilient category was defined as having high levels of the top favorable combination biomarker and low levels of the top unfavorable combination biomarker, while the vulnerable category was defined as having low levels of the top favorable combination biomarker and high levels of the top unfavorable combination biomarker. Indeed, this approach identified highly resilient and highly vulnerable cancer patients (Fig. 3a, Supplementary Fig. S5a, Supplementary Data 3a). For example, resilient patients with papillary or clear cell kidney cancer (KIRP and KIRC) had nearly 90–100% survival, while vulnerable patients had only 10–20% overall survival over a prolonged period of follow-up (~4000 days) (Fig. 3b). This prognostication strategy also corroborated with established independent prognostic biomarkers in these patients63. Specifically, vulnerable patients have more than two-fold higher frequency of repressed CDKN2A expression via copy number loss or DNA methylation when compared to resilient patients (Fig. 3c, Supplementary Data 3b). The same approach in other cancers also yielded added benefit, because in almost all cancers combining the top biomarker of favorable and the top biomarker of unfavorable outcome provided highly accurate prognostication (Supplementary Fig. S5a, Supplementary Data 3a). For instance, the 10-year overall survival rate for resilient patients with bladder cancer is nearly 9 times higher than that of vulnerable patients (HR: 4.6; p < 0.0001; Supplementary Fig. S5a). Such accurate prognostication can have a major influence on clinical decision making and treatment course and can potentially help rationalize less or more aggressive individualized treatment for patients within different vulnerability.Fig. 3: Juxtaposition of favorable and unfavorable prognostic biomarkers offers precise overall survival stratification.a Venn diagrams showing the distribution of patients with kidney renal papillary cell carcinoma (KIRP) (left) and kidney renal clear cell carcinoma KIRC (right) based on expression of the top favorable (labeled in Supplementary Fig. S4e) and the top unfavorable (labeled in Fig. 2e) combination prognostic markers on RNA-seq from tumor samples. b Kaplan–Meier (KM) overall survival curves for 4 subgroups of KIRP or KIRC patients whose tumors express indicated combination of prognostic biomarkers. The four subgroups are highlighted with a bold font in a. Hazard ratio (HR) and logrank test p-values denote comparisons with the lines in red. ns not significant; *p < 0.05; **p < 0.01; ***p < 0.001; ****p < 0.0001. c Plots showing the fraction of KIRP (left) and KIRC (right) tumors which have alterations (methylation or copy number loss) in CDKN2A63 in each group. Source data are provided as a Source Data file.The molecular subtype of tumors, particularly basal and luminal subtypes, is one of the factors implicated in predicting the outcome of bladder cancer. The basal subtype is associated with muscle invasiveness and worse overall prognosis but is more responsive to platinum-based chemotherapy64. To verify if identified prognostic markers corroborate with this known prognostic factor, we sourced three independent cohorts of bladder cancer samples35,64,65 that was annotated with the two major molecular subtypes. Indeed, we found that combination biomarkers of unfavorable prognosis were expressed at high levels in the basal molecular subtype across all three cohorts (Supplementary Fig. S6a, Supplementary Data 3c).We had previously demonstrated using Benchmark that PET is highly robust to noise (Fig. 1e). To formally test this under real-world settings, we obtained RNA-seq datasets from a different species of mammal. We selected canine bladder cancers that included annotations for luminal and basal subtypes66 as these sub-types behave similarly to humans and are considered as having good and poor prognosis, respectively. The genes associated with poor prognosis in human bladder cancer effectively differentiated luminal and basal subtypes of canine bladder cancer as well (Supplementary Fig. S6b), demonstrating the robustness of these biomarkers even across different species. Collectively, these data indicate that our biomarker discovery approach can identify combination prognostic biomarkers with significant clinical value and to stratify molecular subtypes of cancer with differential prognosis.PET-identified biological pathways guide drug repurposingSince prognostic markers are derived from dysregulated pathways, we postulated that manipulating gene expression to promote the biomarkers of favorable prognosis and demote the biomarkers of poor prognosis by re-purposing known drugs might provide clinical benefit. Major consortium efforts have created a valuable library of integrated cellular signatures (LINC) that includes cellular transcriptional responses to large numbers of pharmacologic perturbagens67. Similar signatures have also been derived from mining large numbers of publicly available sequencing read archive (SRA) data68. Recently we showed that computational screening for drugs that normalize expression of genes in highly prognostic and pathogenic pathways using cellular signatures of those drugs provides significant therapeutic potential69,70. Thus, we computationally screened >1600 drugs (mostly FDA-approved) for which the up- or down-regulated drug targets are known68 (Supplementary Data 4a). Drugs were ranked by their ability to suppress unfavorable and to promote favorable prognostic biomarkers. The rank was calculated based on the significance of overlap between genes up/down regulated by a drug and genes from favorable/unfavorable prognostic pathways in cancer (Fig. 4a; see Methods). This in silico screening predicted several drugs already in clinical use for each cancer type, as well as many drugs that have not yet been tested (Fig. 4b and Supplementary Data 4b). Key interactions between highly significant unfavorable prognostic pathways and drugs that can normalize these are summarized for each cancer type in Fig. 4c. By way of validation, the top predictions included drugs currently in clinical practice or under trial for some of these cancers. For example, doxorubicin, a commonly used chemotherapy agent for the management of metastatic breast cancer71,72, was the second top prediction for breast cancer. Celastrol, the top prediction for colon cancer, was recently shown to inhibit proliferation of colorectal cancer cells and migration via the PI3K/AKT pathway73. Ruxolitinib, the top prediction for liver cancer, is known to have marked tumoricidal effects on hepatocellular carcinoma by inhibiting JAK/STAT signaling74. Interestingly, irinotecan, a topoisomerase inhibitor, which regulates the cell cycle, was the top prediction for pancreatic cancer (PAAD), a tumor type in which “cell cycle” is among the top pathways associated with unfavorable prognosis (Fig. 2c). Last year, the FDA approved investigational use of liposomal irinotecan for first line treatment of pancreatic cancer patients75,76,77. Two other pathways strongly associated with unfavorable outcome from pancreatic adenocarcinoma are the Erbb1 (also known as EGFR) and MYC pathways; both of these have previously been associated with worse outcomes from this cancer type78,79 and are also inhibited by irinotecan80. These data indicate that matching drug signatures and prognostic markers can be used for repurposing existing therapeutics.Fig. 4: PET-derived prognostic pathways guide drug repurposing.a Schematic of drug screening aims; computational screen searches for drugs that upregulate expressions of genes from pathways with favorable prognosis and downregulate expressions of genes from pathways with unfavorable prognosis. b Elbow plots showing the significance of overlap between drug up (or down) regulated genes with leading genes from pathways with favorable (or unfavorable) prognostic potential. Drugs are sorted based on the significance (see Methods). Top two drugs for each cancer type are highlighted. c Sankey plot showing top pathways associated with unfavorable prognosis for each cancer type on the left and top two predicted drugs that demote expression of genes in indicated pathway on the right. The thicknesses of links are proportional to the significance of the pathway (E-value > 5; left; see Methods) or the number of genes affected by drug (right). Only drugs that have a significant impact on the top prognosis pathways for each cancer are displayed. As a result, cancers without drugs affecting their top predicted prognosis pathways are not displayed. Source data are provided as a Source Data file.The CDK2/9 inhibition restricts cervical and bladder cells in vitro and/or in vivoWe next sought to test the efficacy of one of our predicted drugs. The top predicted drug for bladder and cervical cancer was a cyclin-dependent kinase (CDK)2/9 inhibitor named 0175029 (Fig. 4b). This was somewhat surprising because inhibitors of other CDKs, specifically CKD4/6 inhibitors, have previously been proposed as monotherapy for bladder cancer. The use of CDK4/6 inhibitors have been justified because elevated expressions of either CDK4 or CDK6 are associated with poor survival in bladder cancer, whereas high expression of CDK2/9 have been associated with better survival (Supplementary Fig. S7a). However, CDK4/6 inhibitors including Palbociclib did not meet the primary endpoint during phase II clinical trials for bladder cancer81,82,83,84. Moreover, taking the conventional DEG-based approach, where one normalizes the most differentially expressed genes between cancer survivors and non-survivors, does not identity CDK2/9 inhibitor as a significant lead in either bladder or cervical cancers (Supplementary Fig. S7b, Supplementary Data 4c). However, it is worth noting that although (CDK)2/9 inhibitor emerges as the top prediction for bladder cancer using DEG-based approach, it does not reach the statistical significance threshold (Supplementary Fig. S7b). We therefore scrutinized whether PET-based drug prediction specifically advocates CDK2/9 or CDK4/6 inhibition in bladder cancer. Since, the drug database used for screening above68 did not include any CDK4/6 inhibitors, we computationally screened an additional 5425 compounds from the LINC1000 database85 that included the CDK4/6 inhibitor Palbociclib. PET-based drug predictions did not advocate Palbociclib (p-value = 1; rank 1601) but returned 0175029, a CDK2/9 inhibitor as the number 1 scoring drug that can significantly normalize genes of prognostic pathways (Supplementary Fig. S7c, Supplementary Data 4d).For experimental validation, we then tested the ability of a panel of CDK inhibitors to restrict the growth of bladder and cervical cancers in vitro, including Palbociclib (a CDK4/6 inhibitor), GW8510 (a CDK2 inhibitor), and Seliciclib (a CDK2/7/9 inhibitor) (Supplementary Fig. S7d). Because 0175029 is not commercially available we used its recently identified analog CCT068127 which is also a potent CDK2/9 inhibitor86. Consistent with PET-based drug prediction, CCT068127, but not Palbociclib, significantly inhibited the growth of both cervical and bladder cancer cell lines, whereas the others either did not, or had minimal effect (Fig. 5a). We also found that in two out of three independent cohorts of bladder cancer samples, the expression of CDK9 (but not CDK2) was higher in the basal subtypes of bladder cancers compared to luminal subtypes (Supplementary Fig. S7e), suggesting that CDK9 inhibition might be more beneficial for the basal subtype of bladder cancer. CDK9 expression also showed a significant correlation (Pearson r = 0.71; p-value = 0.006) with unfavorability scores across bladder cancer cell lines available in the CCLE database22 (see Fig. 5b and Supplementary Data 5a). The unfavorability score for each sample was calculated as the mean (z-score) expression of genes associated with an unfavorable outcome (Supplementary Data 2e). To expand on this, we referenced publicly available data from the GDSC database87 where a range of bladder cancer cell lines were treated with a different CDK9 specific inhibitor (i.e. CDK9_5038) and IC50 values were calculated. Despite the difference in inhibitors, we observed a trend where cell lines with higher unfavorability scores required higher concentrations of the inhibitor to achieve similar inhibition of proliferation (Pearson r = 0.54, p < 0.06) (Fig. 5b). These associations suggest a broader efficacy of CDK9 inhibitors across bladder cancer cell lines, with efficacy correlating inversely with unfavorability scores.Fig. 5: Predicted CDK2/9 inhibition restricts the growth of cervical and bladder cancers in vitro and in vivo.a Plots showing the survival rate, measured by CCK-8 assay, of cervical (Hela) and bladder (T24 and UMUC3) cancer cell lines treated with different concentration of vehicle or indicated drugs for 72 h. The selective targets of CDK inhibitors are shown in parenthesis. Mean +/− sem from n > 3 biological replicates. b Plot showing expression of the CDK9 gene versus the unfavorability score for each bladder cancer cell line in the CCLE database22. The color scale represents the IC50 values reported by the GDSC database87 for the CDK9-specific inhibitor CDK9_5038. Cell lines for which IC50 values are not available are highlighted in gray. The Spearman correlation value and corresponding p-value are shown in red. The center line and error bars indicate the trend line and 95% confidence interval. Selected cell lines are labeled. c Heatmap showing differentially expressed genes (DEGs; 3-fold at adjusted p-value < 0.05) between DMSO or CCT068127 (1uM) treated T24 cells for 48 h in n = 3 independent biological replicates. GSEA plots showing enrichment of genes repressed by the CDK4/6 inhibitor Palbociclib (top) or genes repressed by the CDK9 inhibitor VIP152 (bottom) in T24 cells treated 48 h with the CDK2/9 inhibitor CCT068127 (1 uM) or DMSO control. Palbociclib (SRP404373) or VIP152101 inhibited genes are publicly sourced. d GSEA plot showing enrichment of transcriptomes of CCT068127 (1 uM) or DMSO treated T24 cells in favorable (top) or unfavorable (bottom) PET identified prognosis genes in bladder cancer. Shown are the FDR value and normalized enrichment score (NES). e Schematic of xenograft experiment and treatment plan with 30 mg/kg daily i.p. injections of CCT068127 or DMSO carrier. p-values compare DMSO-treated to CCT068127-treated groups at each time-point. f Plot showing daily tumor volumes during course of experiment. g Bar plots showing the mass of tumors at explant on day 23. Error bars show mean ± s.d. p-values in a and f are by two-way ANOVA with multiple comparisons adjustment and in g by two-sided Wilcoxon rank sum test. ns non-significant; *p < 0.01; **p < 0.001; ***p < 0.0001. Panel e Created with BioRender.com released under a Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International license (https://creativecommons.org/licenses/by-nc-nd/4.0/deed.en). Source data are provided as a Source Data file.At present, the first line of treatment for bladder cancer consists of chemotherapy regimens that feature cisplatin. As expected, cisplatin also significantly inhibited the growth of bladder cancer cell lines (Fig. 5a). It is worth noting that our approach also predicted that cisplatin desirably normalizes prognostic genes (Supplementary Fig. S7c). To further investigate the impact CCT068127 on bladder cancer cells, we performed RNA-sequencing on T24 bladder cancer cells treated with carrier or CCT068127 (see Methods). 822 genes were found to be up-regulated and 490 genes down-regulated (RPKM > 2, Fold change ≥ ±3 at adjusted p-value ≤ 0.05; Fig. 5c; Supplementary Data 5b). CCT068127 significantly reduced the expression of genes associated with CDK9 inhibition, while showing no significant effect on genes associated with CDK4/6 inhibition (Fig. 5c). To test whether CCT068127 can normalize the predicted prognostic signature, we gathered all the genes linked to favorable or unfavorable prognosis from the prognostic pathways identified by PET and conducted a gene set enrichment analysis. As anticipated from the drug predictions, genes that were upregulated by CCT068127 were significantly enriched in those with PET-predicted favorable prognosis. Conversely, genes that were downregulated by CCT068127 were significantly enriched in in those with PET-predicted unfavorable prognosis (Fig. 5c, d). The up-regulated genes included PRDX4 and DIO3, which were determined by PET to be key markers of a favorable prognosis (Supplementary Fig. S4e). Conversely, the down-regulated genes including XAF1, KPNB1 and PAS1 which were found to be significant markers of an unfavorable prognosis (Fig. 2e). These data corroborated the expectation that the top predicted drug normalized identified prognostic genes.We next tested the in vivo efficacy of CCT068127 in mice. Pharmacokinetic studies identified daily administration of up to 30 mg/kg of CCT068127 for two weeks to be well-tolerated by mice, as evidenced by stable weight during that time (Supplementary Fig. S7e) and lack of gross pathology on necropsy. We then implanted bladder cancer cells subcutaneously in immunodeficient mice and after two weeks, when tumors were approximately 100 mm3 in size, randomly assigned mice to control or treatment groups, injecting them with 30 mg/kg/day of DMSO or CCT068127 i.p. for 10 days, respectively (Fig. 5e). CCT068127-treated mice demonstrated a significant reduction in both tumor volumes and weights (Fig. 5f, g), indicating that CCT068127 was effective at restricting tumor growth in vivo.Taken together, these data indicate that, despite being unexpected from traditional methods, CDK9 inhibition predicated on favorable and unfavorable prognostic biomarkers suppress tumor cell growth both in vitro and in vivo. Furthermore, they corroborate the notion that a PET-based approach can guide repurposing drugs that prove effective in the treatment of cancer.
