Precision Combination Therapies Based on Recurrent Oncogenic Coalterations

Targeted therapies such as HER2 and PARP inhibitors have significantly improved clinical outcomes of patients with cancer (1–3). Despite the transformative advances in sequencing technologies and precision medicine, durable responses to genetically matched targeted therapies have remained limited (4, 5). The limited success can be at least partially attributed to treatment with monotherapies targeting individual oncogenic aberrations. Oncogenesis, however, is driven by co-occurring aberrations that lead to coactivated oncogenic processes. A rich repertoire of preclinical and clinical studies has demonstrated that rational combination therapies blocking multiple oncogenic processes are more effective than monotherapies (6–8). Discovery of combination therapies that may benefit large patient cohorts, however, is a daunting task due to the complexity of the actionable oncogenic aberration landscape (9, 10).

There is growing evidence that targeting the co-occurring oncogenic driver aberrations may induce effective and durable therapeutic responses. For example, coalterations in the PI3K–AKT pathway and BRAF affect response to BRAF inhibitors and establish a rationale for combination therapies in melanoma (11). Similarly, coalterations in KRAS, BRAF, and PIK3CA modulate response to EGFR/HER2-targeted therapies in colorectal cancers (8). Genomic studies focusing on codriver events in KRAS-mutated lung cancers (12) and recent combination therapy trials (e.g., I-PREDICT) focusing on coalterations also strongly support this argument (6). Despite the clinical advances, the precision therapies that can target coalterations remain mostly unexplored. Bioinformatics approaches that extract clinically useful knowledge (e.g., recurrent and actionable coalterations) from patient data will accelerate the implementation of precision combination therapies.

Bioinformatics and machine learning methods have demonstrated potential for significantly improving precision medicine. A recently developed method, SELECT, matches patient transcriptomes to targeted therapies based on drug–gene synthetic lethality interactions inferred from CRISPR/RNAi screens in cell lines (13). The method has demonstrated highly promising outcomes for patient selection in clinical trials mainly for monotherapies. The method DrugCell has generated successful predictions of drug responses in preclinical models and stratification of patients who may benefit from CDK4 and mTOR cotargeting (14). It relies on mutational state and drug chemical structure neural networks trained with responses to hundreds of drugs in thousands of cell lines. The OncoTreat method aims to infer protein activity based on mRNA expression data to prioritize drugs that target aberrant protein activities (15). Logical models (e.g., CellNOpt) predict unseen drug responses based on both signaling database information and drug response data (16). We previously developed network modeling methods that nominate precision combination therapies through analysis of the therapy response mechanisms based on proteomics data (17–19). Despite the significant contributions to automating precision therapy selection, a majority of the approaches suffered from one or a combination of shortcomings such that they (i) benefit from a single data type (e.g., mutational or transcriptional) limiting the coverage of alterations, (ii) focus predominantly on monotherapies, (iii) require comprehensive training data sets, especially for drug responses that are costly to generate in the clinical setting, or (iv) cannot be immediately scaled to a large number of samples due to immense data requirements.

Here, we tested the hypothesis that combination therapies that are tailored to recurrent coalteration signatures lead to significant improvements in clinical and preclinical outcomes. We have developed the machine learning algorithm and precision oncology resource Recurrent Features Leveraged for Combination Therapy (REFLECT). First, the REFLECT method generates multiomic coalteration signatures within a patient cohort (e.g., cohort of EGFR-mutated patients). The REFLECT-generated signatures are a collection of recurrent coalterations, which can be in the form of mutational, copy number, RNA or protein expression alterations. Each REFLECT signature defines a patient subcohort whose members are characterized and discriminated from other subcohorts by the recurrent coalterations. Second, REFLECT maps the recurrent signatures to potential combination therapies based on precision oncology knowledge bases (20, 21). We validated the REFLECT method in vivo, in vitro, and in clinical settings. We have demonstrated use cases in patients with immune-checkpoint therapy (ICT) response markers, DNA repair aberrations, and HER2 activation. We provide an online resource for >200 patient cohorts as well as a software package for customized analyses. We expect that the tool and the resource will guide the discovery of combination therapies and the selection of patients who may benefit from combination therapies.

We have generated a map of the recurrent DNA, mRNA, and protein coalteration signatures across diverse cancer types using the REFLECT pipeline (Fig. 1A–C; Supplementary Fig. S1A and S1B). The analysis benefited from mutation, copy number, mRNA expression, and phosphoproteomic analyses of patient tumor tissues, cell lines, and patient-derived xenografts (PDX). The data cover 33 human cancer types [The Cancer Genome Atlas (TCGA); ref. 22], cell lines from 34 lineages [Cancer Cell Line Encyclopedia (CCLE), Catalogue of Somatic Mutations in Cancer (COSMIC) cell lines project (COSMIC-CLP), MD Anderson Cancer Center (MDACC) cell line project (MCLP); refs. 23–25] and PDXs from six cancer types (NIBR PDXE; ref. 26; Fig. 2A; Supplementary Fig. S2A and S2B; Supplementary Table S1). To correct for batch effects across diverse cancer types, the TCGA protein and mRNA expression data had been integrated and batch normalized with replicate-based normalization and EB⁺⁺ methods (27, 28). We also incorporated data from PDXs and cell lines using a Z-score–based normalization (see Methods). The batch effects were not observed across cancer lineages and sources (patients, cell lines, PDXs) within the compendium of mRNA and protein expression data as assessed through principal component analysis and uniform manifold approximation and projection analysis (Supplementary Fig. S2C and S2D). To account for the differences in exome design, depth coverage, and other batch effects in the sequencing data, we used the multisite, multidisease, consensus pan-cancer MAF file generated by the TCGA Multicenter Mutation Calling in Multiple Cancers (MC3) project (29).

We partitioned the pan-cancer meta-cohort into patient cohorts, each of which is defined by the presence of a “master” biomarker that is shared by all members of the cohort (Fig. 1B; Supplementary Table S2). The master biomarkers, which stratify the meta-cohort into cohorts prior to REFLECT analyses, are potential therapeutically actionable genomic alterations (e.g., HER2 amplification). The therapeutic actionability of an alteration is defined in reference to the MDACC Precision Oncology Knowledge Base (N = 182; Supplementary Table S3) and NCI-Molecular Analysis for Therapy Choice (NCI-MATCH) precision therapy trials (N = 15; Supplementary Table S4; refs. 20, 21, 30). We have defined and added an immunotherapy cohort based on CD274 (PD-L1) overexpression and immune cell infiltration (Supplementary Table S5). Additional cohorts were constructed based on the presence of key oncogenic alterations, MYC amplifications, and TP53 and KRAS mutations (Supplementary Table S5). In total, we generated 201 master biomarker-based cohorts, which form the core of the REFLECT resource, and yet can be expanded with new cohorts (Figs. 1B and 2B). According to OncoKB, at least 25 of the master markers are recognized by the FDA and another 14 are recommended by the National Comprehensive Cancer Network as predictive of response to targeted therapy in the clinic. The remaining markers are under evaluation in precision oncology trials (20, 21).

