Integrative Pan-Cancer Genomic and Transcriptomic Analyses of Refractory Metastatic Cancer

Primary untreated tumors have been extensively studied through genomic and transcriptomic profiling, yielding valuable insights into the heterogeneity of tumors and demonstrating the utility of molecular profiling for precision oncology (1). Recently, several studies further investigated the genomic landscape of pan-cancer (2–5) or tumor type–specific (6–8) metastatic cohorts and have shown that genomic differences exist (5). Indeed, mutagenesis may differ in advanced tumors compared with primary tumors due to the impact of some widely used antineoplastic treatments, such as mutagenic platinum-based drugs (9, 10). Moreover, the therapeutic pressure imposes additional constraints that contribute to tumor evolution and acquisition of resistance (11, 12). Consequently, patients who developed resistance to multiple lines of therapies have limited treatment options and dismal prognosis. However, these patients might benefit from phase I clinical trials aiming to find new therapeutic strategies that do not cause severe side effects. An accurate estimation of the expected survival time is vital to determine patients’ eligibility for these clinical trials.

Currently, survival predictions are based on objective risk clinical markers such as lactate dehydrogenase (LDH) levels, serum albumin, neutrophil-to-lymphocyte ratio, or the number of metastatic sites (13, 14). However, these predictions do not consider tumor molecular markers or characteristics of the tumor microenvironment (15, 16), which can also be associated with the patient's survival (17–25).

This work introduces META-PRISM, a pan-cancer cohort of 1,031 tumors that progressed under at least one line of treatment or have no approved therapy options. Somatic variations, germline mutations, and tumor microenvironments were analyzed via whole-exome sequencing (WES) and RNA sequencing (RNA-seq). Identified genetic markers were compared with tumor-type matched untreated primary tumors from The Cancer Genome Atlas (TCGA; ref. 26), focusing on functional pathogenic variants associated with resistance, and validated using metastatic tumors from the MET500 cohort (3). We further investigated the utility of genomic and transcriptomic markers to improve the prediction of survival time based on objective clinical variables in the META-PRISM cohort.

The META-PRISM cohort is composed of 1,031 adult patients who were biopsied at entry into precision medicine trials (MOSCATO, ref. 27; MATCH-R, ref. 28) and subjected to WES (569 biopsies–blood pairs) and transcriptome sequencing (947 biopsies; Fig. 1A and B). The clinical history of patients, drugs received prior to biopsy (and progressed upon), tumor characteristics, and other clinical parameters derived from blood testing and physiologic assessments were tabulated (Supplementary Table S1). Details about the library preparation and sequencing quality of the analyzed samples are also available (Supplementary Table S2). The median age of the patients was 59.8 years, ranging from 43.1 to 68.9 years per tumor type (Fig. 1C). The cohort included 39 cancer types (Supplementary Table S3), with 5 of them represented by more than 60 tumors: 192 lung adenocarcinomas (LUAD), 98 breast carcinomas (BRCA), 95 prostate adenocarcinomas (PRAD), 74 bladder urothelial carcinomas (BLCA), and 61 pancreatic adenocarcinomas (PAAD; Fig. 1A and B). Moreover, 161 tumors were of rare subtypes or unknown primary (Fig. 1B) and included 42 tumors with neuroendocrine differentiation (Supplementary Fig. S1). We focused our genetic analysis on the tumor types (excluding tumors of unknown origin) that were represented by at least 10 sequenced tumors, resulting in 10 types for WES (META-PRISM WES, 83.8% of all DNA samples), 10 types for WES and RNA-seq (META-PRISM WES and RNA-seq, 82.1% of all DNA and RNA samples), and 20 for RNA-seq (META-PRISM RNA-seq, 84.2% of all RNA samples; Supplementary Table S2; Fig. 1B).

At the time of the biopsy, patients had a median number of two metastases distributed at 35 different sites (Fig. 1C; Supplementary Fig. S2). The most frequently biopsied metastatic sites were liver (37%), lung (22%), and lymph nodes (15%; Fig. 1A; Supplementary Fig. S2). Before the biopsy, patients had received a median number of four treatments, with BRCA and ovarian carcinomas (OV) being the most pretreated tumors (medians of six and nine treatments, respectively; Supplementary Table S1). Most treatments were tumor-type specific, such as EGFR inhibitors in LUAD or hormone therapies in BRCA or PRAD. Other therapies, particularly platinum-based drugs, were applied to multiple tumor types (Fig. 1D; Supplementary Table S1).

The median survival time from the biopsy date was 7.8 months (n = 909 deceased patients out of 1,029; Fig. 1C) and was highly heterogeneous across tumor types. PAAD, BLCA, carcinomas of unknown primary and lung squamous cell carcinomas (LUSC) had very poor prognosis (<6 months median survival), whereas cervical squamous cell carcinomas and endocervical adenocarcinomas (CESC), OV, sarcomas (SARC), and kidney renal clear-cell carcinomas had longer survival (>12 months; Fig. 1C and E).

In order to provide a comparative assessment of the types and extent of genetic variations in META-PRISM, we took advantage of the TCGA cohort representing primary untreated tumors (Supplementary Fig. S3A and S3B). Furthermore, we reanalyzed the previously published MET500 pan-cancer cohort of metastatic tumors to validate discoveries made on META-PRISM (Supplementary Fig. S3C and S3D).

WES analysis revealed a significant increase of somatic mutations in META-PRISM versus TCGA for 6 of 10 studied tumor types. Tumors from the MET500 cohort demonstrated a tumor mutational burden (TMB) similar to that of META-PRISM tumors. We observed the most significant increase of TMB in low-burden tumor types BRCA (2.1× fold in META-PRISM, 2.7× fold in MET500), PRAD (2.1× fold in META-PRISM, 2.2× in MET500), and PAAD (1.7× fold in META-PRISM), whereas no such increase was observed in high-burden tumor types, namely, BLCA, LUAD, and LUSC (Fig. 2A).

Deconvolution of mutational signatures revealed similar signature compositions between metastatic META-PRISM, MET500, and primary TCGA cancers in all studied tumor types (Supplementary Fig. S4 and Methods). However, a notable and consistent difference in META-PRISM versus TCGA across tumor types was the presence of signatures associated with platinum treatments (SBS31 and SBS35), reflecting that the majority of META-PRISM tumors (691 of 1,011 with known drug history) were pretreated with platinum compound therapies. Signatures SBS31 and SBS35 were detected in more than 50% of tumors in six META-PRISM tumor types [LUAD, BLCA, LUSC, HNSC, cholangiocarcinomas (CHOL), and adrenocortical carcinomas (ACC)] and contributed significantly more mutations compared with TCGA in four of them (Fig. 2B). Tumor types varied in frequency and types of received platinum drugs. For example, BRCA and PRAD rarely received platinum treatments and concordantly demonstrated a very low platinum signature. We then investigated the association of SBS31 and SBS35 with three main platinum drugs in our cohort: cisplatin, carboplatin, and oxaliplatin. Among these drugs, cisplatin had the strongest association with SBS31 and had an association with SBS35 that was comparable with carboplatin, as revealed by logistic regression (Fig. 2B). More specifically, among tumors harboring at least 50 somatic mutations (see Methods) and treated with cisplatin and no other platinum compound (n = 95), 40% harbored SBS31 and 17% SBS35, whereas in patients treated only with carboplatin (n = 80), 17% had a detectable activity of SBS35. In contrast, these two signatures were rarely detectable in tumor types (COAD and PAAD) predominantly treated with oxaliplatin (Fig. 2B).

Whole-genome dupli­cation (WGD) is a frequent event in cancer involving doubling the chromosome complement. We detected WGDs in 53.7% of META-PRISM WES tumors, which was comparable with MET500 (50.2%) but significantly higher than in TCGA (39.9%; P = 0.01, U test; see Methods). The fraction of WGDs varied between tumor types ranging from 25.0% in CHOL to 82.6% in LUSC in META-PRISM. The most striking increase of WGD events in metastatic cancers compared with TCGA was observed in PRAD (META-PRISM 48.5%, 12× fold increase, P < 0.001; MET500 38.6%, 9.4× fold increase, P < 0.001; Mann–Whitney U tests) and PAAD (META-PRISM 47.4%, 3.1× fold increase, P < 0.001; MET500 37.5%, 2.5× fold increase, not significant due to small size; Fig. 2C).

We next investigated the landscape of somatic copy-number alterations (SCNA) in META-PRISM WES tumors without WGD. Using FACETS (29), we categorized SCNAs in META-PRISM, MET500, and TCGA as copy gains and losses. Copy gains were subdivided into either low-level (LLG), middle-level (MLG), or high-level (HLG) gains. Copy losses included loss-of-heterozygosity (LOH), copy-neutral LOH (cn-LOH), and focal homozygous deletions (HD). LLGs, MLGs, LOH, and cn-LOH often spanned large regions or covered full chromosome arms, whereas HDs and HLGs were almost always focal. The fraction of the genome covered by low- to medium-level gains or LOH was most drastically increased in metastatic PRAD tumors. Additionally, LOH events were significantly more frequent in META-PRISM BRCA and PAAD tumors (Fig. 2D). The SCNA profiles of META-PRISM tumors were overall similar to that of TCGA tumors for both the tumor type–adjusted dataset (Fig. 2E) and specific tumor types (Supplementary Fig. S5A–S5E). The frequency of the majority of large SCNAs did not differ significantly.

In META-PRISM WES tumors without WGD, an average of 10.7% and 19.7% of the genome were covered by copy gains and copy losses, respectively. No significant increase in this type of instability was observed in most studied tumor types except for PRAD, which demonstrated a dramatic increase in metastatic tumors compared with primary tumors, and for PAAD, in which the increase was limited to copy losses (Fig. 2D; Supplementary Fig. S5D and S5E). Three chromosome arm gains (5p, 7p, and 8q) and seven losses (6q, 8p, 9p, 13q, 17p, 18p, and 18q) were observed in more than 20% of the non-WGD META-PRISM cohort. However, their frequency was not significantly different from TCGA. The majority of these chromosome regions enriched with gains and losses events were shared between several tumor types, whereas some chromosome regions were tumor-type specific: for instance, +16p, −16q, −22q in BRCA, +20p, +20q, −9q, and −11p in BLCA (Supplementary Fig. S6A and S6B). META-PRISM tumors with WGD did not demonstrate significant differences in arm-level copy-number losses, focal amplifications, focal deletions, or number of driver mutations as compared with tumors without WGD in the five main tumor types of the WES subcohort (Supplementary Fig. S7A–S7C).