For each of these 201 cohorts, each defined by a master marker, we generated the multiomic/multigene signatures, termed REFLECT signatures. In order to capture the coalteration signatures relevant to cancer and therapeutic targets, we focused on the genes listed within the COSMIC Cancer Gene Census (CGC; ref. 31) for the genomic and transcriptomic alterations. For signature detection, we used a sparse hierarchical clustering algorithm (32) followed by statistical assessment of recurrent alterations (Fig. 1C). The sparse hierarchical clustering algorithm iterates LASSO-based feature detection and hierarchical clustering until convergence to a final signature (see Methods). In the succession of iterations, the feature detection selects a small set of alterations that contribute most to discriminating patient subcohorts, which are identified by the clustering. Next, the significantly recurrent events from the signature are extracted with a statistical test (see Supplementary Methods). The statistical test filters out the features whose observed recurrence is likely by chance in a large patient cohort. The test preserves the ones that are significantly recurrent and thus potentially functional. The resulting signature carries a combination of the master biomarker (i.e., the identifier of the associated cohort) and a small set of REFLECT-detected recurrent features with nonzero contribution (W ) to patient classification. The distinct signatures discriminate subcohorts within a patient cohort. To demonstrate the robustness of REFLECT signatures, we quantified the concordance between the features that were calculated with complete versus subsampled data with varying coverage (60%–95%; see Supplementary Methods; Supplementary Fig. S2E). The concordance between partial and complete data sets had an average of 83% for 70% coverage across the seven cohorts tested. The concordance increased monotonically with increased coverage. In summary, we have acquired a compendium of robust multiomic REFLECT signatures spanning patient cohorts and subcohort-specific coalterations selected by our algorithm.

Based on the REFLECT signatures, we built a comprehensive resource of recurrent coalterations that can guide selection of combination therapies tailored to complex molecular signatures (Supplementary Fig. S2F and S2G). For this purpose, we extracted recurrent and coactionable alterations from REFLECT signatures of 201 patient cohorts. Next, we identified the patient subcohorts that carry potentially coactionable events based on the recurrent coalterations. The coactionable targets include (i) the master biomarker that defines the cohort and/or (ii) REFLECT-selected recurrent and therapeutically actionable coalterations.

To select actionable mutation and copy-number events in subcohorts, we matched the oncogenic coalterations in the signatures to the MDACC and OncoKB precision oncology knowledge bases. To match protein and mRNA expression aberrations to therapeutic agents, we developed a preclinical resource based on the drug-target information from the MDACC Precision Oncology Knowledge Base, OncoKB, Genomics of Drug Sensitivity in Cancer (GDSC), and Cancer Therapeutics Response Portal (CTRP) drug-target annotations (see Supplementary Methods; Supplementary Tables S6 and S7; refs. 20, 21, 23, 24). The mRNA/proteomics-based actionable events are selected to match drugs to expression aberration of their prospective targets. For example, the increased phosphorylation and total levels of signaling molecules in key pathways [e.g., AKT/PI3K, receptor tyrosine kinase (RTK), MAPK, cell cycle, DNA repair] are linked to their well-characterized inhibitors. Our resulting proteomics- and transcriptomics-based resource matches 82 proteomic and 89 mRNA aberrations to targeted therapy options (Supplementary Tables S6 and S7).

We quantified the number of REFLECT-matched precision combination therapies and recurrent coalteration signatures (P-recurrence < 0.05). Each signature-matched combination involves at least one FDA-approved agent for the treatment of one or more cancer (sub)types. The selected drug combinations target a diverse battery of potentially actionable multiomic alterations such that no individual event dominates the multiomic target space (Supplementary Fig. S2H). As expected, the most frequently altered oncogenes (e.g., PIK3CA) were among the most common targets (e.g., targeting the PIK3CA alterations is included in 6% of all combinations). The mutation and copy-number coalteration signature analysis led to 325 unique drug combinations. Based on mutational and copy-number signatures, 18% of the patients in the meta-cohort are matched to combination therapies. For protein and mRNA aberrations (|Z-expression| > 2, P-recurrence < 0.05), 1,235 and 1,087 unique signatures were matched to FDA-approved drugs, respectively. Based on protein and mRNA analyses, 32% and 29% of patients were matched to at least one combination therapy, respectively. A total of 2,166 combinations with one or more FDA-approved agents were matched to mutation/copy number, mRNA, and protein coalteration signatures. In total, 45% of the patients in the meta-cohort were matched to at least one combination therapy based on the multiomic data.

The majority of the REFLECT-selected coalterations involved different genes such that therapeutic strategies to target the coalterations would likely involve drug combinations. We have detected only 50 coalterations that converged on the same gene (i.e., a gene in a REFLECT signature is identical to the corresponding stratification marker). Given that there are thousands of actionable coalterations in the REFLECT resource, such events can be considered rare. Meanwhile, within the same gene coalterations, we observed an enrichment of events involving concurrent amplification and overexpression of the same gene (60% of the same gene coalteration cases). In the remaining cases, the mutational alterations usually lead to increased phosphorylation and likely hyperactivation of the corresponding protein products (e.g., mutations in EGFR, BRAF, and AKT that led to increased protein phosphorylation). The exceptions were increased KIT or EGFR total levels in tumors that carry KIT or EGFR mutations. We also observed that, in the copy-number analysis, nearly 60% of the gene amplification and deletion events were mapped to the same chromosomal regions and resided in the identical Genomic Identification of Significant Targets in Cancer (GISTIC) regions. These observations suggest arm-level copy-number alterations (CNA) play a dominant role in defining the coalterations space; however, coalteration of focal chromosomal events residing on distinct locations is still common (40%).

Next, we addressed whether the REFLECT-generated pan-cancer signatures are solely driven by cancer types/cell of origin or tumor proliferation levels as opposed to the landscape of multiomic alterations (Supplementary Fig. S2I). To assess tumor proliferation levels, we calculated the proliferation scores based on the expression levels of key proliferation genes reported previously (ref. 33; see Methods). We observed a weak correlation (Spearman correlation = 0.14) between clustering patterns (i.e., ranking of samples across dendrograms) and proliferation scores. The median concordance between cancer types across the REFLECT dendrograms was also moderate across all REFLECT analyses (N = 201 stratifications) with a concordance score of 0.39 (range, 0–1; see Methods). More importantly, the distribution of cancer-type concordance scores across all REFLECT analyses follows a bimodal distribution such that for a large fraction of stratified cohorts, the median concordance was low (0.25), whereas a smaller population of cohorts had a higher concordance (median = 0.6). The results suggest that in a larger fraction of the cohorts, clustering patterns were driven by the distribution of the molecular features that are affected by factors other than tumor type; in a smaller fraction of cohorts, tumor type is critical but likely not the only factor determining the clustering patterns. The analyses suggest the REFLECT signatures are not dominated by cell of origin or proliferation levels of the tumors. In summary, the outcome of the REFLECT analyses is a comprehensive resource of precision combination therapies tailored to recurrent coalteration signatures (Fig. 2C and REFLECT portal for coalteration maps).

We then tested whether REFLECT-matched drug combinations can induce significant therapeutic benefit in a PDX preclinical trial. The preclinical trial recapitulates a clinical trial design and captures the interpatient heterogeneity through the “one animal per model per treatment” approach (1 × 1 × 1; ref. 26). We analyzed the therapeutic effects including therapy efficacy, survival, and synergy of the drug combinations in PDXs. The validations compare the therapeutic effects of REFLECT-matched combination therapies versus the combinations that are not matched by REFLECT.

The in vivo data are available for a large array of PDXs from NIBR PDXE (277 PDXs from six cancer types treated with 22 targeted drug combinations in 1,168 independent treatments; ref. 26). The PDXs are profiled for mutation, CNA, and mRNA expression levels (Supplementary Fig. S2A and S2B). The accuracy, robustness, and clinical relevance of the resource had been established through replication of the drug response–genotype associations, mechanisms of resistance from clinical studies, and >2,000 single-drug response data points from historical PDX experiments (26). The established quality and clinical relevance of the screen suggest that the data are well suited for validation of the REFLECT pipeline.

In the validations, we used a merged data set of patient and PDX samples to execute the REFLECT pipeline and identify the predictions in PDXs. The mRNA data are normalized with the Z-score–based normalization (see Methods). We calculated the multiomic REFLECT signatures and matched drug combinations to recurrent coalteration signatures (Fig. 3A; Supplementary Fig. S3A; Supplementary Table S8). FDA-approved drug combinations were REFLECT-matched to 36% of the PDX models (100 of 277 PDXs). Fifteen of the REFLECT-predicted combinations had been experimentally tested in matched PDXs within the NIBR PDXE resource (Supplementary Table S9). We compared the therapeutic benefit by REFLECT-matched drug combinations [REFLECT(+)] to unmatched combinations [REFLECT(−); Fig. 3B–D]. The median tumor volume change with respect to the baseline was a 35.4% decrease for REFLECT(+) and a 5.1% increase for REFLECT(−) cases (two-tailed t test, P = 0.017; Fig. 3B). The progression-free survival (PFS), defined as the amount of time that PDXs take to double their volume from the baseline (26), of REFLECT(+) PDXs was significantly longer than that of REFLECT(−) PDXs (Fig. 3C for Kaplan–Meier curves and survival statistics). Finally, we quantified the relative drug synergy using the highest single agent (HSA) model, which provides a simple and unbiased framework to compare drug synergy differences between REFLECT(+) and REFLECT(−) cases. The HSA model quantifies the improvement in efficacy of a drug combination over the agent that is a part of the combination and carries the highest impact as a monoagent. The drug combinations in the REFLECT(+) cohort had a significantly higher synergy (median HSA score: 28.02) compared with the REFLECT(−) cohort (median HSA score: −0.24; two-tailed t test P = 0.005; Fig. 3D). There is a significant enrichment of more synergistic interactions in the REFLECT(+) cohort such that 75% of the REFLECT(+) and 27% of the REFLECT(−) combinations had HSA scores greater than 20 (Fisher exact test P = 0.003; Supplementary Fig. S3B). On the other hand, we observed a shift to more antagonistic interactions in the REFLECT(−) cohort such that 30% of the REFLECT(−) combinations had HSA scores less than −20, in contrast to only 8% in the REFLECT(+) group.

We have also tested an alternative validation strategy, which involved matching of PDXs to drug combinations using the precomputed REFLECT signatures based on the TCGA samples only (Fig. 2C) and without rerunning the REFLECT with the merged data set (Supplementary Fig. S3C). Compared with the approach that involved rerunning REFLECT (Fig. 3A), this strategy led to a larger volume for REFLECT(+) cases (N = 237) as it lacks the recurrence constraint for the alterations in the PDXs. With this alternative approach, the therapeutic improvement with respect to the REFLECT(−) group is still significant (Supplementary Fig. S3C). As the matching is not constrained by the recurrence criterion, the overall therapeutic gain, however, is moderate compared with the integrated REFLECT analysis with the merged data (Fig. 3A–D). In conclusion, the validations in PDX models demonstrated that REFLECT introduced significant therapeutic benefits with improved tumor shrinkage, PFS, and drug synergy.