High-level amplifications and HD were rare in the META-PRISM WES tumors, spanning an average of 0.15% and 0.05% of the genome, respectively. However, highly amplified and homozygously deleted genes were detected on average 63 and 13 times per tumor. The most frequent of these events included amplification of CCND1 (8.2%), AR (4.4%), 19q13 genes (2.7%), EGFR (2.7%), MYC (2.7%), and KRAS (2.5%; Supplementary Fig. S7D) and losses of CDKN2A (13%), FAM106A/LGALS9C (5.0%), Killer Ig–like receptors genes (3.4%), RHD/RSRP1 (2.9%), and PTEN (2.7%; Supplementary Fig. S7E).

The discovery of significantly mutated genes with Mutpanning (30) confirmed previously reported cancer drivers (Supplementary Table S4; Supplementary Fig. S8). Interestingly, this analysis highlighted genes that were reported as drivers in advanced tumors but not in primary tumors, namely, EP300 in META-PRISM LUAD, HERC2 in META-PRISM BRCA, ESR1 in META-PRISM and MET500 in BRCA, and AR in META-PRISM and MET500 in PRAD. Mutations in these drivers may have been selected by the therapies or have arisen in the late stage of tumor evolution. We next selected a list of 360 cancer genes by intersecting COSMIC census “Tier 1” and “Tier 2” (v92) and OncoKB-annotated genes (ref. 31; OncoKB 07/2021; Supplementary Table S5) and created a catalog of driver mutations and SCNAs (point mutations, small indels, amplifications, and deletions) using OncoKB annotations on these genes only (Supplementary Tables S6 and S7). These types of driver events were observed in 96% of the META-PRISM WES tumors. The most frequently altered driver genes in META-PRISM were TP53 (55% of samples), KRAS (25%), CDKN2A (18%), and EGFR (14%; Fig. 3A). On the whole cohort level, 10 oncogenes and 3 tumor suppressor genes were significantly enriched compared with TCGA; 40% and 100% of those genes, respectively, were also enriched in MET500. Some driver genes were significantly enriched in specific tumor types: for example, EGFR and CTNNB1 in LUAD; TP53, AR, PTEN, and RB1 in PRAD; ESR1 and CCND1 in BRCA; TP53 and KRAS in PAAD; and FGFR3 in BLCA (Fig. 3A, Fisher–Boschloo tests, Benjamini–Hochberg correction; see Methods). The number of WES-derived driver events (mutations and SCNAs) was significantly higher in metastatic tumors as compared with primary tumors at the cohort level (means of 3.9 and 3.6 in META-PRISM and MET500 vs. 3.2 in TCGA, P < 0.0001) and for 5 of 10 tumor types included in the META-PRISM WES cohort (Fig. 3B).

Major tumor suppressor genes frequently underwent biallelic inactivation. Such inactivations were observed in 92% of TP53-hit tumors, 95% for CDKN2A, 81% for PTEN, 88% for SMAD4, 41% for ARID1A, 78% for APC, and 92% for RB1. However, the predominant mechanisms of biallelic inactivation differed from one gene to another: in TP53, it was mutation followed by LOH; in CDKN2A and PTEN, it was HD; and other genes demonstrated a combination of mechanisms (Supplementary Fig. S9A–S9F). Few oncogenes also underwent multihit events, most notably EGFR (57%), AR (12%), PIK3CA (7.7%), and KRAS (4.2%). Multihit events in EGFR in META-PRISM were significantly more frequent than in TCGA, likely reflecting the effect of EGFR inhibitors in LUAD (Supplementary Fig. S10A–S10F).

WES of germline DNA from 569 META-PRISM patients was used to identify pathogenic cancer-predisposing variants. We focused on the germline variants annotated as pathogenic or likely pathogenic in the ClinVar database or protein-disrupting and residing in genes strongly associated with cancer predisposition as described in Huang and colleagues (ref. 32; see Methods). We identified 75 patients in META-PRISM (13.1%) harboring at least one such variant. The fraction of patients with cancer-predisposing variants was similar to that in MET500 patients and was 1.75 times higher than in TCGA (8 cancer types, P = 0.0012, Fisher exact test) or in ExAC non-Finnish Europeans (OR = 1.75, P = 0.0002; Fig. 3C). 73% of variants were in DNA repair pathway genes (Fig. 3D). The most frequent genes with cancer-predisposing variants in META-PRISM patients were MUTYH (1.9%), BRCA2 (1.8%), NF1 (0.9%), ERCC2 (0.9%), ERCC3 (0.9%), RAD51D (0.7%), and FANCG (0.7%; Fig. 3A). MUTYH, NF1, and RAD51D were also significantly enriched compared with TCGA. We detected an increase of germline cancer-risk variants in most cancer types in META-PRISM. However, it reached significance only for PRAD (P = 0.03) and LUSC (P = 0.004) cancer types. Mutations in the homologous recombination pathway were the most frequent in BRCA, PRAD, LUAD, and BLCA; base excision repair pathway in PAAD and LUSC; and mismatch repair pathway in COAD (Fig. 3E). Thirty-seven percent of genes with germline cancer-risk variants harbored a somatic second-hit event, including somatic mutations in 9% and LOH resulting in the retention of the pathogenic variant in 27%.

The tumor microenvironment (TME) plays a significant role in clinical outcomes and response to therapy (17, 18, 33, 34). Tumors from META-PRISM, MET500, and TCGA were classified into four tumor microenvironment subtypes—immune-enriched, fibrotic (IE/F); immune-enriched, nonfibrotic (IE); fibrotic (F); and immune-depleted (D)—using functional gene expression signatures scores as described in Bagaev and colleagues (ref. 16; see Methods). No consistent difference in the distribution of TMEs across tumor types was observed when comparing META-PRISM to TCGA. However, the TME subtypes in some individual tumor types showed striking variation between the cohorts. For instance, immunosuppressive subtypes D and F were significantly increased in PRAD and BLCA, respectively, in META-PRISM compared with TCGA (63% vs. 41%, P < 0.001; 32% vs. 17%, P = 0.008; P values are not adjusted due to the lack of independence; Fig. 4A). Enrichment of the D subtype in PRAD was also significant in MET500 versus TCGA (78% vs. 41%, P < 0.001).

We next investigated gene fusions in the 795 META-PRISM RNA-seq tumors that were successfully processed by the fusion pipeline. The calling of fusions was restricted to known oncogenic fusions and fusions involving a cancer driver, as only these two types of fusions are most likely relevant to cancer and could reliably be detected in all three cohorts (see Methods).

A total of 432 known oncogenic gene fusions (Chimer KB v4.0, ref. 35; Chitars v5.0, ref. 36; COSMIC v95; TIC v3.3, ref. 37) were identified in META-PRISM RNA-seq tumors (34%; Supplementary Table S8). As previously described, well-known oncogenic gene fusions were often tumor-type specific—TMPRSS2–ERG 28% in PRAD, EML4–ALK 6% in LUAD, DNAJB1–PRKACA 25% in LIHC, MYB–NFIB 64% in head and neck adenoid cystic carcinomas (HNAC, rare subtype)—and some were significantly enriched in META-PRISM versus TCGA (Fig. 4B). In addition to known oncogenic fusions, we identified 329 fusions involving cancer driver genes and promiscuous partners (29% of META-PRISM RNA-seq tumors; Supplementary Table S8).

Among the most recurrently fused driver genes, we identified tumor suppressors PTEN (1.9%), TP53 (0.9%), and RB1 (0.8%). In these genes, 72% of fusion breakpoints were recurrent (Fig. 4C). PTEN fusions were observed in several tumor types in META-PRISM, but they were most prevalent in PRAD (Fig. 4B). The oncogene most frequently involved in fusions across different tumor types was FGFR2 (1.5%). Interestingly, ESR1 and AR, both known to be important factors of hormone therapy resistance in BRCA and PRAD, were involved in fusions in 2.2% and 4.7% of the respective tumor types in META-PRISM (Fig. 4B and C). Sixteen of 18 tested fusions (89%) were validated through RT-PCR and Sanger sequencing (Fig. 4D; Supplementary Table S9).

Biomarkers of resistance and response to anticancer therapies were annotated in META-PRISM, MET500, and TCGA by mapping somatic events to the OncoKB and CIViC databases (refs. 31, 38; Supplementary Tables S6–S8). ESCAT guidelines (39) were used to classify biomarkers as Tier 1, Tier 2, or Tier 3. Tier 1 is standard-of-care (SOC) markers that are used in routine clinical practice for treatment indications; Tier 2 is investigational markers supported by clinical trials but for which additional data are needed; and Tier 3 is hypothetical targets supported by scarce or off-label data.

We observed SOC resistance and sensitivity biomarkers in 9.6% and 47.5% of META-PRISM WES and RNA-seq tumors, respectively (Fig. 5A and B), whereas biomarkers of any level were found in 74.9% and 88.4% of these tumors, respectively. Tier 1 resistance biomarkers were detected in only three tumor types: LUAD (EGFR 17%), LUSC (EGFR 5%), and COAD (KRAS 65%, NRAS 6%), whereas investigational Tier 2 and hypothetical Tier 3 biomarkers were rather frequent in the majority of tumor types (Fig. 5A). META-PRISM WES and RNA-seq tumors harbored significantly more resistance biomarkers of all three tiers compared with TCGA: Tier 1, common odds ratio (cOR) 7.5 [confidence interval (CI) 95%, 3.7–15.2, Cochran–Mantel–Haenszel test stratified by tumor type]; Tier 2, cOR 1.7 (CI 95%, 1.3–2.2); and Tier 3, cOR 2.2 (CI 95%, 1.7–2.8). The increase was also replicated in MET500 for all three tiers (Tier1 cOR 4.7; CI 95%, 1.3–17.3; Tier2 cOR 2.2, CI 95%, 1.6–2.9; Tier3 cOR 3.2, CI 95%, 2.4–4.3).