We devised an independent validation scheme based on three comprehensive drug combination screens in cancer cell lines (AstraZeneca, NCI ALMANAC, and O'Neil et al. resources; refs. 34–36), which had been previously merged and used within the DREAM challenges (Supplementary Fig. S3D; ref. 34). The merged screen includes 75,612 unique data points for 3,694 combination targets across 142 cell lines. The mutation, copy number, mRNA expression, and phosphoproteomic data are available from the COSMIC-CLP, CCLE, and MCLP projects (23–25). The drug synergy (Loewe) scores and target annotations were already incorporated within the framework of DREAM challenges (34). The reproducibility and data quality had been tested through comparisons of replicates across the screens and assessment of curve-fitting for drug responses using the Combenefit analysis suite (37). The data set has already been evaluated and analyzed by >160 research teams for benchmarking drug response prediction algorithms (34). The well-tested drug response resource serves as an unbiased platform for testing REFLECT predictions with high statistical power.

First, we identified the subset of drug combinations that target the REFLECT-nominated coalterations in matched cell lines [REFLECT(+); Fig. 3E; Supplementary Fig. S3D; Supplementary Table S10]. We compared the therapeutic effect as quantified by the synergy scores between the REFLECT(+) combinations versus all other drug combinations (Fig. 3F and G; Supplementary Table S11). The synergy scores for REFLECT(+) combinations (mean = 25.8, median = 30) are significantly higher compared with scores for all other combinations (mean = 10.5, median = 7.3; two-tailed t test, P < 0.001; Fig. 3F). For each individual data modality (protein, mRNA, mutation/CNA), REFLECT(+) synergy scores are also significantly higher than the scores for the nonmatched cases (Fig. 3G; Supplementary Fig. S3E–S3G; Supplementary Table S11). The therapeutic benefit monotonically increased with recurrence (low P-recurrence) and expression levels (high Z-scores). The results from the systematic validation scheme suggest that REFLECT predictions may guide selection of more synergistic drug combinations, regardless of the data modality. Meanwhile, recurrence of coalterations and protein/gene overexpression are strong predictors of synergistic drug responses.

We tested REFLECT in the clinical setting using the patient response and molecular alteration data from the I-PREDICT trial (NCT02534675; Fig. 3H; Supplementary Fig. S3H; ref. 6). The trial included patients with metastatic cancers (stage IV disease) treated with targeted therapy, in most cases in combinations. The cohort included patients with diverse cancer types, but predominantly having gastrointestinal, gynecologic, or breast cancers. The drug combinations were matched to patients based on next-generation sequencing and immunotherapy response markers (PD-L1 IHC and tumor mutation burden) of their malignancies. A tumor board consisting of physicians and scientists matched drug combinations to the response markers. We had access to therapeutically actionable alteration annotations, immunotherapy markers, and survival data of each patient (6). As our analyses focused on alterations in multitargets and combination therapies, we excluded monotherapies and combination therapies targeting the same targets (e.g., trastuzumab and pertuzumab targeting ERBB2/3 alterations). We identified patients whose tumors carried REFLECT-determined coalterations and were treated with the REFLECT-matched combination therapies, the REFLECT(+) cohort (N = 12), and patients who did not match to a prediction, the REFLECT(−) cohort (N = 17; Fig. 3H; Supplementary Fig. S3H; see Supplementary Table S12 for the patient characteristics). The treatments in the REFLECT(+) and REFLECT(−) groups were highly diverse with minimal overlap. In total, there were 27 unique treatments across 29 patients, with only one treatment administered to two patients in the REFLECT(−) group (trametinib and bevacizumab), and there was only one overlapping combination (trametinib and everolimus) between the two groups. Although the drug combinations in both groups involved a total of 25 unique agents, only eight monoagents (32% of all agents) were matched in both REFLECT(−) and REFLECT(+) combinations. The remaining 17 monoagents were administered as part of either REFLECT(+) or REFLECT(−) combinations but not both. The REFLECT(+) cohort had significantly longer PFS and overall survival than the REFLECT(−) cohort (Fig. 3I for Kaplan–Meier curves and statistics). There was no observable enrichment of a particular cancer type in either group, suggesting that disease type was not confounding the predictions (Supplementary Fig. S3I). Through a multivariate Cox regression analysis, we validated that the REFLECT status was the only statistically significant predictor of survival against a set of potential confounding factors (PD-L1 expression, tumor mutation burden, and BRCA1 mutation; Supplementary Fig. S3J and S3K). Microsatellite instability (MSI) status most likely had no impact on the survival predictions, as only one patient carried an MSI-high annotation. Thus, the data from the I-PREDICT trial provided strong evidence that REFLECT can provide significant therapeutic benefits in patients with cancer.

Precision combination therapies that improve responses to immunotherapies (e.g., PD-1/PD-L1 inhibitors) are highly desired. Here, we determined the multiomic signatures and co-occurring actionable targets within cohorts of likely responders to immunotherapies. First, we focused on PD-L1 expression as it has been established as a potential marker of PD-1/PD-L1 inhibitor response in select cancers. In addition to the PD-L1 expression, another response predictor is the immune cell infiltration into the tumor microenvironment (TME; ref. 38). To establish a robust PD-L1–high cohort, we included the tumors whose PD-L1 mRNA and PD-L1 protein levels are in the top quartile across the pan-cancer meta-cohort. Next, we selected the PD-L1–high tumors whose tumor-infiltrating leukocyte (TIL) fraction is in the top quartile across the pan-cancer meta-cohort (Fig. 4A and B). The TIL fraction had already been quantified with a mixture model of tumor subpopulations based on DNA methylation markers with the greatest differences between pure leukocytes and nonleukocyte tissues (39). The resulting high PD-L1 and high TIL cohort contains 406 cases from 23 disease types (Fig. 4C; Supplementary Table S5).

Using REFLECT, we identified the significantly recurrent mutational, mRNA, and phosphoproteomic signatures that co-occur with high PD-L1/TIL (Fig. 4D; Supplementary Fig. S4A). The mutation analysis of the PD-L1/TIL cohort identified a signature of eight genes (P-recurrence < 0.05) covering the commonly mutated oncogenes in cancer. The most recurrent signature genes (P-recurrence < 0.0001) discriminating patient subcohorts are TP53 (50% of the cohort), PIK3CA (15%), and CDKN2A (10%) followed by RB1, BRAF, and KRAS. The mRNA REFLECT signature is characterized by a combination of significantly recurrent events (P-recurrence < 0.05) in immune genes as well as DNA damage response and repair markers. The strongest contributing mRNA feature is the overexpression of S100A7 (P-recurrence < 0.0001). In the context of cancer, the S100A7 functions in the TME with roles in inflammation-driven carcinogenesis and recruitment of immune cells, including tumor-associated macrophages. S100A7 is highly expressed in a series of cancers including head and neck cancers (40–42). The S100A7-enriched subcohorts also carried increased expression of TP63, HRAS, and FGFR3, particularly in head and neck, bladder, and squamous lung cancers.

The transcriptomic signature is also marked by genes related to DNA repair through chromatin remodeling and epigenetic control (SETD2, NSD1, TRRAP, CHD2, ZMYM3, and BRD4), G₂–M DNA damage checkpoint (ATM and ATR), and nucleotide excision repair (XPC and ERCC2). We observed two patient subcohorts with enrichment of epigenetic regulators of DNA damage response (DDR) and repair. The first subcohort, composed of head and neck cancer cases, is characterized by low expression of the DNA damage checkpoint kinase ATM, as well as epigenetic regulators SETD2 and ZMYM3 genes, which function in recruitment of key enzymes (e.g., BRCA1 and RAD51) in both homologous recombination (HR) and mismatch repair (43). A subgroup of this patient group is also characterized by the concurrent loss of nucleotide excision repair (NER) enzyme XPC and homologous recombination repair (HRR) protein XRCC2. The second subcohort was composed of patients with diverse cancer types (breast, stomach, lung, and colorectal cancers) and characterized by high expression of epigenetic regulators NSD1 (encoding SETD1 protein) and TRRAP, which links chromatin remodeling to DNA double-strand break (DSB) repair (44). A functional enrichment analysis using the Web-based Gene Set Analysis Toolkit (45) revealed highly significant enrichment of gene sets involved in histone methylation (FDR P < 0.001), DNA repair–associated processes, replicative senescence (FDR P = 0.03), and telomere maintenance (FDR P = 0.01; Fig. 4E). Another subcohort populated by lung squamous cancer cases was characterized by coexpression of the immune chemotactic receptor CCR4 and TNFRSF17, a marker of B cells (46, 47). We demonstrated, based on a CIBERSORT analysis, that tumors in the high CCR4/TNFRSF17 coexpression subcohort are significantly enriched with naïve and plasma B-cell markers (Fig. 4F).

In the phosphoproteomics analysis, depletion of RAD50, a key protein for initial processing of DSB prior to HRR and tethering the two broken ends of damaged DNA, has the highest contribution to discriminating the PD-L1/TIL subcohorts (Fig. 4D; Supplementary Fig. S4A). The RAD50 protein loss was recurrent in two patient subcohorts: one is populated by diverse cancer types (N = 111, P-recurrence < 0.0001), whereas the other is enriched by bladder cancers (N = 48, P-recurrence < 0.0001). Increased CHK1 protein level (P-recurrence < 0.0001, N = 25) and CHK2 phosphorylation at T68 (P-recurrence < 0.0001, N = 52), a marker of DDR and G₂–M DNA damage checkpoint, also co-occurred within a fraction of the RAD50-depleted subcohort, suggesting aberrant HR-based DSB repair in such tumors.

We increased the scope and depth of the immune analysis to identify potentially actionable and recurrent alterations that are robustly associated with diverse immunotherapy markers. This effort is based on the rationale that actionable and recurrent alterations that consistently overlap with diverse immune markers may inform potentially effective combination therapies. In addition to the PD-L1/TIL cohort, we defined 10 patient cohorts marked with combinatorial enrichment of potential immunotherapy response markers (tumor immune infiltration, tumor mutation burden, T-cell inflammatory signature), major immune cell types (CD8⁺ T cells, CD4⁺ T cells, B cells, and Tregs), and expression of targets of immune-checkpoint therapies that are in clinical use or trials (CTLA4, ICOS, TIGIT, and LAG3; Supplementary Fig. S4B; refs. 48–55). The proportion of immune cell types in each patient sample had been inferred by a CIBERSORT analysis of transcriptomics data and imported from Thorsson and colleagues (39). The cohorts are generated through the inclusion of tumors whose levels of corresponding markers are in top quartiles across the pan-cancer meta-cohort as in the case of the PD-L1/TIL cohort (Fig. 4A–C). To increase the depth of the REFLECT analysis beyond the COSMIC CGC, we expanded our gene set by the addition of the genes in the IFNγ gene signature, the expanded immune gene signature and the major histocompatibility complex (MHC)—the gene sets that are associated with immunotherapy responses (Supplementary Fig. S4C; ref. 51). We generated the REFLECT signatures and analyzed the mRNA (based on the expanded gene list) and protein signatures across the 11 immune cohorts.

The transcriptomic signatures within the PD-L1/TIL cohort remained enriched by mediators of DDR when we added the immune-related genes to the COSMIC CGC set (Fig. 4D; Supplementary Fig. S4D). A gene set enrichment returned DDR and histone methylation among the top significantly enriched processes (Supplementary Fig. S4E). The mRNA coalteration signatures, however, were also enriched by a set of inflammatory features, most striking being IFNG expression (Supplementary Fig. S4D). In the PD-L1/TIL cohort, we observed a highly recurrent and strong IFNG expression signature that spans diverse tumor types (melanoma, sarcoma, bladder, breast, and cervical cancers). The IFNG overexpression was mutually exclusive with expression of the inflammatory protein S100A7 (mainly in head and neck cancers; Fisher exact test, P < 1.0e−6), whereas the patients whose tumors express the LAG3 immune-checkpoint receptor were overlapping with both IFNG-enriched and S100A7-enriched subcohorts. LAG3 and S100A7 coexpression also overlaps with expression loss of a set of DNA repair genes such as DNA damage checkpoint kinases (ATR and ATM), telomere maintenance gene POT1, epigenetic regulator ZMYM3, and the base excision repair gene XPC, suggesting DNA repair and inflammatory events coexist in this subcohort. Next, we performed an integrated analysis of the discriminant RNA alterations that were shared across the different immune cohorts (Supplementary Fig. S4F). This analysis generated a consensus signature across the cohorts each marked by a different set of immune markers (immune-checkpoint expression, cell type, and mutation burden). The consensus signature was highly enriched by the recurrent inflammatory alterations involving the overexpression of genes encoding IFNG, CCR7, and CCR4, B-cell lymphocyte lineage markers PAX5 and CD79A–CD79B (B-cell receptor complex), costimulatory immune receptor CD28, and MHC complex.

We calculated the consensus proteomic alteration signatures by quantifying the discriminating features that were shared across different immunotherapy response cohorts (Supplementary Fig. S4G). Consistent with the coalteration landscape of the PD-L1/TIL cohort, the REFLECT analyses showed that the recurrent loss of the HR repair protein RAD50 was the most discriminant proteomic alteration in nearly all immune cohorts (Supplementary Fig. S4G). Other DNA repair proteins, single-strand break repair protein XRCC1, G₂–M checkpoint kinase pCHK2, and epigenetic regulator of DNA repair SETD2, were among the most recurrent coalterations with the immune markers. Finally, although not as strong as the DNA repair and inflammatory signals, we observed signatures associated with activation of the oncogenic growth pathways (e.g., RAF1, pEGFR, mTOR, AKT, and cMET) in the proteomic and transcriptomic modalities (Supplementary Fig. S4F and S4G).

In summary, we have consistently observed coalterations in DNA repair, particularly the HR repair mechanisms (e.g., loss of RAD50), DNA checkpoint enzymes (CHK1/2, ATM), and epigenetic regulators of DNA repair (e.g., SETD2) in combination with enrichment of diverse inflammation markers. The coalteration signatures involving DNA repair, oncogenic signaling and inflammation signatures warrant further studies to investigate therapeutic benefits first in preclinical models and next in patients with cancer, whose tumors carry the therapy-matched coalterations.

PARP inhibitors have provided significant therapeutic benefit, particularly for treatment of tumors with DNA repair defects. Here, we have nominated combination therapies in which PARP and other DNA repair inhibitors serve as the anchor compound. We applied REFLECT to patient cohorts marked with alterations in 22 DNA repair genes (Supplementary Table S13). The selected genes are implicated in BRCAness and potentially synthetically lethal with inhibition of PARP (or other DNA repair enzymes), as listed in the MDACC Precision Oncology Knowledge Base (20). The DNA repair genes were selected based on the annotations in MDACC and MSKCC precision oncology knowledge bases with refinement of the BRCAness and PARP inhibitor sensitivity information through literature curation (56). All of the alterations involved genes that function in DNA DSB repair with the exception of the mismatch-repair gene MSH3. A majority of the genes play direct roles in HRR. For example, the genes encoding members of the MRN complex (i.e., RAD50 and MRE11) tether DNA double strands prior to HRR. BRCA partners such as BARD1 and PALB2 are included. The gene encoding the key enzyme of the HRR mechanism RAD51, which catalyzes the strand transfer, as well as various RAD51 partners (e.g., RAD51B/D, RAD52, RAD54, and XRCC2), are included. The list also includes the genes encoding the G₂–M DNA damage checkpoint enzymes (ATM, ATR, and CHK1/2). Despite being a mismatch-repair gene, MSH3 was included in the list as BRCA1 is a downstream effector of the MSH3, and it was included as a BRCAness gene that is potentially associated with PARP inhibitor sensitivity (56). First, we generated the multiomic signatures for each of the cohorts to identify recurrent coalterations. Next, we selected potential combination therapies using the forward target selection criterion (Supplementary Fig. S1B; Supplementary Table S14).

To achieve an integrated analysis of the signatures associated with the DNA repair alterations, we mapped the proteomic coalterations paired with cohort markers on a bipartite graph (Fig. 5A). On the graph, each recurrent feature within a cohort is linked to the respective co-occurring master biomarker of the cohort (i.e., the DNA repair aberration that defines a cohort; Fig. 5A). We observed an enrichment of aberrant expression in DNA repair proteins (RAD51, ATM, and CHK1/2), RTKs (pEGFR, pHER2, MET, KIT, and VEGFR), and downstream signaling proteins (pAKT, mTOR, p70S6K, pJNK, pSRC, pP38/MAPK14, pPDK1, pRAF, and pMAPK) that co-occur with the master biomarkers in the DNA repair family. Guided by the graph enrichment of coalterations in DNA damage markers, RTKs, and downstream signaling, we experimentally cotargeted the PARP and co-occurring aberrations (Fig. 5B–D; Supplementary Fig. S5A–S5C). In the drug perturbation experiments, we used the cell lines molecularly matched to (i.e., cocluster with) the specific patient subcohorts carrying the coalterations. We prioritized molecularly matched gynecologic and breast cancer patient subcohorts because PARP targeting has demonstrated strong, albeit transient, therapeutic benefit for some patients with high-grade serous ovarian or basal-like breast cancer. Using the forward target selection (Supplementary Fig. S1B) in the RAD50-, PALB2-, and ATM-altered cohorts, we identified subcohorts characterized by statistically significant recurrent CHK2 phosphorylation, AKT phosphorylation, and MET expression, respectively (P-recurrence < 0.25). For each of these subcohorts, we identified a cell line with a matched co-occurrence pattern, namely MFE-296 (RAD50 mutation with recurrent CHK2 phosphorylation), HEC50B (PALB2 deletion with recurrent AKT phosphorylation), and DU4475 (ATM deletion with recurrent MET protein overexpression). The drug combinations produce additive to synergistic effects in each of the tests in the matched cell lines, as quantified by the combination index (Fig. 5B–D). We have also tested the same drug combinations in unmatched cell lines that do not carry the corresponding coalterations (i.e., REFLECT criterion for unmatched: P-recurrence > 0.25; Supplementary Fig. S5A–S5C). We did not observe therapeutic gains by the combinations in the negative control cases. Overall, the REFLECT calculations determined the landscape of coalterations involving DNA repair aberrations. When coupled with preclinical and clinical translation, the landscape may inform selection of effective combination therapies tailored to coalterations.

Therapeutic targeting of HER2 has improved clinical outcomes particularly in HER2⁺ breast cancers, yet resistance to therapy remains a major clinical problem. Drug combinations that will improve the therapy response rates among HER2⁺ tumors are highly needed. We used the REFLECT pipeline and a reverse search strategy (Supplementary Fig. S1B) to identify patient cohorts that may benefit from combination therapies involving HER2 inhibitors. In the reverse search strategy, we scanned proteomic signatures within 201 cohorts, each defined by the actionable biomarkers. We identified patient subcohorts carrying actionable oncogenic aberrations (master biomarkers) and HER2 activation markers [increased total protein or phosphorylation (pY1248) levels]. The analysis led to the identification of 69 subcohorts that carry a HER2 activation marker (P-recurrence < 0.05) alongside an actionable master biomarker (Fig. 6A). As a result, 338 patients (4% of the TCGA meta-cohort) were matched to 51 unique combinations involving HER2 targeting (Supplementary Table S15). For proof-of-principle validation experiments, we identified cell lines that map to subcohorts of increased HER2 activity. The molecularly matched cell lines were sampled within the cohorts of (i) AURKA mutation/amplification (cell line: BT474), (ii) MAP3K4 mutation/amplification (cell line: SKOV3), (iii) CDKN1B mutation/loss (cell line: HCC1954), and (iv) PTEN loss/mutation (cell line: MDAMB468; Fig. 6B–E; Supplementary Fig. S6A–S6D). These cell lines were perturbed with HER2 inhibitors and agents targeting the cohort biomarkers. The drug combinations targeting the coalterations substantially improved responses compared with single agents. The drug combinations carried moderate to strong synergies as quantified by the combination index (Fig. 6B–E). The same drug combinations involving HER2 inhibitors did not lead to significant therapeutic gains in the unmatched cell lines that serve as negative controls (REFLECT criterion for unmatched: P-recurrence > 0.25; Supplementary Fig. S6A–S6D). More importantly, we have mapped the coalteration signatures that may guide selection of patients with HER2⁺ tumors and match to 51 unique drug combinations involving HER2 inhibitors. The REFLECT calculations generated a precision combination therapy resource for the treatment of HER2-activated tumors through targeting coalterations.

Here, we present a resource and a bioinformatic analysis method, REFLECT, to match combination therapies to combinatorial aberration signatures in patients’ tumors. The REFLECT method includes a machine learning algorithm for the detection of recurrent coalterations and an informatics component to match the coalterations to therapies. With a REFLECT-based analysis, we have generated a resource that contains combination therapies matched to >10,000 patients across >200 patient cohorts. Our comprehensive validations using cell lines, PDXs, and clinical data demonstrate that REFLECT consistently provides significant therapeutic benefits in both preclinical and clinical settings. The REFLECT predictions led to more desired drug responses in terms of efficacy, synergy, and survival outcome compared with control groups for all data sets (preclinical and clinical) and modalities (DNA, mRNA, and protein). This is particularly striking given that the predictions rely on signatures in therapy-naïve patients or preclinical models and are not constrained by any drug response data. The ability to nominate drug combinations based on data from pretreatment tissues is critical, as it suggests REFLECT can be easily implemented to currently used or emerging precision oncology protocols. Based on the consistently superior outcomes from REFLECT(+) cases, we conclude that REFLECT has the potential to facilitate precision medicine applications.

REFLECT is built on a set of core principles in cancer genomics and therapeutics. First, oncogenesis and cancer progression are driven by multiple coexisting driver aberrations (12, 57, 58). Such coexisting aberrations may be causally coupled to each other, and targeting such coupled events may lead to significantly synergistic, strong responses. Alternatively, coalterations may drive parallel and relatively independent events without causal links in the same or separate tumor subclones. Such independent events may contribute to cancer progression in an additive manner. Subsequently, targeting the independent events may lead to an additive response with increased durability, a highly desired therapeutic outcome. Second, we emphasized recurrence in the analysis of coalteration signatures. Oncogenic alterations that are significantly recurrent across patient cohorts are likely functional, as such recurrence across samples suggests selection to drive oncogenic hallmarks (59). The significant recurrence of actionable features also suggests the resulting therapy regimes may benefit substantially large patient populations. Third, most precision therapies are usually matched to functional variants of an oncogene (e.g., BRAF^V600E). Here, when selecting potentially actionable mutations, we have not filtered out variants of unknown significance (VUS); instead we have included all variants of a master biomarker if the variant exists in more than two cases in the meta-cohort. This approach involves the risk that some recurring passenger mutations are included in the analysis. This approach, however, has the advantage that we detect the REFLECT signatures that are shared by both known functional/actionable variants and VUS, which may guide therapeutic options to target tumors with important VUS. Finally, the stratification into patient cohorts based on actionable biomarkers emulates the ongoing basket trials such as NCI-MATCH and NCI-ComboMATCH. The approach also carries the added value that it enables matching combination therapies to subcohorts with shared molecular compositions. In summary, the success of our approach relies on the observation that recurrent coalterations that coactivate multiple oncogenic processes may prevent benefit from monotherapies, yet such coalterations that are shared across patient subcohorts create vulnerabilities to matched combination therapies. The REFLECT resource is highly dynamic and is open to exponential growth with the addition of new cohorts. We anticipate further growth and increased usability of the REFLECT resource as new drug targets are discovered and translated through clinical trials.

To demonstrate the capability of REFLECT, we focused on the signatures of three clinically relevant events: potential immunotherapy response markers, DNA repair/damage response aberrations, and HER2 activation. The analyses of HER2 and DNA repair aberrations have identified patient subcohorts matched to the combination therapy candidates that we validate in matched cell lines. The experimental validations reveal either novel drug synergies, such as cotargeting Aurora kinase and HER2 in characterized patient groups, or patients who potentially would benefit from previously tested drug combinations, such as those cotargeting PARP with CHK2 (60), AKT (61, 62), or MET (63) and cotargeting HER2 with CDK4 (64) or PI3K (65). Indeed, a large number of synergistic combination therapies have been previously discovered in diverse preclinical studies. It remains elusive, however, if there is a substantial number of patients with cancer with a well-defined signature who may benefit from a previously discovered combination. Therefore, a key contribution of REFLECT to future precision medicine applications is in the mapping of recurrent coalteration signatures to known drug combinations, as well as discovery of novel combination therapies tailored to recurrent coalteration signatures of patients with cancer.

Selection of ICT response markers particularly in a pan-cancer fashion is not a trivial task. Despite their success in select cancers, neither PD-L1 expression nor tumor mutation burden is a robust predictive marker across all cancer contexts (66, 67). Acknowledging these shortcomings of individual biomarkers, we used a compendium of events involving immune-checkpoint expression, tumor immune cell infiltration, and immune cell–type proportions to define a set of cohorts that are most likely to be primed for ICT response. The resulting REFLECT-generated signatures of potential response markers involve enrichment of specific inflammation (e.g., IFNG, S100A7, and LAG3) and immune cell (B cells) markers as well as DNA repair alterations, particularly the HR repair and G₂–M DNA damage checkpoint mechanisms. The coenrichment of DNA damage and immunotherapy response markers establishes a rationale for testing PD1/PD-L1 inhibitors (and likely other ICTs) in combination with DNA repair inhibitors in molecularly matched patients. Clinical trials combining PARP inhibitors and immunotherapies are already in progress with promising results in breast, prostate, and ovarian cancers (10, 68). However, the therapeutic opportunities are not limited to PARP inhibitors, as combinations of ATR and PD-L1 inhibitors have recently progressed to clinical trials (69). The subcohort-specific co-occurrence patterns (e.g., recurrent RAD50 loss) may guide the selection of patients who are likely to benefit and inform which specific DDR-targeting agent may optimally combine with anti–PD-1/PD-L1 for a specific indication. Associated with the immune markers, we have also identified a series of coalterations in epigenetic factors (SETD2, NSD1, TRRAP, CHD2, and ZMYM3) that link chromatin remodeling to DNA repair through loading of key enzymes (e.g., BRCA1 and RAD51) to damaged DNA (43, 44). Preclinical studies are warranted to determine whether loss of such epigenetic factors is synthetically lethal with PARP (or other DDR) inhibitors, or alterations involving their overexpression are indicative of therapeutic benefit from bona fide epigenetic inhibitors. An equally important question in the context of immunotherapy trials is on the interactions of such epigenetic alterations with immunotherapy responses. Indeed, a recent study has already suggested that SETD2 mutations are associated with high tumor mutation burden and sensitivity to ICTs possibly through impaired DNA mismatch repair (70). The immunotherapy opportunities beyond use of anti–PD-1/PD-L1 are also compelling. For example, the overlapping LAG3 expression and loss of diverse DDR mechanisms (DNA damage checkpoint, epigenetic regulation, and NER) may warrant further investigation of the vulnerabilities to LAG3 targeting in the matched patients. In summary, our studies nominated diverse DNA repair pathway alterations that co-occur with immune activities. Such alterations may guide selection of patients for emerging clinical trials or enable design of novel translational studies with combinations of DNA repair inhibitors and immunotherapies.

The contributions from bioinformatics tools to precision oncology efforts have become more needed as multiomic data-driven clinical studies (e.g., WINTHER trial) are emerging (71). The clinical implementation of REFLECT for generation of therapy recommendations tailored to patients’ multiomic signatures is a critical future goal. We expect implementation of REFLECT with an impact on clinical decision-making will be possible with the increasing availability of DNA and RNA sequencing data in the clinical setting. DNA sequencing in various forms (e.g., whole–exome and targeted gene panel sequencing) is already in clinical use. Data sets from DNA sequencing across large patient cohorts are immediately amenable to REFLECT analysis. Importantly, REFLECT does not necessitate whole-genome data and could even be implemented with large-scale panel sequencing data. Whole transcriptome data from RNA-seq are also becoming available in accredited clinical laboratories, and their use as predictive biomarkers of response is in progress (e.g., WINTHER trial; ref. 71). Although currently not approved for clinical applications, we anticipate that proteomics will also become a critical component of precision oncology pipelines in the future. This is particularly critical as most targeted therapies modulate phosphoproteomic states, which are likely better predictors of drug responses than DNA or mRNA alterations. For example, increased EGFR/HER2-activating phosphorylation levels are a predictive marker for the EGFR/HER2 inhibitors. Nearly half of ERBB2/EGFR-mutated tumors, however, do not lead to high levels of EGFR/HER2-activating phosphorylation as confirmed by our calculations. By assessing multiple data modalities at once, we expect tools such as REFLECT will improve the fraction of benefiting patients while providing more reliable target matching in precision oncology. Indeed, the validations of REFLECT showed that targeting recurrent DNA, mRNA, and protein coalterations all have beneficial outcomes, suggesting all data modalities add value. There is, however, a great need for improvement in the accuracy and completeness of the precision oncology knowledge bases to leverage multiomic data for therapy selection. We anticipate that our knowledge of how to target mRNA and protein expression alterations will improve significantly in the near future. There is also a need for a better assessment of drug combination toxicities through the use of emerging informatics resources. Therefore, it is imperative that strategies such as REFLECT remain agile and adapt to evolving multiomic data, precision oncology resources, and the clinical landscape to improve treatment selection for patients with cancer.

With expanding data from precision therapy trials and genomics studies, it will become increasingly feasible to match combination therapies to complex molecular signatures. As lack of durable responses to genomically matched monotherapies significantly limits the population of benefiting patients with cancer, REFLECT is potentially a critical addition to the current precision medicine paradigm. Our analyses suggest >40% of patients may potentially be eligible for precision combination therapies when multiomic data are utilized to identify coactionable alterations. We expect this eligibility rate will be reduced due to other critical considerations such as drug toxicity. Nevertheless, this already corresponds to a substantial improvement over the conventional paradigm that relies on only DNA alterations, as recent estimates suggest 15% of all patients with cancer are eligible for precision therapies (72). Therefore, we expect REFLECT will be a commonly used tool and resource to select precision combination therapies in both preclinical and clinical studies. We conclude that a revolution in precision medicine to significantly improve benefiting patient populations and response durations is possible through clinical bioinformatic analyses of multiple data types with predictive tools such as REFLECT.

We used phosphoproteomics data sets of TCGA patients [The Cancer Proteome Atlas (TCPA)] and cancer cell lines (MCLP) measured by reverse-phase protein array (RPPA), which were available from the TCPA Portal (https://tcpaportal.org/). Level 4 data (batch normalized across disease types) were downloaded for both patients and cell lines. To remove missing values, we filtered both data to maximize the number of samples and the number of proteins with no missing values. As a result, the patient data set had 7,544 samples by 179 (phospho-)proteins from 32 cancer types (no RPPA data were available for the LAML cancer type), and the cell line data set had 397 samples by 179 (phospho-)proteins from 18 lineages. Because of their overlapping coverages of diverse disease types, we assumed that the patient data set and the cell line data set follow the same distribution of protein expressions across samples. To remove the batch effects and merge the two data sets, we standardized each data set using Z-score transformation so that each protein had a zero mean and a unit variance. Then, the two standardized data sets were combined to form an integrated phosphoproteomics data set.

For mRNA expression, we used RNA-seq data sets of TCGA patients (N = 8,887) and CCLE cell lines (N = 1,156). The batch effects within each data set from diverse cancer types had already been corrected for and validated using the EB⁺⁺ (28). Using the same integration strategy for the phosphoproteomics data, we Z-score standardized the log₂-transformed mRNA expressions in each RNA-seq data set and then combined the two data sets together. To focus on essential genes relevant to cancer, in the current implementation of REFLECT features, we used genes from the COSMIC CGC (31) v86. The oncogenic and tumor suppressor roles of genes are assigned based on the COSMIC CGC annotations. For actionable target selection, we used |Z-score| > 2 as a threshold to detect protein and mRNA alterations and match to therapy.

For genomic alterations (mutations and CNAs), we used patient data sets from TCGA (N = 8,764) and cell line data sets from the COSMIC-CLP (N = 996). The cancer-critical genes from the COSMIC CGC v73 were used for REFLECT calculations. For the mutation data, we filtered out hypermutated samples. Then, mutations were discretized to binary values (mutant = 1, wild-type = 0). The CNAs were discretized to {−1, 0, 1}, representing deletion, wild-type, and amplification, respectively. To account for differences in exome design, depth coverage, and other batch effects, we used the multisite, multidisease, consensus pan-cancer MAF file generated by the TCGA Multicenter Mutation Calling in Multiple Cancers (MC3) project (29). For copy-number analysis, we used the significantly altered copy-number events based on GISTIC2.0 analysis of pan-cancer data (Affymetrix SNP 6.0) that had passed through the standardized TCGA data processing pipeline for all cancer types.

We selected actionable genomic alterations as the stratification markers from three sources. First, we selected 15 genomic vulnerability markers from NCI-MATCH (30). Then, from potentially actionable alterations of the MD Anderson Precision Oncology Knowledge Base (20), we extracted 182 genomic alteration markers that were not overlapped with the NCI-MATCH markers. As a result, these markers also covered the actionable genes listed in OncoKB. For each marker, the alterations included mutations, gene amplification/deletions, and specific amino acid substitutions (e.g., BRAF^V600 mutation). Depending on specific genes, activating mutations were combined with amplification, and inactivating mutations were combined with gene deletions. Finally, we constructed a cohort of patients carrying potential immunotherapy response markers and three other cohorts based on the existence of key oncogenic alterations: MYC amplifications, TP53 hotspot mutations, and KRAS mutations.

We termed these markers as “master” biomarkers (N = 201). Annotations based on Gene Ontology (GO) biological process terms were used to categorize the master biomarkers. The top three GO terms were used; terms were ranked by evidence types and counts as provided by the GO database. The code is provided at https://github.com/cannin/extract_go_slim_annotation.

Sparse hierarchical clustering identifies clusters driven by the most discriminant features (32). The input to the algorithm is X, an n × p data matrix with n samples (e.g., patients) by p features (e.g., genes). Let x_i be a sample that is a p-dimensional vector indexed by i, where 1 ≤ i ≤ n, and d(x_i, x_i′) be a measure of the distance between samples x_i and x_i′. We assume that the distance measure is additive in features; that is, |$d( {{{{\bf x}}_i},{{{\bf x}}_{i'}}} )\ = \ \mathop \sum \limits_j^p {d_{i,i',j}}$| where |${d_{i,i',j}}$| is the distance measure between samples x_i and x_i′ along feature j.

The distance is defined as either the squared Euclidean distance |${d_{i,i',j}}\ = \ {( {{X_{ij}} - {X_{i'j}}} )^2}$| or the absolute difference |${d_{i,i',j}}\ = \ | {{X_{ij}} - {X_{i'j}}} |$|⁠. The squa­red Euclidean distance, due to the square operation, assigns a relatively higher weight to larger values, whereas the absolute difference treats all values uniformly. In REFLECT, to smoothen small differences and strengthen big differences, the squared Euclidean distance was used as the distance measure for continuous data types (i.e., protein and mRNA expressions). For discrete data types (i.e., mutations and CNAs), where only a small finite set of values exists, that is, {−1, 0, 1}, we used the absolute difference as their distance measures.

The sparse hierarchical clustering can be derived from the standard hierarchical clustering formulation. Note that the standard hierarchical clustering solves a constrained maximization problem

which leads to the optimal solution |$U_{i,i'}^* \propto \mathop \sum \limits_j {d_{i,i',j}}$|⁠. The sparse hierarchical clustering can be formulated as the solution to a weighted constrained optimization problem

where w_j is a weight for feature j, and s is a tuning parameter that controls sparsity of w. Let U^** optimize the above criterion, from which we have |$U_{i,i'}^{**} \propto \mathop \sum \limits_j {w_j}{d_{i,i',j}}$| a weighted distance with feature j having a weight w_j. Due to the constraints, a small value of s makes w sparse, so U^** depends only on a small subset of features. The optimization is iterated until convergence (w is fixed whereas U is optimized, and next U is fixed whereas w is optimized). After it is converged, hierarchical clustering is performed on the optimal solution U^**. The output is hierarchical clusters of samples defined by a small set of discriminant features corresponding to nonzero elements of the converged w. In REFLECT, we used the sparcl R package (implemented by Witten and Tibshirani; ref. 32) to execute the sparse hierarchical clustering.