Some common sensitivity biomarkers were frequent in META-PRISM and often enriched versus TCGA, most notably EGFR L858R/Exon 19 del (first- and second-generation EGFR inhibitors), EGFR T790M (third-generation EGFR inhibitors), and ALK oncogenic mutations or fusions (ALK inhibitors) in LUAD; PTEN fusions or loss-of-function mutations (mTOR inhibitors) in PRAD; and FGFR3 p.R248C, p.S249C, p.G372C, and p.Y375C mutations (FGFR inhibitors) in BLCA. On top of that, we detected 13% of hypermutated tumors (>10 mut/Mb) and 3% of microsatellite-unstable tumors in META-PRISM, which are indications for treating patients with immunotherapies (Fig. 5B).

We next used the alterations described above to uncover the associations between ESCAT resistance markers and the treatments received before the biopsy in META-PRISM WES and RNA-seq tumors. Additionally, we retrieved “emerging” resistance markers that are supported by the literature but are not yet included in the versions of OncoKB and CIViC databases used for this analysis (Supplementary Table S10).

We could only identify resistance mechanisms for a minority of treatments (Fig. 6A–G). SOC resistance events were detected only for first- and second-generation EGFR inhibitors in LUAD for 40% of treated patients and EGFR antibodies in COAD for one in five treated patients (Fig. 6A and D). In contrast, investigational and hypothetical resistance markers substantially increase the fraction of observed treatment resistances that may be explained. In LUAD, Tier 2 resistance alterations were identified in 8 of 28 patients treated with third-generation EGFR inhibitors, in 2 of 17 patients who progressed on first- or second-generation ALK inhibitors, and in 5 of 27 patients who progressed on immunotherapies (Fig. 6A). BRCA patients treated with hormone therapy harbored Tier 3 resistance through ESR1 mutations in 33% of cases (Fig. 6B). PRAD patients treated with hormone therapy harbored Tier 2 resistance through AR alterations in 31% of cases, 25% of which were the overexpression of AR-V7 splice isoform (Fig. 6C). In COAD, Tier 2 resistance alterations could explain resistance to EGFR antibodies in two of five treated patients (Fig 6D). In PAAD, we did not identify any resistance markers for the drugs patients received.

Additionally, we were able to associate “emerging” resistance mechanisms with a considerable fraction of resistance to targeted and hormonal therapies. For example, 13%, 32%, and 67% of LUAD patients treated with first- or second-generation EGFR inhibitors, third-generation EGFR inhibitors, and BRAF inhibitors (dabrafenib), respectively, harbored “emerging” resistance mechanisms with literature support (Fig. 6A; Supplementary Table S10). Interestingly, one out of the two patients who progressed on lorlatinib (third-generation EGFR inhibitor) harbored an EML4–ALK fusion and a double ALK mutation (p.F1174L. and p.G1202R) as previously described (40). In HNSC, we associated driver mutations in EGFR, NF1, PIK3CA, and PTEN to the EGFR inhibitor cetuximab (Fig. 6E; Supplementary Table S10). In BLCA and CHOL, we were able to associate mutations in FGFR1, FGFR2, PIK3CA, KRAS, and TSC1 with resistance to FGFR inhibitors (Fig. 6F and G; Supplementary Table S10).

The observations above show that SOC resistance biomarkers were associated with only 1.6% of received treatments, whereas investigational, hypothetical, and “emerging” mechanisms could further explain 2.7%, 2.3%, and 7% of resistance, respectively. Strikingly, except for BRAF V600E mutation in COAD, no known or emerging mechanism of resistance was observed in our cohort for chemotherapies or antiangiogenic treatments, though these two types of treatments account for 54% of treatments administered.

The association of molecular markers with overall survival measured from the day of diagnosis has been the subject of extensive work (19–25), some of which has become SOC. However, the association of these markers with survival at the very late stage has been scarcely explored, partly because it is not current clinical practice to systematically profile advanced tumors and partly because the clinical value of such sequencing is still an open question.

We investigated how driver genes are univariably associated with metastatic cancer patient survival in META-PRISM WES and RNA-seq tumors. Twenty-one genes that were altered (SCNA, mutations, fusions) in more than 5% of tumors were considered for the analysis. Only TP53 and CCND1 genes showed an association with survival (adjusted HR = 1.44; CI 95%, 0.99–2.02 and HR = 1.82; CI 95%, 0.98–3.38, respectively) after adjusting for the tumor-type composition (Fig. 7A, upper row). Analyses per tumor type revealed an association of SMARCA4 and STK11 events with poor prognosis in LUAD and of ERG events with favorable prognosis in PRAD (Fig. 7A, middle row).

TME classifications, which have already been shown to be prognostic (15, 16), were associated with survival in META-PRISM. Indeed, immune-cold subtypes (F and D) had worse survival compared with immune-enriched subtypes (IE/F and IE). Strikingly, the F subtype in BLCA and D subtype in PRAD, both enriched in META-PRISM versus TCGA, were also associated with the poorest prognosis in the corresponding tumor type (P = 0.002, F vs. non-F BLCA; P = 0.02, D vs. non-D in PRAD; Fig. 7A, bottom row).

Patients with advanced metastatic cancer are characterized by severe physiologic deterioration as measured by LDH levels, serum albumin, or neutrophil-to-lymphocyte ratio. The Gustave Roussy Immune Score (GRIM score; ref. 13), which combines these physiologic markers, could predict 6-month survival in the META-PRISM patients with an average concordance index of 0.67 (Fig. 7B and Methods). We then investigated if the prediction of survival at this stage of the disease could be improved by considering genetic markers engineered from the previous analyses. The list of genetic markers retrieved from WES was composed of summary statistics, including SCNAs, WGD status, microsatellite instability (MSI) score, TMB, and the presence of oncogenic alterations aggregated into genes or pathways discriminating between ESCAT levels. The list of markers retrieved from RNA-seq included TME subtypes and gene expression signatures measuring the activity of immune pathways, activation of general pathways (41), or main transcription factors (TF; ref. 42). In order to compare the added prognostic value from different categories of markers to the objective clinical markers, we ran multiple Cox proportional hazard regressions on the META-PRISM WES and RNA-seq subcohort and each of the five main tumor types considering only the samples with both WES and RNA-seq (see Methods).

As expected, the models’ performances based on discrete GRIM score (M1 model) were significantly improved by using the continuous metrics underlying it and by adding baseline clinical variables (age, gender, tumor type, clinical subtypes; M2 model). Nevertheless, the incorporation of WES-derived markers (M3 and M4 models), RNA-seq–derived markers (M5 and M6 models), or both (M7 model) resulted in a further increase of the C-index over the continuous GRIM score model (M2) in all analyzed tumor types and for the full cohort except for LUAD, in which the discrete GRIM score classification was the most prognostic (Fig. 7B).

Strikingly in META-PRISM BRCA tumors, for which the GRIM score alone (M1, C-index = 0.469) had no prognostic value and the GRIM score components along with clinical variables including IHC subtypes was moderately prognostic (M2, C-index = 0.620), the inclusion of genetic markers engineered from WES or RNA-seq considerably improved the 6-month survival prediction (C-index = 0.784). Analysis of the model's coefficients show that multiple markers have significant prognostic value independently of the baseline clinical variables, which include the IHC subtypes, namely: (i) the total number of Tier 1 alterations [log hazard ratio (LHR) 3.75; CI 95%, 1.07–6.53], the fraction of the genome covered by (ii) focal deletions (LHR −4.05; CI 95%, −8.73, −0.34), by (iii) focal amplifications (LHR 7.28; CI 95%, 2.63–13.55), and by (iv) middle- to low-level gains (LHR −3.65, CI 95%, −8.18, −1.22), and (v) the TME subtypes (LHR −3.18 for F vs. D, −2.72 for IE vs. D, −4.80 IE/F vs. D; Supplementary Fig. S11A and S11B).

The predicted risk scores from this model were significantly higher in patients with a survival time greater than 6 months than in patients with a survival time lower than 6 months from the date of the biopsy (Fig. 7C). By discriminating between patients having a positive risk score (poor prognosis, 8/51) and a negative risk score (good prognosis, 43/51), the model was able to split BRCA patients into two groups with very different clinical outcomes (P < 0.0001 log-rank test; Fig. 7D).

In order to validate the survival model on an external cohort, we took advantage of the large TCGA BRCA cohort and considered all patients with survival, mutation, gene expression, and gene fusion data available (n = 670). Due to the unavailability in TCGA of clinical variables underlying the GRIM score, we considered one extra model (M7bis) that was identical to M7 except that the GRIM score components were excluded. This model was trained on META-PRISM BRCA and validated on TCGA BRCA. It was found to be predictive of survival in META-PRISM (C-index = 0.71, cross-validation; P < 0.0001 log-rank test) and in TCGA (C-index = 0.65 at 1 year and C-index = 0.63 at 6 years; P = 0.0068; Fig. 7E and F).

The META-PRISM project provides genetic and transcriptomic variation for a large cohort of refractory tumors from 39 tumor types, 20 of which were analyzed in depth. This cohort is characterized by a short survival time after the biopsy date and by a high proportion of multiresistant tumors or rare tumors with no approved therapy options. Consequently, the genomes of these tumors represent a much-advanced evolutionary stage containing the footprints of mutagenic treatments and therapeutic pressure. To characterize the genetic traits specific to this cohort and assess how these may inform the tumor's aggressiveness and resistance to therapies, we compared META-PRISM to >10,000 primary untreated tumors from TCGA (1, 43) and validated all results using an external cohort of 500 metastatic tumors (3).

Our analysis reveals that several types of genomic instability were strongly enriched in refractory cancers, particularly the mutation rate, the frequency of WGDs, and the fraction of the genome covered by focal copy-number alterations. These results are consistent with previous studies that reported increased genomic instability and mutational burden in metastases of different cancer types (2, 5, 44–47). Correlative analyses between the mutational profiles of tumors and the history of treatments received have shown increased activity of signatures SBS31 and SBS35 in tumors treated with cisplatin and to a lesser extent with carboplatin. The high mutation rate of these two signatures supports previous observations of strong mutational footprints caused by platinum treatment in metastatic cancers (10) and relatively low mutagenic effects of oxaliplatin (48). The driver fusions and SCNAs represented 9.4% and 18.8% of all detected variation in driver genes, respectively, and were strongly enriched in META-PRISM tumors compared with TCGA tumors. However, these two types of events were associated with treatment resistance in 0% and 5.2% of META-PRISM patients, respectively, reflecting the scarcity of the current knowledge about genetic events conferring drug resistance.