To evaluate the association of cancer types with the clustered hierarchy of samples, represented as a dendrogram, we calculated the concordance of clusters with cancer types as

where i denotes indices of samples along the clustered dendrogram, I is the indicator function, L_i is the label of the cancer type to which the sample i belongs, n is the number of samples, and C is the number of cancer types. The denominator n - C normalize the range of the concordance to be within 0 and 1. A large concordance implies a strong association between clustering and cancer types. When successive clustered samples always belong to different cancer types, the concordance becomes 0; when all samples of the same cancer types are tightly connected, the concordance is 1.

To evaluate the association of proliferation with the clustered hierarchy of samples, represented as a dendrogram, we calculated the Spearman correlation of sample rank along the dendrogram with proliferation scores. To calculate proliferation scores, we obtained marker genes of proliferation from Whitfield and colleagues (33), which includes MYBL2, BUB1, PLK1, CCNE1, CCND1, CCNB1, AURKA, E2F1, FOXM1, MKI67, MCM2, MCM3, MCM4, MCM5, and MCM6. Using log₂-transformed RNA-seq data of TCGA samples, the proliferation score of a sample was calculated as the mean of expressions of its proliferation genes.

The cell lines used in experiments included MFE-296, HEC50B, DU4475, BT-474, HCC1954, MDA-MB-468, SK-OV-3, KM12, HCC2998, SW480, HCT15. They were obtained through the MD Anderson Characterized Core Cell Line Facility and the MDACC Colorectal Cancer Moon Shot program. Cell lines were authenticated by fingerprinting using short tandem repeat testing. The absence of Mycoplasma contamination was also verified. Cells were maintained in RPMI 1640 supplemented with 10% heat-inactivated FBS and 1% penicillin with streptomycin. The targeted compounds [PARP inhibitor olaparib; CHK2 inhibitor BML-277; AKT inhibitor MK-2206; MET inhibitor JNJ-38877605; HER2 inhibitor neratinib; AURKA inhibitor alisertib; CDK4/6 inhibitor ribociclib; MAP3K4 inhibitor doramapimod (directly inhibiting downstream signaling partner, p38/MAPK14), PIK3CA inhibitor buparlisib] were obtained from Selleckchem.