Our research also reveals that current annotations of resistance mechanisms can only account for a small proportion of all observed treatment resistance, as resistance markers were found in 74.9% of patients, but for patients who received a given treatment, SOC markers could explain only 1.6% of all resistance; investigational, hypothetical, and “emerging” mechanisms could further explain 2.7%, 2.3%, and 7% of resistances. These low percentages were observed even though patient inclusion in MOSCATO and MATCH-R trials was driven by the possibility of using innovative treatments. These data highlight the unmet need for large-scale efforts that combine molecular profiling with exhaustive clinical annotations to fill our current lack of understanding of resistance mechanisms in cancer.

Interestingly, MSI and high TMB were detected in 3% and 13% of META-PRISM tumors, respectively, suggesting that a considerable fraction of metastatic cancers refractory to conventional treatments could benefit from immune-checkpoint inhibitor therapies, such as pembrolizumab (49, 50). Additionally, a higher fraction of pathogenic cancer-risk germline variants was observed in META-PRISM than in the general population or the comparison cohort of primary tumors. The increased incidence of predisposing mutations suggests that patients with metastatic cancer could benefit from germline testing and genetic counseling.

Markers of physiologic deterioration, such as levels of albumin, neutrophils, lymphocytes, and LDH, are used in objective risk scoring systems to predict patient survival, which is essential for assessing the eligibility of these patients for phase I clinical trials (51, 52). This study shows that tumor genomic and transcriptomic features can be used to improve the accuracy of predictions based on objective risk factors at this stage of the disease. The added value of genetic markers to current prognostic scores is currently limited at the pan-cancer level, likely due to the high heterogeneity of mechanisms driving tumorigenesis in each tumor type. However, models incorporating genetic markers are significantly more accurate in predicting 6-month survival in BRCA refractory tumors, showing that WES and RNA-seq are important for accurately establishing the individual risk profile of patients with late-stage cancer.

This study validates the previous combined genomic and transcriptomic descriptions from pretreated pan-cancer metastatic diseases (53, 54). It also highlights the feasibility of precision medicine in clinical routine and the benefit of access to new drug protocols for patients without standard treatment (54, 55).

The present cohort advances translational cancer genomics by providing a unique resource combining detailed clinical data with exome and transcriptome profiling.

Patients 18 years or older with unresectable or metastatic cancer that had progressed during at least one line of prior therapy or lacked treatment options, that had at least one site accessible to biopsy, and that were considered noncurable by a multidisciplinary board were included as part of the MATCH-R (NCT02517892) and MOSCATO (01/02; NCT01566019) clinical studies. All patients provided written informed consent, and the studies were approved by the Institutional Review Board of Gustave Roussy and were executed in accordance with the Declaration of Helsinki and Good Clinical Practice guidelines. In this study, we included biopsies of patients enrolled in MATCH-R, MOSCATO, or both (302, 697, and 32 biopsies, respectively). In cases in which multiple biopsies were available (84 cases), the latest sequenced biopsy was selected for the study (with the exception that preference was given to samples with DNA and RNA data). All patients have given consent for the trial, genomic/transcriptomic analyses, and data sharing for cancer research purposes. Needle biopsies were sampled from a metastatic lesion whenever possible or the primary tumor otherwise and were frozen in liquid nitrogen. A senior pathologist assessed the tumor cellularity on the hematoxylin and eosin slide from the biopsy that was used for DNA and RNA extraction. Samples with tumor cellularity >30% were selected for genetic analysis.

Clinical variables were manually collected by clinicians or extracted from electronic health records stored on the Dr. Warehouse database (56) using semiautomatic methods. Semiautomatic data collection consisted in extracting text chunks containing specific patterns of interest that were further processed by automatic scripts and subjected to manual inspection and correction by three different oncologists. Through these techniques, we were able to collect data regarding patients’ demographics (age at biopsy, gender), tumor characteristics (primary tumor site, tumor histology, date of first cancer diagnosis, metastatic status at diagnosis, number of metastatic sites, and localization at the time of MOSCATO/MATCH-R inclusion), blood test results performed within 1 month of the tumor biopsy date, and history of treatments received before the biopsy. All anticancer drugs were classified into either conventional chemotherapy (further subdivided into platinum salts, alkylating agents, and nonalkylating agents), targeted therapies (further subdivided on the basis of the drug target, e.g., EGFR inhibitors, ALK inhibitors, BRAF inhibitors), hormone therapies, or immunotherapies (further subdivided into immune-checkpoint inhibitors and non–immune-checkpoint inhibitors). Patients were considered resistant to a systemic treatment if they had progressed under or after the treatment and before the biopsy date. All clinical variables are tabulated in Supplementary Table S1.

The GRIM score (52) was calculated based on the blood test results measured within 1 month of the biopsy date (LDH levels, albumin concentration, neutrophil counts, and lymphocyte counts). Survival time was calculated from the date of biopsy to the date of death from any cause. Out of 1,031 patients, dates of death and dates of last news were available for 1,029 patients (99.81%). At the time of the present study, 909 patients had died (88.17%). Follow-up time ranged from 2 days to 2,667 days (Q1 127–Q3 616 days).

Tumor DNA and whole-blood germline DNA were extracted using the AllPrep DNA/RNA Mini Kit (Qiagen) and DNeasy Blood and Tissue Kit (Qiagen) according to the manufacturer's instructions. Library preparation, exome capture, sequencing, and data analysis were performed by IntegraGen SA.

Genomic DNA was captured using Agilent Human All Exon V5, CR, and CR2 capture kits (for samples processed before the year 2019) or Twist Human Core Exome Enrichment System (Twist Bioscience; for samples processed after the year 2019). Sequence capture, enrichment, and elution were performed according to the manufacturer's instructions and protocols without modifications. For library preparation, at least 150 ng of each genomic DNA was fragmented by sonication and purified to yield fragments of 150 to 200 bp. Paired-end adapter oligonucleotides from the NEB kit were ligated to repaired and a-tailed fragments, and then purified and enriched by 4 to 7 PCR cycles. Then 500 ng of these purified libraries was hybridized to the SureSelect or Twist oligo probe capture libraries. After hybridization, washing, and elution, the eluted fraction was PCR-amplified with 8 to 10 cycles, purified, and quantified by qPCR to obtain enough DNA templates. Each DNA sample was then sequenced on an Illumina HiSeq 2000 as paired-end 75 bp reads (before the year 2019) or on an Illumina NovaSeq as paired-end 100 bp reads (after the year 2019). Image analysis and base calling were performed using Illumina Real-Time Analysis with default parameters. Supplementary Table S2 provides details about the kits used for each DNA sample.

Quality control of paired-end reads was accomplished using FastQC v0.11.8. Fastp v0.20 (57) was subsequently used for trimming adapters and polynucleotide tracts from reads that were longer than 25 nucleotides. Afterward, the resulting cleaned FASTQ files were aligned to the reference human genome GRCh37 using BWA-MEM v0.7.17 (58). Intermediate BAM files were further processed for deduplicating reads using MarkDuplicates from Picard v2.20.3, sorting coordinates using SAMtools v1.9 (59), and finally recalibrating their base qualities using BaseRecalibrator and ApplyBQSR. All these tools are included in the GATK bundle v4.1.8.1 (60). Alignment quality was controlled using three different algorithms: mosdepth v0.2.5 (61), flagstat from SAMtools v1.9 (59), and CollectHsMetrics from GATK v4.1.8.1.

Tumor RNA was extracted using the AllPrep DNA/RNA Mini Kit (Qiagen) according to the manufacturer's instructions. RNA libraries were prepared with the TruSeq Stranded mRNA kit (for samples processed before the year 2019) or with the NEBNext Ultra II Directional RNA Library Prep Kit for Illumina protocol (for samples processed after the year 2019) according to supplier recommendations.

The purification of PolyA-containing mRNA molecules was performed using oligo(dT) coupled to magnetic beads from 1 μg total RNA (with the Magnetic mRNA Isolation Kit from NEB). Fragmentation was performed using divalent cations under elevated temperature to obtain 300 to 400 bp fragments, followed by double-strand cDNA synthesis, Illumina adapter ligation, and cDNA library amplification by PCR for sequencing. Paired-end read sequencing was then carried out on Illumina NextSeq 500 (75 bp; before the year 2019) or on Illumina NovaSeq (100 bp; after the year 2019). Image analysis and base calling were performed using Illumina Real-Time Analysis (3.4.4) with default parameters. Supplementary Table S2 provides details about the kits used for each RNA sample.

Quality control of reads was performed with Trim galore v0.4.4. Reads were then pseudoaligned to the human transcriptome from GENCODE v27 using Kallisto v0.44.0 (62), and transcript-level estimates were aggregated to the gene level (58,288 genes) using TxImport v1.16.0 (ref. 63; see Data Tables 1 and 2).