The code for the REFLECT and data analysis is available at https://github.com/korkutlab/reflect and https://codeocean.com/capsule/8642040/tree/v1. The REFLECT portal (https://bioinformatics.mdanderson.org/reflect/) contains all the analysis results.

P.G. Pilie reports personal fees from AstraZeneca and Dendreon outside the submitted work. T.A. Yap reports grants and personal fees from AstraZeneca, Artios, Bayer, Bristol Myers Squibb, Clovis, EMD Serono, Forbius, F-Star, ImmuneSensor, Merck, Seattle Genetics, Pfizer, and Repare, grants from Ionis, Beigene, BioNTech, Eli Lilly, Glaxo­SmithKline, Genentech, Haihe, Ipsen, Jounce, Karyopharm, KSQ, Kyowa, Ribon Therapeutics, Regeneron, Rubius, Sanofi, Scholar Rock, Novartis, Tesaro, and Vivace, and personal fees from Almac, Aduro, Athena, Atrin, Axiom, Calithera, Cybrexa, GLG, Guidepoint, Ignyta, I-Mab, Janssen, Roche, Schrodinger, Varian, Zai Labs, and ZielBio outside the submitted work. S. Kopetz reports other support from Sanofi, Biocartis, Guardant Health, Array BioPharma, Genentech/Roche, EMD Serono, MedImmune, Novartis, Amgen, Lilly, and Daiichi Sankyo during the conduct of the study, as well as other support from Genentech, EMD Serono, Merck, Holy Stone, Novartis, Lilly, Boehringer Ingelheim, Boston Biomedical, AstraZeneca/MedImmune, Bayer Health, Pierre Fabre, EMD Serono, Redx Pharma, Ipsen, Daiichi Sankyo, Natera , HalioDx, Lutris, Jacobio, Pfizer, Repare Therapeutics, Inivata, GlaxoSmithKline, Jazz Pharmaceuticals, Iylon, Xilis, AbbVie, Amal Therapeutics, Gilead Sciences, Mirati Therapeutics, Flame Biosciences, Servier, Carina Biotechnology, Bicara Therapeutics, Endeavor BioMedicines, Numab Pharma, Johnson & Johnson/Janssen, Genomic Health, Frontier Medicines, Replimune, and Taiho Pharmaceutical outside the submitted work. A. Korkut reports grants from the NIH/NCI, the Cancer Prevention & Research Institute of Texas, the MDACC Colorectal Cancer Moon Shot program, the OCRA Collaborative Research Award, and the Innovation in Cancer Informatics (ICI) Fund during the conduct of the study, as well as personal fees from Inno Advisors, the Collaborative Fund, and Baylor College of Medicine Multiomics Core, and other support from the ICI Fund outside the submitted work. No disclosures were reported by the other authors.