Germline single-nucleotide variants (SNV) and indels were called using HaplotypeCaller (https://doi.org/10.1101/201178). After calling, putative germline variants were filtered using hard thresholding on QualByDepth (QD > 2), genotype quality (QUAL > 30), FisherStrand (FS < 60 for SNV, < 200 for indels), ReadPosRankSumTest (ReadPosRankSum > −8 for SNVs, > − 20 for indels), RMSMappingQuality (MQ > 40 for SNVs only), and MappingQualityRankSumTest (> −12.5 for SNVs only) as recommended in GATK best practices (see https://gatk.broadinstitute.org/hc/en-us/articles/360035890471-Hard-filtering-germline-short-variantsh). Variants that passed all filters were annotated using ANNOVAR (64).

Somatic point mutations and small indels were detected using Mutect2 (65). In order to remove artifacts and false positives, a panel of normal was created from normal samples and used at the Mutect2 calling step as specified in the GATK best practices (see https://gatk.broadinstitute.org/hc/en-us/articles/360035890631-Panel-of-Normals-PON-). Putative variants were then analyzed for read orientation artifact and sample contamination by running GATK LearnReadOrientationModel and CalculateContamination again as recommended in the best practices.

The following set of filters was applied:

• Not filtered by Mutect2 (MUTECT_FILTERS).

• Minimum variant allele frequency (VAF) of 5% (LOW_VAF).

• Minimum sequencing coverage of 20× in the tumor sample (LOW_COVERAGE_TUMOR).

• Minimum sequencing coverage of 10× in the normal sample (LOW_COVERAGE_NORMAL).

• Located inside exonic regions as defined by the set of canonical transcripts used by VEP v104 on GRCh37 assembly (NOT_EXONIC).

• Allele frequency across all gnomAD v2.1 exome subpopula­tions is <0.04% (COMMON_VARIANT). This rule is not applied for variants annotated by OncoKB (see section “Cancer Gene List and Catalog of Oncogenic Events”).

• Localized within the META-PRISM target region, which is defined by the intersection of the capture regions of all four different kits used (OFF_TARGETS_INTERSECTION). This region spans approximately 36.6 Mb.

Supplementary Fig. S12 summarizes the effect of each filter individually and in combination with other filters. A total of 117,747 somatic point mutations and small indels were used for analysis (see Data Table 3). All mutations were annotated using VEP release 104 (66) on canonical transcripts. We took advantage of the plugin feature of VEP to enhance the annotations with different metrics from the dbNSFP v4.2 database.

MSI analysis was performed with MANTIS v1.0.3 following the same procedure that was described in the original study (67). In order to run MANTIS, a list of microsatellite loci has to be compiled in a 6-column BED file, which includes not only the genomic coordinates but a column with the motif and another column with its count on the reference genome. In our analysis, we used the same BED file that was used in the TCGA study (68) and includes 2,530 loci. The method scrutinizes repetitive regions in aligned reads from the tumor and normal BAM files, each locus at a time. Per-locus read counts for each repeat motif are calculated in the tumor and normal and used for computing an instability score. Subsequently, the average of all locus instability scores is computed to produce a final score for the tumor/normal pair. Scores reported range from 0.0 (entirely stable) to 2.0 (entirely unstable). MSI-high was called in samples that had a final score exceeding the default threshold of 0.4.

SCNAs, tumor purity, and average tumor ploidy were identified with the FACETS R package v0.5.14 (29) run with parameters cval_pre = 25 and cval_pro = 500 (see Data Tables 5 and 6). For each copy-number segment identified, FACETS provides allele-specific estimations of the copy numbers, namely, the lower copy number (LCN) and total copy number (TCN).

The number of WGDs was determined to be the lowest positive integer k such that at least 11 autosomes have undergone k duplications. An autosome is identified to have undergone k duplications if the major allele ploidy (calculated as TCN-LCN) is strictly superior to 1.5 × 2^k-1 on at least half of the chromosome length. If this condition cannot be fulfilled with k = 1, we assume that no WGD has occurred. This rule mirrors the simple heuristic used in Priestley and colleagues (2) for identifying WGD-positive samples (regardless of the number k of duplications).

SCNAs in all chromosomes were then categorized into six classes (similarly to what was done by Priestley and colleagues; ref. 2) taking into account the estimated number of WGDs (0 or k ≥ 1) and the patient's gender, as detailed in Supplementary Fig. S13A.

Chromosome arm events were identified by running facetsSuite v2.0.8 on the segments generated by FACETS (see Data Table 7). Chromosomal SCNAs were then classified into either gain or loss with rules depending on the estimated number of WGDs (0 or k ≥ 1) and the patient's gender, as detailed in Supplementary Fig. S13B.

As WGD is a strong confounding factor for the calling of LLG, MLG, LOH, and cn-LOH events, the fractions of these types of SCNAs at each locus in Fig. 2E were estimated only among nearly diploid tumors, whereas high-level SCNAs were computed considering the full cohort.

Point mutations from all 569 samples with WES data available were summarized into a matrix of mutation counts with 96 rows for the 96 trinucleotide mutation categories (disregarding the strand), totaling 109,235 point mutations (see Data Table 8). Only samples with at least 50 detected somatic point mutations were used for the mutational signature analysis, resulting in the selection of 450 out of 569 samples with WES data available. Mutational signature activities in each sample were determined by projecting this matrix onto the database of SBS signatures COSMIC v3.2 (78 signatures; ref. 69) using the MutationalPatterns R package v3.2.0. We then made an effort to mirror the original algorithm used for deconvoluting mutational signatures (SigProfiler MATLAB package, also available in Python on Alexandrov's lab Github page: https://github.com/AlexandrovLab/SigProfilerExtractor). For that, we assessed the individual contribution of each signature to the mutation profile of each sample and removed signatures that contributed to less than 1% of the profile as measured by the decrease in cosine similarity by removing the signature. Mutations from removed signatures were redistributed among the other signatures using a nonnegative least-squares procedure. We used the part of the Julia code (https://bitbucket.org/bbglab/sigprofilerjulia), which reimplements the SigProfiler algorithm, from Pich and colleagues (10) to run this sparsity-inducing step (see Data Table 9).

Gene fusions were called on the GRCh38 reference genome by filtering and combining fusions predicted by Arriba v1.2.0 (70), EricScript v0.5.5 (71), Pizzly v0.37.3 (https://doi.org/10.1101/166322), and StarFusion v1.8.1 (72), all included in the nf-core/rnafusion pipeline v1.2.0.

A total of 707,027 fusions were identified by any of the four callers listed above (see Data Table 10). For each individual caller, the set of putative gene fusions was limited to exclude fusions that have been previously reported in studies of normal tissues (blacklists) and to retain only the fusions that have been reported in studies of cancer tissues (whitelists) or that involve a cancer driver partner (see section “TCGA Data—Gene Fusions” for all the details of the lists used). Additionally, only gene fusions seen by both Arriba and EricScript or by both Pizzly and StarFusion, independently of the predicted breakpoints, were retained (see section “TCGA Data—Gene Fusions” for detailed explanations about this specific combination). After filtering, 851 gene fusions (disregarding breakpoints, see Supplementary Fig. S14), detected in 445 out of 944 samples with RNA-seq data (3 samples repeatedly failed on one or more algorithms; see Supplementary Table S2), were considered. Supplementary Table S8 lists the 761 gene fusions from the 797 samples of the RNA-seq subcohort that were used in the analyses.

RNA from tumor samples was extracted using the AllPrep DNA/RNA Mini Kit (Qiagen) combined with TRIzol LS (Thermo Fisher). cDNA was synthesized with 50 ng of RNA using the High-Capacity cDNA Reverse Transcriptase Kit (Life Technologies). PCR amplification was performed using Fast Start PCR Master (Roche). We used the following program PCR amplification: hot start 95°C for 2 minutes, then 38 cycles of 95°C for 30 seconds, 58°C for 1 minute and 72°C for 1 minute. A final extension was done at 72°C for 7 minutes (VeritiPro Thermal Cycler, Thermo Fisher Scientific). Electrophoretic migration was performed in TBE1 × 2% agarose gel. PCR products were purified according to size with a QIAquick PCR purification kit for single bands or a QIAquick gel extraction kit (Qiagen) for multiple bands. Eurofins Genomics synthesized the primers and did the sequencing. Out of 18 gene fusions tested, 16 exhibited a positive PCR (Fig. 6E).

To find germinal cancer-predisposing variants, we used a conservative approach aiming to maximize the representation of variants with strong effects in our dataset. After annotation with ANNOVAR (version from April 2018), we extracted variants that met all the following criteria:

• The variant is in a curated list of 130 cancer-predisposing genes. The list was established starting from the list of 152 genes curated by Huang and colleagues (32), from which we additionally excluded 22 genes after literature curation showing that these genes have only weak support for being cancer-predisposing. These 22 excluded genes are ABCB11, ELANE, FAH, FH, GBA, GJB2, HFE, HMBS, HNF1A, ITK, MTAP, PHOX2B, PMS1, PRF1, SBDS, SERPINA1, SLC25A13, SRY, STAT3, TNFRSF6, UROD, and WAS.

• The variant has a maximum general population allele frequency of 5% in the gnomAD v2.1.1 exome database (73).

• The variant is annotated in the ClinVar database (74) as “Pathogenic,” “Likely_pathogenic,” or “Pathogenic/Likely_pathogenic” or is located in a tumor suppressor gene and has a protein-disrupting role (“splicing,” “frameshift deletion,” “frameshift insertion,” “stopgain,” or “stoploss”) except for the POLD1 and POLE tumor suppressors for which only variants annotated by ClinVar were retained.

Candidate variants were then manually reviewed in the BAM files. A total of 93 cancer-predisposing variants were detected in 73 out of 569 META-PRISM samples with germline data available and included in analyses (15 samples harbored two or more such variants; see Data Table 2).

MutPanning v2.0 (30) was run on each of the five main tumor types represented in META-PRISM with default options.

We compiled a list of 360 cancer driver genes by intersecting the list of driver genes (Tiers 1 and 2) from COSMIC census v92 and the list of genes annotated in the OncoKB database as of July 2021 (31).

Oncogenic events were identified by intersecting point mutations, indels, focal high-level gains (HLG from segments spanning less than 10 Mb), focal HD from segments spanning less than 10 Mb, and gene fusions with the OncoKB and CIViC (38) databases. OncoKB constitutes a literature- and knowledge-based database that is accessible through the Web application programming interface (API) after having registered an account and requested a token. In accordance with the type of event to be annotated, different scripts from oncokb-annotator (available at https://github.com/oncokb/oncokb-annotator) were used when operating with the OncoKB API:

• MafAnnotator.py for point mutations and indels

• CnaAnnotator.py for SCNAs

• FusionAnnotator.py for gene fusions

CIViC is also a literature and knowledge-based database, but no API was available yet at the time of our analyses. As a consequence, we have developed in-house scripts to annotate point mutations, indels, gene fusions, HLGs, and HDs using the table of clinical evidence summaries 01-Jan-2022-ClinicalEvidenceSummaries.xlsx from the January 2022 release. A minor number of errors were manually curated from the table. Additionally, missing genomic coordinates for point mutations and indels were manually filled where possible.

Importantly, the majority of OncoKB and CIViC annotations are tumor-type specific. Consequently, the annotation algorithms require an extra file containing the tumor type of each sample to allow for precise on-label annotations. OncoKB uses the Memorial Sloan Kettering Cancer Center's (MSKCC) oncotree nomenclature, whereas CIViC uses designations that do not follow specific rules and have a varying degree of specificity. A thorough work of tumor-type matching was performed by a panel of oncologists to navigate between TCGA, MSKCC's oncotree, and CIViC designations.

All HLGs, HDs, and gene fusions found in either OncoKB or CIViC databases were retained. However, point mutations and indels that intersected with the CIViC database but not the OncoKB database were discarded given that manual inspection of these CIViC-only events revealed false matching arising from an unspecific variant description in CIViC. Among the point mutations and indels that were annotated by oncokb annotator, all events with MUTATION_EFFECT as “Likely Neutral,” “Neutral,” or “Unknown” were discarded unless the ONCOGENIC field was “Likely Oncogenic” or “Predicted Oncogenic.” Lastly, only point mutations and indels that were classified as either “Missense_Mutation, “Frame_Shift_Del,” “Frame_Shift_Ins,” “In_Frame_Del,” “In_Frame_Ins,” Nonsense_Mutation,” “Splice_Site,” or “Translation_Start_Site” were retained.

After the catalog of oncogenic events was defined as described in the previous section, we took advantage of these databases to assess how these events may inform treatment decisions. OncoKB provides indications of how events may be used to inform treatment indications (sensitivity) or counterindications (resistance) on 5-level and 3-level confidence scales (version from July 2021), respectively. Similarly, CIViC provides the same type of indications but uses a 5-level confidence scale for both types. Sensitivity and resistance confidence scales were harmonized against a common 3-level reference scale following ESCAT guidelines (39).

ESCAT Tier 1 level (SOC) is for biomarkers that are used in clinical practice. These were identified as:

• CIViC Level A (proven/consensus association in human medicine)

• OncoKB Level 1 (FDA-recognized biomarker predictive of response to an FDA-approved drug in this indication)

• OncoKB Level 2 (SOC biomarker recommended by the National Comprehensive Cancer Center or other professional guidelines predictive of response to an FDA-approved drug in this indication)

• OncoKB Level R1 (SOC biomarker predictive of resistance to an FDA-approved drug in this indication)

ESCAT Tier 2 level (investigational) encompasses investigational targets that likely define a patient population that benefits from a given drug but for which additional data are needed. These were identified as:

• CIViC Level B (clinical trial or other primary patient data support association)

• OncoKB Level 3A (compelling clinical evidence supports the biomarker as being predictive of response to a drug in this indication)

• OncoKB Level R2 (compelling clinical evidence supports the biomarker as being predictive of resistance to a drug)

ESCAT Tier 3 level (hypothetical) includes all targets that have demonstrated a clinical impact on other tumor types or that are supported by scarce data (often from case reports only). These were identified as:

• CIViC Levels C, D, and E (supported by case-study, preclinical, and inferential data, respectively)

• OncoKB Level 3B (SOC or investigation biomarker predictive of response to an FDA-approved drug in another indication)

• OncoKB Level 4 (compelling biological evidence supports the biomarker as being predictive of response to a drug)

The transcriptomic profile of each tumor was used to assess the expression level in a set of 29 functional gene expression signatures representing the functional pathways or cell-type composition (immune, stromal, and other cellular populations) of the TME (16). The list of genes defining each of these signatures was downloaded from the Github repository of the article (https://github.com/BostonGene/MFP/tree/master/signatures). The single-sample gene set enrichment analysis activity score of each signature was computed from the gene-level transcript per million (TPM) tables (Kallisto-TxImport pipeline, see section 3) using the Python function “ssgsea_formula” defined in the Github (see Data Table 11). These scores were then used to classify tumors into one of the four tumor TME subtypes defined in the article.

In order to assign a TME to a new transcriptomic profile, we trained a model using TCGA scores and labels provided by the authors in their Github repository (see https://github.com/BostonGene/MFP/tree/master/Cohorts/Pan_TCGA). To begin with, we checked that we were able to reproduce the normalized signature scores for the same TCGA dataset as used by Bagaev and colleagues (16). For that purpose, signature scores of TCGA samples were recomputed starting from the TPM expression tables that we have used in this study (see section “TCGA Data—Gene Expression”) and were then median-centered and scaled by the mean-absolute deviation independently for each tumor, mirroring the method described in the article. Normalized scores were compared with the scores provided by the authors and showed very high concordance, as depicted in Supplementary Fig. S15.

We next trained a multiclass classification model on TCGA data to predict to which TME each sample belongs using the 29 normalized scores as predictors. We compared the accuracy of different machine learning models trained on the normalized scores recomputed by us and the normalized scores provided in the article. Best hyperparameters were fit using a 5-fold cross-validation procedure, and the scores reported in Supplementary Fig. S16 present the average of the five test scores from the five internal cross-validation splits.

Overall, we saw a very small decrease in classification performance when substituting the scores computed by the authors with the scores recomputed by us, in line with the high concordance described above. In the Github, the authors recommend the usage of KNeighbors classifiers with k = 35, but these models were the worst performing in all five models that we have tested [with best hyperparameters (tested values ranged from 2 to 80) at k = 37 and k = 36 for scores provided by the authors and recomputed by us, respectively]. For our analyses, we used predictions from the logistic regression model trained on the scores recomputed by us as it showed excellent performance (90% weighted F1-score) and was the easiest to interpret.

Per-tumor-type normalization parameters learned from the full TCGA dataset were reapplied when performing predictions for the transcriptomic profiles used in our analysis. This normalization procedure allows for the use of the model in a predictive setting and for comparison between different cohorts. For that reason, only tumor types that were analyzed by Bagaev and colleagues (16) could have a TME assignment. Most notably, the SARC tumor type was excluded by the authors, and therefore we did not analyze SARC tumors, nor did we analyze rare tumors.

To profile the relative activity of TFs across the samples, we filtered out lowly expressed genes from the raw read counts table (median number of reads less than 10 per gene) and used the TMM method of edgeR package (75) to normalize the counts and extract logCPM values. The logCPM values were used as input to the DoRothEA regulons package (42) based on the statistical method VIPER (76) to get the relative activity of each TF in every sample as a normalized enrichment score (nes = T) taking into account pleiotropy of the gene activity (pleiotropy = T; see Data Table 12). The DoRothEA regulons is a database of direct targets of TFs with different levels of evidence (A–E). We used only TFs with at least 10 targets with a high level of evidence (A, B, and C), which were annotated as important parts of oncogenic signaling pathways (FOXA1, AR, MYC, MYCN, GLI1, GLI2, LEF1, TCF7L1, TCF7L2, TCF7, SMAD2, SMAD3, TEAD1, TEAD2, TEAD3, TEAD4, E2F1, E2F2, FOXO1, FOXO3, FOXO4, FOXO6, STAT1, STAT2, STAT3, STAT4, STAT5A, STAT5B, STAT6, ELK1, ELK4, ETS1, ETS2, FOS, JUN, RBPJ, RBPJL, NFKB, RELA, RELB, TP53, HIF1A, ARNT, MYB; refs. 1, 77).

Patients and sample attributes from the MET500 study were downloaded from the database of Genotypes and Phenotypes study ID phs000673.v4.p1 with permission and from the Supplementary Tables S1 and S2 of Robinson and colleagues (3). All relevant variables were merged into curated tables that were used for analysis.

Raw sequencing files were downloaded with permission and processed with our internal pipeline as described in the previous sections. Only samples included from patients in the initial MET500 article were used. A total of 106,341 somatic mutations and small indels were used for analysis after applying the filtering procedure described in the section “Germline and Somatic Calling of Point Mutations and Small Indels.” The filtering that resulted in this list of somatic calls is summarized in Supplementary Fig. S17.

The “OFF_TARGETS_INTERSECTION” does not appear in Supplementary Fig. S17 because for these samples, the bed file containing only the positions in the intersection of all capture kits used on META-PRISM samples was provided as input to the –intervals parameter of Mutect2.

Raw sequencing files were processed with our internal pipeline and germline mutations were called using HaplotypeCaller as described in the previous sections. A total of 71 germline cancer-predisposing variants were detected in the 500 MET500 samples after applying the filtering procedure described in the section “Germline Cancer-Risk Variant Calling.”

Raw sequencing files were processed with our internal pipeline and MSI was called using MANTIS as described in the previous sections.

Raw sequencing files were processed with our internal pipeline, and SCNAs, WGD status, chromosome arm SCNAs, purity, and ploidy were identified using our internal SCNA pipeline based on FACETS.

RNA-seq data files were downloaded for 497 samples. The majority of samples had multiple files available produced from either polyA+ selection, hybridization capture, or both. As all RNA-seq files were produced from polyA+ in META-PRISM, we prioritized polyA+ whenever possible, resulting in 386 polyA+ and 111 hybridization capture RNA-seq libraries. Gene expressions were then quantified using the Kallisto + TxImport pipeline as used on META-PRISM samples.

Raw sequencing files were processed with our internal pipeline, and gene fusions were identified using our internal gene fusion pipeline. A total of 731 gene fusions (disregarding breakpoints; see Supplementary Fig. S18) were detected in 308 out of 497 samples with RNA-seq data (all samples were successfully processed by each fusion-calling algorithm).

Patients and sample attributes for the TCGA study were downloaded from the Genomic Data Commons (GDC) data portal using the R package GenomicDataCommons and from Supplementary Tables publicly available on the PanCanAltas page (https://gdc.cancer.gov/about-data/publications/pancanatlas). In short, attributes from all 11,315 patient samples and all 34,815 aliquots (23,346 WES and 11,469 RNA-seq) of interest for our study were downloaded. Considering the data that were available to us and the latest update of patients’ annotations, a total of 875 patients were excluded. The reasons and numbers of patients excluded for each reason are detailed in Supplementary Fig. S19A. A total of 10,440 patients for whom we had one or multiple types of molecular data available and had no reason for exclusion were considered for analysis. As detailed hereafter, data from each modality (somatic mutations, germline mutations, MSIs, SCNAs, gene expressions, and gene fusions) were limited to these 10,440 patients and further reduced to exclude samples that did not pass exclusion criteria or quality control checks (Supplementary Fig. S19B–S19F).

TCGA somatic mutation catalog was taken from the controlled-access MAF file mc3.v0.2.8.CONTROLLED.maf.gz (https://gdc.cancer.gov/about-data/publications/pancanatlas), down­loaded with permission. All filters described in the original “FILTER” column were applied as well as the following additional filters:

• Minimum VAF of 5% (LOW_VAF).

• Minimum sequencing coverage of 20× in the tumor sample (LOW_COVERAGE_TUMOR).

• Minimum sequencing coverage of 10× in the normal sample (LOW_COVERAGE_NORMAL).

• Localized within the META-PRISM target region (OFF_TARGETS_INTERSECTION).

• Located inside exonic regions as defined by the set of canonical transcripts used by VEP v104 on GRCh37 assembly (NOT_EXONIC).

• Allele frequency across all gnomAD v2.1 exome subpopulations is <0.04% (COMMON_VARIANT). This rule is not applied for variants annotated by OncoKB (see section “Cancer Gene List and Catalog of Oncogenic Events”).

• Only SNPs and multi-nucleotide polymorphisms (MNP) mutations identified by at least two of the five callers used by the MC3 consortium were retained. Likewise, only INDELs identified by Indelocator or Varscan (INDELOCATOR or VARSCANI tags in the “CENTERS” column from the controlled-access MAF file) were retained (INDEL/SNP_CALLING_ALIGNMENT).

These filters were determined after running our internal pipeline on a test set of 58 TCGA raw WES files downloaded with permission from the GDC data portal. Agreement results from the comparison of the variants detected using our pipeline with the variants reported in the mc3.v0.2.8.CONTROLLED.maf.gz file after applying the set of optimal filters are summarized in Supplementary Fig. S20A, in which DSC stands for Dice-Sorensen coefficient, TP for true positives, and FP for false positives. True mutations designate mutations from the MC3 MAF file. DSCs of 0.89 for SNPs/MNPs and 0.84 for INDELs indicate very good agreement. Supplementary Fig. S20B (resp. Supplementary Fig. S20C) provides details about the number of predicted SNPs (resp. INDEL) versus true SNPs (resp. INDEL) for each of the 58 samples analyzed.

A total of 2,109,671 somatic mutations and small indels were used for analysis in the 8,688 patients for whom somatic mutation data were available and no reason for exclusion existed as detailed in Supplementary Fig. S19B. The filtering that resulted in this list of somatic calls is summarized in Supplementary Fig. S21. All PASS mutations from the refiltered MC3 MAF file were split into individual VCF files for each tumor/normal pair and were subsequently annotated using the same annotation pipeline as used on META-PRISM data (see section “Germline and Somatic Calling of Point Mutations and Small Indels”).

The list of cancer-risk germline variants detected in TCGA samples was taken from the file PCA_pathVar_integrated_filtered_adjusted.tsv available at https://gdc.cancer.gov/about-data/publications/PanCanAtlas-Germline-AWG (32). Reduc­ing the original list of 1,393 mutations to only mutations in the 10,440 patients considered for analysis and excluding samples flagged as contaminated in GDC release notes v32.0 (Supplementary Fig. S19C) results in 1,342 pathogenic germline variants detected in 1,253 patients. After running the filtering procedure described in the section “Germline and Somatic Calling of Point Mutations and Small Indels,” 1,082/1,253 variants were retained and used for analysis.

MSI scores were retrieved from supplementary data of Bonneville and colleagues (68). Samples flagged as contaminated in GDC release notes v32.0 were excluded, resulting in available MSI score for 9,298 patients among the 10,440 TCGA patients we considered (Supplementary Fig. S19D).

SCNAs, WGD status, chromosome arm SCNAs, purity, and ploidy were identified by running, with permission, our internal SCNA pipeline on all available TCGA WES files (21,987 BAM files from 10,332 patients with gender available and at least one pair tumor/normal WES files; 12,129 pairs in total). Samples flagged as contaminated in GDC release notes v32.0, that failed FACETS, or that failed quality control of FACETS profiles were excluded, resulting in available SCNAs for 9,570 patients (Supplementary Fig. S19E). This reanalysis was a necessary step to avoid large batch effects arising when comparing data generated from different copy-number variation assessment techniques. Execution of the pipeline was performed via the Google Cloud Engine with support from the Institute for Systems Biology Cancer Genomics Cloud.

Gene and transcript expression tables for all TCGA RNA-seq samples were retrieved from the supplementary data of Zheng and colleagues (78). This pipeline is identical to the pipeline that we have used on META-PRISM samples, thereby minimizing technical differences between the cohorts. Samples flagged with quality control warnings in GDC aliquot-level notes were excluded, resulting in available RNA-seq profiles for 9,298 patients among the 10,440 TCGA patients we considered (Supplementary Fig. S19F).

Three independent and publicly available lists of TCGA gene fusions were retrieved from the following sources:

• Hu and colleagues (PRADA; ref. 79) Supplementary Table nar-02671-data-e-2017-File007.xlsx

• Gao and colleagues (StarFusion; ref. 80) Supplementary Table S1

• Dehghannasiri and colleagues (DEEPEST; ref. 81) Supplementary Table pnas.1900391115.sd01.xlsx

In order to adjust the filtering criteria and compare the nf-core fusion pipeline (see section “Identification of Gene Fusions”) used for META-PRISM fusions to the three different pipelines used for the aforementioned lists of TCGA fusions, we downloaded with permission RNA-seq FASTQ files for 69 TCGA samples from the GDC data portal and analyzed them using the nf-core fusion pipeline. All 69 samples were successfully processed by the six different callers available in the nf-core fusion pipeline, namely, Arriba (“AR,” v1.2.0), EricScript (“ES,” v0.5.5), FusionCatcher (“FC,” v1.20), Pizzly (“PZ,” v0.37.3), Squid (“SQ,” v1.5), and StarFusion (“SF,” v1.8.1). All detected fusions were then annotated with FusionAnnotator that connects with datasets of fusions detected in cancer or normal tissues (see https://github.com/FusionAnnotator).

Supplementary Fig. S22 summarizes the overlap between the fusions detected by these six algorithms. The fusions predicted by Squid were very different from that predicted by any of the five other callers (Supplementary Fig. S22) and were therefore removed from all analyses. Additionally, even though FusionCatcher was successful on all 69 TCGA samples, it repeatedly failed on some META-PRISM samples and was therefore not considered further. Consequently, only fusion calls from Arriba, EricScript, Pizzly, and StarFusion were used.

As our study only aimed at describing the variations relevant to cancer, only fusions known in cancer or involving at least one oncogenic partner (see section “Cancer Gene List and Catalog of Oncogenic Events”) were analyzed. We therefore limited the lists of fusions as detailed hereafter.

First, we removed all fusions that have been previously reported in studies of normal tissues. More specifically, fusions were removed if they met any of the following criteria:

• Are in Babiceanu_Normal list (82).

• Are in ChimerSeq_Normal_v4.0 list (available upon request to the authors; ref. 35), which was established from the analysis of 1,144 TCGA normal samples and curated in order to remove well-known fusions (e.g., TMPRSS2–ERG) sometimes seen in normal samples.

• Are in GTEX_V6 Supplementary Table S3 (83).

• Have at least one “Red Herring” flag (FusionAnnotator annotations) among the following: GTEx_recurrent_StarF2019, BodyMap, DGD_PARALOGS, HGNC_GENEFAM, Greger_Normal, ConjoinG.

• One of the partners is not protein coding.

Second, only fusions that met one of the following criteria were retained:

• Are in COSMIC v95 list of fusions (see https://cancer.sanger.ac.uk/cosmic/fusion).

• Are in ChimerKB v4.0 list (35).

• Are in Chitars Cancer v5.0 list (ref. 36; available upon request to the authors).

• Are in TIC v3.3 list (see https://genetica.unav.edu/TICdb/).

• One of the partners is a cancer driver (see section “Cancer Gene List and Catalog of Oncogenic Events”).

Lastly, after filtering the fusion calls for the 69 samples in each of the three TCGA published lists and in each of the lists of fusions predicted by the four callers considered, we looked for the combinations of calls that showed the best agreement. The best agreement was obtained between the following combinations:

• For the three published TCGA lists: fusions seen by StarFusion or by both DEEPEST and PRADA

• For the four callers in our pipeline: fusions seen by both Arriba and EricScript or by both Pizzly and StarFusion

Agreement results are summarized in Supplementary Fig. S23, in which DSC stands for Dice-Sorensen coefficient, TP for true positives, and FP for false positives. True fusions designate fusions from the combination of published TCGA lists. The DSC was 0.90 if fusion calls were assessed to be concordant regardless of the breakpoint prediction and was 0.77 if the predicted breakpoint was required to be identical.

Applying the above rules, we retained TCGA fusions seen by StarFusion alone (8,194 fusions), by DEEPEST and PRADA exactly (2,604 fusions), or by all three methods (9,989 fusions). Supplementary Fig. S24 details the overlap between these three gene fusion lists and highlights in red the fusions that we have retained.

ExAC release 1.0 vcf file was downloaded from the gnomAD download page (https://gnomad.broadinstitute.org/downloads#exac-variants) and annotated with ANNOVAR using the same pipeline as used for META-PRISM vcf files. The same filtering procedure as described in the section “Germline Cancer-Risk Variant Calling” was applied, resulting in 2,002 cancer-predisposing germline variants among the non-Finnish European population.

The tumor purity of DNA samples for META-PRISM, MET500, and TCGA tumors was estimated using FACETS. We observed a very good tumor purity in all three cohorts, with an average purity above 50% in each. Moreover, the average tumor purity of META-PRISM samples was lower than that of TCGA, which rules out the possibility that the observed increased variation in META-PRISM is due to better tumor purity.

All univariate analyses are based on Kaplan–Meier estimates comparing two or more groups of patients defined by values of a categorical variable. The survival curves and log-rank tests comparing these groups of patients were all obtained using the survminer R package. Univariate survival analyses were applied to assess the association with the survival of (i) the most frequently altered driver genes (only genes altered in at least 12 patients if considering only one tumor type, 20 patients otherwise); (ii) the TME subtypes; (iii) the number of metastatic sites and the number of resistant drugs before biopsy; and (iv) the GRIM score in the five most represented tumor types in our cohort (BLCA, BRCA, LUAD, PAAD, PRAD) and in the META-PRISM WES and RNA-seq subcohort (10 tumor types).

All multivariate analyses are based on Cox multivariate regressions incorporating parts of all possible covariates depending on the effects that we were trying to measure. To provide a comparative assessment of the added and independent value of biomarkers derived from WES or RNA-seq as compared with clinical biomarkers (age, gender, global tumor characteristics, blood test results, treatment history), we experimented with multiple models that incorporated different combinations of predictors. The baseline model (M2) to be improved upon was composed of all standard clinical variables plus the continuous components of the GRIM score (LDH, albumin, neutrophil, and lymphocyte levels). A separate model was run for each of the five most represented tumor types in our cohort (BLCA, BRCA, LUAD, PAAD, and PRAD) and the META-PRISM WES and RNA-seq cohort (10 tumor types). In order to ensure a fair comparison with more complex models, only samples with complete molecular profiles (WES and RNA-seq) were used in the different Cox models.

For each combination of predictors and samples, we repeated the feature selection and coefficient estimation steps 1,000 times to assess the selection procedure's stability, if any, and to provide robust estimates and CIs for the effect of each variable on survival. Each of the 1,000 repeats consisted in running the selection procedure, if any, and training the model on a random 80% subsample. The remaining 20% were used to assess the model quality (C-index). The C-index values reported in Fig. 7B are averages of the 1,000 estimates of the C-index on the 20% test subsamples. Estimates of the covariates’ coefficients were computed by averaging the coefficient fitted on each 80% subsample across the 1,000 repeats. In case a selection procedure is active at each repeat, covariates were assigned an estimated coefficient of zero every time the selection procedure did not retain them. Confidence intervals at the 95% level were estimated by computing the empirical 2.5% and 97.5% quantiles from the 1,000 estimates.

A selection procedure was applied in case the number of covariates was unreasonably high given the number of available observations. We experimented with univariate–fit univariate Cox model for each candidate predictor and selected only predictors for which adjusted P values fell below a certain threshold—and multivariate selection procedures—fit Cox model with lasso penalization and selected only predictors with nonzero coefficients to identify a small set of predictors that, hopefully, was a superset of all important predictors. The Cox regressions on predictors selected by the lasso were always more prognostic (higher value of C-index) than when predictors were selected through the univariate procedure.

Each set of 1,000 repeats of Cox regressions was preceded by preprocessing steps (run once for each combination of samples and features) in which covariates were formatted, imputed, min–max transformed, and analyzed for redundancy. The imputation relied on the MICE R package (84) that performs multiple imputations of the same dataset by iteratively learning a predictive model of each covariate with missing data from all other covariates and randomly sampling from observed values through the predictive mean matching procedure. Data tables with at least one missing value were imputed five times, and coefficient estimates across each of the 1,000 repeats were calculated by averaging estimates across the five imputed tables.

The C-index was estimated in two different ways, either using the C-statistic proposed by Harrell and colleagues (85) or using the inverse probability weighting estimate presented by Uno and colleagues (86) to account for the fact that the C-statistic proposed by Harrell and colleagues depends on the censoring distribution that in practice is rarely independent of the covariates used in the Cox regression. As survival times are subject to censoring, and because we were mainly interested in discriminating between patients having a survival longer or shorter than 6 months, we used an estimate of the C-index truncated at 6 months by considering only pairs of patients for which one of the two patients had a survival shorter than 6 months. The estimator proposed by Uno and colleagues (see formula 2.3 from their work), as well as the C-statistic proposed by Harrell and colleagues, was implemented in C using a truncation time of 6 months. A comparison of the estimated values from both estimators showed very little difference due to the fact that only a minority of patients were right-censored in our cohort (<10%). Values reported in the main text use the C-statistic from Harrell and colleagues truncated at 6 months.

For each analysis in which the frequency of a binary event was compared at the cohort level, we used the stratified Cochran–Mantel–Haenszel test stratifying by tumor type to control for different tumor-type composition between the cohorts. In cases in which we analyzed the distribution of a continuous variable, the cohort-level comparison was made through a standard Mann–Whitney U test regardless of the tumor-type composition.

For analyses performed for each individual tumor type, we used the Fisher–Boschloo test when comparing the frequency of a binary event (e.g., chromosome arm gains/losses, WGD status, mutation status) and the Mann–Whitney U test when comparing the medians of a continuous variable (TMB, number of driver events, number of mutations contributed by each mutational signature, genome fraction covered by each type of SCNA).

Log odd ratio estimates and CIs measuring the associations between one or multiple predictors and a binary outcome were computed by fitting standard logistic regression models (statsmodels python package and stats R package).

Comparisons of Kaplan–Meier survival curves were performed using the standard log-rank test implemented in the ggsurvplot and survdiff R packages. Adjusted HRs were estimated by multivariate Cox regression modeling incorporating the variable of interest and the variable to adjust for (e.g., absence/presence of oncogenic alterations adjusted for tumor type in Fig. 7).

Whenever multiple independent tests were performed within a given analysis, we applied the Benjamini–Hochberg multiple testing correction controlling for the FDR and only reported adjusted P values < 0.05 (so-called q-values) unless explicitly mentioned otherwise.

All the codes supporting the analysis, numbers, and figures presented in this study are available at https://github.com/gustaveroussy/MetaPRISM_Public. The code is also made available on CodeOcean at https://codeocean.com/capsule/2014781/tree/v1.

A list of all external data sources and files used for this study is available in Supplementary Table S11.

Experimental data generated in this study have been deposited to the European Genome-phenome Archive under the accession number EGAD00001009684. Please refer to the forms and README file from https://github.com/gustaveroussy/MetaPRISM_Public/tree/master/data for instructions on how to access the data.

Standardized data tables are made publicly available on the study page of Gustave Roussy's cBioPortal available at https://cbioportal.gustaveroussy.fr/study/summary?id=metaprism_2023. Additionally, 12 data tables (Data Tables 1–12) of processed sequencing data are made publicly available on Gustave Roussy's nextcloud server. The link is included in the description of the study on the cBioPortal page.

All the other data supporting the findings of this study are available within the article and its Supplementary Information files or from the corresponding authors upon reasonable request. A reporting summary for this article is available as a Supplementary Information file.

L. Verlingue reports personal fees from Adaptherapy, grants from Bristol Myers Squibb, and other support from Pierre Fabre and Servier outside the submitted work. S. Michiels reports personal fees from Amaris, Roche, IQVIA, Sensorion, Biophytis, Servier, and Yuhan outside the submitted work. A. Hollebecque reports grants from AstraZeneca during the conduct of the study, as well as personal fees from Basilea, Debiopharm, Eisai, Incyte, QED Therapeutics, Relay Therapeutics, Servier, and Taiho and grants from Incyte outside the submitted work. L. Mezquita reports research grant/funding from Bristol Myers Squibb, Boehringer Ingelheim, Amgen, Stilla, Inivata, and AstraZeneca; advisory/consultancy for Roche and Takeda; educational activities: Bristol Myers Squibb, Takeda, Foundation Medicine, and Janssen; and travel/accommodation/expenses from Roche, Takeda, Bristol Myers Squibb, and Janssen outside the submitted work. Y. Loriot reports grants from Sanofi and Janssen during the conduct of the study, as well as personal fees from Janssen and Bristol Myers Squibb, personal fees and nonfinancial support from Astellas, MSD, Gilead, Pfizer, and Merck KGaA, and grants, personal fees, and nonfinancial support from Roche outside the submitted work. B. Besse reports grants from 4D Pharma, AbbVie, Amgen, AstraZeneca, BeiGene, Blueprint Medicines, Celgene, Cergentis, Chugai Pharmaceutical, Da Voltera, Daiichi Sankyo, Eli Lilly, Ellipse Pharma, Eisai, F-Star, Genmab, Genzyme Corporation, GSK, Hedera Dx, Inivata, Ipsen, Janssen, MSD, Onxeo, OSE Immunotherapeutics, Pfizer, Pharmamar, Roche/Genentech, Sanofi, Socar Research, Taiho Oncology, Takeda, Tolero Pharmaceuticals, and Turning Point Therapeutics during the conduct of the study. L. Friboulet reports grants from Debiopharm, Incyte, Relay Therapeutics, Sanofi, and Nuvalent and nonfinancial support from Illumina outside the submitted work. F. Andre reports grants from Roche, Lilly, Novartis, Pfizer, Daiichi Sankyo, and AstraZeneca outside the submitted work. No disclosures were reported by the other authors.

Y. Pradat: Resources, data curation, software, formal analysis, investigation, visualization, methodology, writing–original draft, writing–review and editing. J. Viot: Resources, data curation, investigation, writing–review and editing. A.A. Yurchenko: Data curation, formal analysis, visualization. K. Gunbin: Software, formal analysis, visualization. L. Cerbone: Resources, data curation, investigation, writing–review and editing. M. Deloger: Resources, software. G. Grisay: Resources, data curation, investigation. L. Verlingue: Conceptualization, data curation, investigation. V. Scott: Validation. I. Padioleau: Resources, software. L. Panunzi: Resources, software. S. Michiels: Resources, methodology, writing–review and editing. A. Hollebecque: Resources, writing–review and editing. G. Jules-Clement: Resources, software. L. Mezquita: Resources, investigation. A. Laine: Formal analysis. Y. Loriot: Resources, data curation, investigation, writing–review and editing. B. Besse: Resources, investigation. L. Friboulet: Resources, data curation, validation, investigation, writing–review and editing. F. Andre: Conceptualization, resources, supervision, funding acquisition, investigation, project administration, writing–review and editing. P.-H. Cournede: Conceptualization, supervision, funding acquisition, methodology, project administration, writing–review and editing. D. Gautheret: Resources, supervision, investigation, methodology, project administration, writing–review and editing. S.I. Nikolaev: Conceptualization, supervision, funding acquisition, investigation, writing–original draft, project administration, writing–review and editing.

This work is supported by Prism—National Precision Medicine Center in Oncology funded by the France 2030 programme and the French National Research Agency (ANR) under grant number ANR-18-IBHU-0002. S.I. Nikolaev was supported by grants from Foundation ARC 2017, Foundation Gustave Roussy, and The French National Cancer Institute (RPT21145LLA). The authors are thankful to Fabian Seidl, Bill Longabaugh, and David Pot from the Institute for Systems Biology Cancer Genomics Cloud for providing technical and financial support in the analysis of TCGA WES files.

The publication costs of this article were defrayed in part by the payment of publication fees. Therefore, and solely to indicate this fact, this article is hereby marked “advertisement” in accordance with 18 USC section 1734.