X. Li: Conceptualization, data curation, software, formal analysis, validation, investigation, visualization, methodology, writing–original draft. E.K. Dowling: Data curation, software, formal analysis. G. Yan: Validation, methodology, writing–original draft. Z. Dereli: Conceptualization, methodology, writing–original draft, writing–review and editing. B. Bozorgui: Conceptualization, validation, methodology, writing–original draft, writing–review and editing. P. Imanirad: Data curation, formal analysis, validation. J.H. Elnaggar: Data curation, software, formal analysis, methodology, writing–review and editing. A. Luna: Resources, software, formal analysis, methodology, writing–review and editing. D.G. Menter: Conceptualization, resources, formal analysis, supervision, funding acquisition, writing–review and editing. P.G. Pilie: Conceptualization, resources, formal analysis, supervision, methodology, writing–review and editing. T.A. Yap: Conceptualization, resources, data curation, software, formal analysis, supervision, funding acquisition, validation, investigation, visualization, methodology, writing–original draft, project administration, writing–review and editing. S. Kopetz: Conceptualization, resources, formal analysis, supervision, funding acquisition, writing–original draft, writing–review and editing. C. Sander: Conceptualization, resources, formal analysis, supervision, methodology, writing–review and editing. A. Korkut: Conceptualization, resources, data curation, software, formal analysis, supervision, funding acquisition, validation, investigation, visualization, methodology, writing–original draft, project administration.

We thank Jordi Rodon Ahnert, Funda Meric-Bernstam, Scott Woodman, Gordon Mills, Rehan Akbani, Ziyi Li, Caishang Zheng, and Murat Cokol for invaluable comments. We also thank the MDACC DQS-IT team led by Chris Wakefield for deploying the REFLECT portal. This work is supported by grants from MDACC Support Grant P30 CA016672 (the Bioinformatics Shared Resource), OCRF Collaborative Research Award, ICI Fund, CPRIT High-Impact/High-Risk Award (RP170640), and funding from the MDACC Colorectal Cancer Moon Shot program, U01 CA253472-01A1, 5UL1TR003167-02, 7R01CA206025-06, and W81WH-18-1-0678, as well as NRNB-National Resource for Network Biology, NIGMS, and GM103504 (C. Sander). This manuscript was edited by Life Science Editors.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.