Pan-cancer whole-genome comparison of primary and metastatic solid tumours – Nature

0
37


Cohort gathering and processing

We have matched tumour–normal WGS data from patients with cancer from two independent cohorts: Hartwig and PCAWG. A detailed description of the Hartwig and PCAWG cohort gathering and processing as well as comprehensive documentation of the PCAWG sample reanalysis with the Hartwig somatic pipeline is described in the Supplementary Note 1.

Tumour clonality analysis

Each mutation in the .vcf files was given a subclonal likelihood by PURPLE. Following PURPLE guidelines, we considered mutations with subclonal scores equal or higher than 0.8 to be subclonal and mutations below the 0.8 threshold to be clonal. For each sample, we then computed the average proportion of clonal mutations by dividing the number of clonal mutations by the total mutation burden (including SBS, multinucleotide variants and IDs). Finally, for each cancer type, we used Mann–Whitney test to assess the significance of the clonality difference between the primary and metastatic tumours. P values were adjusted for false discovery rate (FDR) using the Benjamini–Hochberg procedure. An adjusted P < 0.05 was deemed as significant.

In addition, we leveraged biopsy site data in patient reports to further investigate differences in metastatic tumour clonality according to the metastatic biopsy site (see also Supplementary Note 2). If the metastatic biopsy site was in the same organ or tissue as the primary tumour, we considered them as ‘local’, whereas if the metastatic biopsy site was reported in the lymphoid system or other organs or tissues, they were dubbed as ‘lymph’ and ‘distant’, respectively. Cancer types for which there was a minimum of five samples available for each of the biopsy groups were selected and Mann–Whitney test was used to compare the clonality between the biopsy groups.

Karyotype

Chromosome arm level and genome ploidy was estimated as previously described20.

First, for each chromosome arm, tumour purity and ploidy-adjusted copy number (CN) segments (as determined by PURPLE) were rounded to the nearest integer. Second, arm coverage of each integer CN was calculated as the fraction of chromosome arm bases with the specific CN divided by the chromosome arm length (for example, 60% of all chromosome 5p segments have a CN of 2, 30% have a CN of 1 and 10% have a CN of 3). We defined the arm-level ploidy level as the CN with the highest coverage across the whole arm (in the example above it would be 2). Third, we computed the most recurrent chromosome arm ploidy levels across all chromosome arms per sample (that is, observed genome ploidy).

Next, we estimated the true genome ploidy by taking WGD status (given by PURPLE) into account. If a sample did not undergo WGD, its total expected genome ploidy was deemed to be 2n. If a sample did undergo WGD and its observed genome ploidy was less than six, the estimated genome ploidy was deemed to be 4n, and 8n if the observed genome ploidy was six or more. An observed genome CN of more than eight was not found in our dataset.

Then, for each chromosome arm in each sample, we defined the normalized arm ploidy as the difference between the arm-level ploidy level and the expected genome ploidy. The resulting value was classified as 1 for differences higher than or equal to 1 (representing arm gains), as −1 for differences lower than or equal to −1 (representing arm losses) or as 0 (no difference). Normalized arm ploidy values were averaged across all samples from a cancer type in a cohort-specific manner (that is, separating primary and metastatic samples). A Mann–Whitney test was performed per cancer type and chromosome arm to assess the mean difference in arm gains or losses at the cancer-type level. The resulting P value was FDR adjusted across all arms per cancer type. Finally, q < 0.01 and a normalized arm ploidy difference higher than 0.25 was deemed to be significant.

Genomic instability indicators

To compare the differences in aneuploidy scores and the LOH proportions in each group, a Mann–Whitney test was performed per cancer type. The aneuploidy score represents the number of arms per tumour sample that deviate from the estimated genome ploidy as previously described20. The LOH score of a given sample represents the sum of all LOH regions divided by the GRCh37 total genome length. A genomic region is defined as LOH when the minor allele CN < 0.25 and major allele CN ≥ 0.8.

To compare the fraction of samples with a driver mutation in TP53 as well as the fraction of WGD samples per cohort, a Fisher’s exact test was performed per cancer type. Any TP53 driver alteration (non-synonymous mutation, biallelic deletion and homozygous disruption) was considered in the analysis. Multiple driver mutations per sample in a single gene were considered as one driver event. WGD was defined as present if the sample had more than 10 autosomes with an estimated chromosome CN of more than 1.5. P values were FDR corrected across all cancer types. A q < 0.01 was deemed to be significant for all statistical tests.

Mutational signature analysis

Signature extraction

The number of somatic mutations falling into the 96 SBS, 78 DBS and 83 ID contexts (as described in the COSMIC catalogue51; https://cancer.sanger.ac.uk/signatures/) was determined using the R package mutSigExtractor (https://github.com/UMCUGenetics/mutSigExtractor, v1.23).

SigProfilerExtractor (v1.1.1) was then used (with default settings) to extract a maximum of 21 SBS, 8 DBS and 10 ID de novo mutational signatures. This was performed separately for each of the 20 tissue types that had at least 30 patients in the entire dataset (aggregating primary and metastatic samples; see Supplementary Table 3). Tissue types with less than 30 patients as well as patients with metastatic tumours with unknown primary location type were combined into an additional ‘Other’ group, resulting in a total of 21 tissue-type groups for signature extraction. To select the optimum rank (that is, the eventual number of signatures) for each tissue type and mutation type, we manually inspected the average stability and mean sample cosine similarity plot output by SigProfilerExtractor. This resulted in 440 de novo signature profiles extracted across the 21 tissue-type groups (Supplementary Table 3). Least squares fitting was then performed (using the fitToSignatures() function from mutSigExtractor) to determine the per-sample contributions to each tissue-type-specific de novo signature.

Aetiology assignment

The extracted de novo mutational signatures with high cosine similarity (≥0.85) to any reference COSMIC mutational signatures with known cancer-type associations51 were labelled accordingly (288 de novo signatures matched to 57 COSMIC reference signatures).

For the remaining 152 unlabelled de novo signatures, we reasoned that there could be one or more signatures from one cancer type that is highly similar to those found in other tissue types, and that these probably represent the same underlying mutational process. We therefore performed clustering to group likely equivalent signatures. Specifically, the following steps were performed:

  1. (1)

    We calculated the pairwise cosine distance between each of the de novo signature profiles.

  2. (2)

    We performed hierarchical clustering and used the base R function cutree() to group signature profiles over the range of all possible cluster sizes (minimum number of clusters = 2; maximum number of clusters = number of signature profiles for the respective mutation type).

  3. (3)

    We calculated the silhouette score at each cluster size to determine the optimum number of clusters.

  4. (4)

    We grouped the signature profiles according to the optimum number of clusters. This yielded 27 SBS, 7 DBS and 8 ID de novo signature clusters (see Supplementary Table 3).

For certain de novo signature clusters, we could manually assign the potential aetiology based on their resemblance to signatures with known aetiology described in COSMIC51, Kucab et al.26 and Signal (access date 1 February 2023)52. Some clusters were an aggregate of two known signatures, such as SBS_denovo_clust_2, which was a combination of SBS2 and SBS13, both linked to APOBEC mutagenesis. Other clusters had characteristic peaks of known signatures, such as DBS_denovo_clust_4, which resembled DBS5 based on having distinctive CT>AA and CT>AC peaks. Finally, DBS_denovo_clust_1 was annotated as a suspected POLE mutation and MMRd, as samples with high contribution (more than 150 mutations) of this cluster are frequently microsatellite instable (MSI) or have POLE mutations. Likewise, DBS_denovo_clust_2 was annotated as a suspected MMRd as the aetiology, as samples with high contribution (more than 250 mutations) of this cluster were all MSI. See Supplementary Table 3 for a list of all the manually assigned aetiologies.

After applying the aetiology assignment, the de novo extraction resulted in 69 SBS, 13 DBS and 18 ID representative mutational signatures (Supplementary Table 3). Most of these (42 of 69 SBSs, 7 of 13 DBSs and 8 of 18 IDs) mapped onto the well-described mutational signatures in human cancer35,53.

Comparing the prevalence of mutational processes between primary and metastatic cancer

We then compared the activity (that is, the number of mutations contributing to) of each mutational process between primary and metastatic tumours. For each sample, we first summed the contributions of signatures of the same mutation type (that is, SBS, DBS or ID) with the same aetiology, hereafter referred to as ‘aetiology contribution’. Per cancer type and per aetiology, we performed two-sided Mann–Whitney tests to determine whether there was a significant difference in aetiology contribution of primary and metastatic tumours. Per cancer type and per mutation type, we used the p.adjust() base R function to perform multiple testing correction using Holm’s method. Next, we added a pseudocount of 1 to the contributions (to avoid dividing by 0) and calculated the median contribution log2 fold change, that is, log2((median contribution in metastatic tumours + 1)/(median contribution in primary tumours + 1)). We considered the aetiology contribution between primary and metastatic tumours to be significantly different when q < 0.05, and log2 fold change ≥ 0.4 or log2 fold change ≤ −0.4 (= ± ×1.4).

Relative contribution

Relative aetiology contribution was calculated by dividing aetiology contribution by the total contribution of the respective mutation type (that is, SBS, DBS or ID). To determine the significant difference in relative aetiology contribution, we performed two-sided Mann–Whitney tests as described above. We also calculated the median difference in contribution (that is, median relative contribution in metastatic tumours − median relative contribution in primary tumours). We considered the relative aetiology contribution between primary and metastatic tumours to be significantly different when q < 0.05 and median difference was 0.01 or more.

We also determined whether there was an increase in the number of samples with high aetiology contribution (that is, hypermutators) in metastatic versus primary cohorts. For each signature, a sample was considered a hypermutator if the aetiology contribution was 10,000 or more for SBS signatures, 500 or more for DBS signatures or 1,000 or more for ID signatures. For each cancer type, for each aetiology, we performed pairwise testing only for cases in which there were five or more hypermutator samples for either metastatic or primary tumours. Each pairwise test involved calculating P values using two-sided Fisher’s exact tests, and effect sizes by multiplying Cramer’s V by the sign of the log2(odds ratio) to calculate a signed Cramer’s V value that ranges from −1 to +1 (indicating enrichment in primary or metastatic, respectively). We then used the p.adjust() base R function to perform multiple testing correction using Bonferroni’s method.

SBS1–age correlations in primary and metastatic tumours

To count the SBS1 mutations, we relied on the definition from ref. 54 that is based on the characteristic peaks of the COSMIC SBS1 signature profile: single-base CpG > TpG mutations in NpCpG context. To ensure that these counts and the downstream analyses are not affected by differential APOBEC exposure in primary and metastatic cohorts, we excluded CpG > TpG in TpCpG, which is also a characteristic peak in the COSMIC SBS2 signature profile. In addition, for skin melanoma, CpG > TpG in [C/T]pCpG, which overlaps with SBS7a, was excluded. To obtain the SBS5 and SBS40 counts, we relied on their exposures derived from the mutational signature analyses performed in this study (described above).

To assess the correlation between SBS1 burden and the age of the patient, at biopsy we performed a cancer-type and cohort-specific linear regression (that is, separate regression for primary and metastatic tumour samples). To avoid spurious effects caused by hypermutated tumours, samples with a TMB greater than 30,000 as well as those with SBS1 burden greater than 5,000 were excluded.

For each cancer type and cohort, we then computed 100 independent linear regressions by randomly selecting 75% of the available samples. We selected the median linear regression (based on the regression slope) as representative regression for further analyses. Similarly, confidence intervals were derived from the 1st and 99th percentile of the computed regressions.

To evaluate the significance of the differences between primary and metastatic representative linear regressions (hereafter referred to as linear regression for simplicity), we first filtered out cancer types that failed to show a positive correlation trend between SBS1 burden and age at biopsy in both primary and metastatic tumours (that is, Pearson’s correlation coefficient of primary and metastatic regression greater than 0.1). Next, for each selected cancer type, we computed the regression residuals of primary and metastatic SBS1 mutation counts using, in both cases, the primary linear regression as baseline. The primary and metastatic residual distributions were then compared using a Mann–Whitney test to evaluate significance. Cancer types with a Mann–Whitney P < 0.01 were deemed as significant. Finally, to ensure that the differences were uniform across different age ranges (that is, not driven by a small subset of patients), we only considered significant cancer types in which the metastatic linear regression intercept was higher than the primary intercept.

SBS5/SBS40 correlations were computed following the same procedure and using the sum of SBS5 and SBS40 exposures for each tumour sample. If none of the mutations were attributed to SBS5/SBS40 mutational signatures, the aggregated value was set to zero. In the ploidy-corrected analyses, we divided the SBS1 mutation counts (and SBS5/SBS40 mutation counts for the SBS5/SBS40 ploidy-corrected regression, respectively) by the PURPLE-estimated tumour genome ploidy.

For each cancer type, the mean fold change (fc) was defined as \(\underline{{\rm{fc}}}=\frac{1}{40}{\sum }_{i=40}^{80}\frac{{{\rm{M}}{\rm{Pred}}}_{i}}{{{\rm{P}}{\rm{Pred}}}_{i}}\) where MPredi and PPredi are the estimated number of SBS1 mutations for a given age ith according to the metastatic and primary linear regressions, respectively. Similarly, the mean estimated SBS1 burden difference (SBS1diff) was defined as: \(\underline{{\rm{SBS1diff}}}=\frac{1}{40}{\sum }_{i=40}^{80}{{\rm{MPred}}}_{i}-{{\rm{PPred}}}_{i}\).

Clonality of clock-like mutations

SBS1 individual mutations were identified as described in the previous section. For SBS5 and SBS40 mutations, we used a maximum likelihood approach to assign individual mutations to the SBS5 and SBS40 mutational signatures in a cancer-type-specific manner.

For every SBS1 (and SBS5/SBS40 mutation), we then assign the clonality according to the PURPLE subclonal likelihood estimation, in which only mutations with subclonal (SUBCL) likelihood ≥ 0.8 were considered as such (see above).

For each tumour sample, the SBS1 clonality ratio (or respectively SBS5/SBS40 clonality ratio) was defined as the ratio between the proportion of clonal SBS1 mutations (\(\frac{{\rm{SBS\; 1\; clonal\; mutations}}}{{\rm{SBS\; 1\; mutations}}}\)) divided by the total proportion of clonal alterations in the sample (\(\frac{{\rm{Total\; clonal\; mutations}}}{{\rm{Total\; mutations}}}\)).

Primary SBS1 mutation rate and metastatic SBS1 age-corrected enrichment

We computed for each primary cancer type the average number of SBS1 per year as the number of SBS1 mutations divided by the age of the patient at biopsy (only considering primary samples and excluding hypermutated samples as described above). We then used a Spearman’s correlation to assess its association with the estimated mean SBS1 mutation rate fold change in metastatic tumours (see above). In addition, to exclude potential biases in our primary cohort, we repeated the same analysis relying on an independent measurement of primary cancer SBS1 yearly accumulation. Specifically, we used the best-estimated accumulation of SBS1 per year from ref. 30 (Supplementary Table 6) and regressed it to the fold-change estimates for the matching cancer types present in both datasets.

SV analysis

Definitions of SV type

LINX55 chains one or more SVs and classifies these SV clusters into various event types (‘ResolvedType’). We defined deletions and duplications as clusters with a ResolvedType of ‘DEL’ or ‘DUP’ whose start and end breakpoints are on the same chromosome (that is, intrachromosomal). Deletions and duplications were split into those less than 10 kb and 10 kb or more in length (small and large, respectively), based on observing bimodal distributions in these lengths across cancer types (Extended Data Fig. 5b). We defined complex SVs as clusters with a ‘COMPLEX’ ResolvedType, an inversion ResolvedType (including: INV, FB_INV_PAIR, RECIP_INV, RECIP_INV_DEL_DUP and RECIP_INV_DUPS) or a translocation ResolvedType (including: RECIP_TRANS, RECIP_TRANS_DEL_DUP, RECIP_TRANS_DUPS, UNBAL_TRANS and UNBAL_TRANS_TI). Complex SVs were split into those with less than 20 and 20 or more SVs (small and large, respectively), based on observing similar unimodal distributions in the number SVs across cancer types whose tail begins at approximately 20 breakpoints (Extended Data Fig. 5b). Finally, we defined long interspersed nuclear element insertions (LINEs) as clusters with a ResolvedType of ‘LINE’. For each sample, we counted the occurrence (that is, SV burden) of each of the seven SV types described above. In addition, we determined the total SV burden by summing counts of the SV types.

Comparing SV burden between primary versus metastatic cancer

We then compared the SV-type burden between primary versus metastatic tumours as shown in Fig. 3a. First, we performed two-sided Mann–Whitney tests per SV type and per cancer type to determine whether there was a statistically significant difference in SV-type burden between primary versus metastatic. The Bonferroni method was used for multiple testing correction on the P values from the Mann–Whitney tests (to obtain q values). Next, we calculated relative enrichment as follows: log10(median SV-type burden in metastatic tumours + 1) − log10(median SV-type burden in primary tumours + 1); and calculated fold change as follows: (median SV-type burden in metastatic tumours + 1) / (median SV-type burden in primary tumours + 1). When calculating relative enrichment and fold change, the pseudocount of 1 was added to avoid the log(0) and divide by zero errors, respectively. Fold changes are displayed with a ‘>’ in Fig. 3a when the SV burden for primary tumours is 0 (that is, when a divide by 0 would occur without the pseudocount). We considered the SV-type burden between primary versus metastatic to be significant when: q < 0.05, and fold change ≥ 1.2 or fold change ≤ 0.8

Identifying features associated with increased SV burden in metastatic cancer

To identify the features that could explain increased SV burden, we correlated SV burden with various tumour genomic features. This included: (1) genome ploidy (determined by PURPLE); (2) homologous recombination deficiency (determined by CHORD33) and MSI (determined by PURPLE) status; (3) the presence of mutations in 345 cancer-associated genes (excluding fragile site genes that are often affected by CN alterations5), hereafter referred to as ‘gene status’; and (4) treatment history, including the presence of radiotherapy, the presence of one of the 79 different cancer therapies as well as the total number of treatments received. All primary samples and all metastatic samples without treatment information were considered to have no treatment. Genome ploidy and total number of treatments received were numeric features, whereas all of the remaining were boolean (that is, true or false) features. In total, there were 429 features.

SV-type burden was transformed to log10(SV-type burden + 1) and was correlated with the 429 features using multivariate linear regression models (LMs). This was performed separately for each of the seven SV types, and for each cancer type (or subtype). In the SV main analysis (Fig. 4b–f), there were 23 cancer types, resulting in a total of 161 (23 cancer types × 7 SV types) LM models.

Each LM model (that is, per SV type and cancer type) involved training of three independent LMs with (1) both metastatic and primary samples (primary + metastatic), (2) only Hartwig samples (metastatic only), and (3) only PCAWG samples (primary only). This was done to filter out correlations between features and increased SV-type burden solely due to differences in feature values between primary and metastatic tumours. We then required features that positively correlated with SV-type burden in the primary + metastatic LM to independently show the same association in the metastatic-only or primary-only LMs. Only genomic features that independently showed positive correlation with the SV burden were further considered as significant (that is, represented in the lollipop plots).

Each of the three LMs was trained as follows:

  1. (1)

    Remove boolean features with too few ‘true’ samples.

  1. (i)

    For the primary + metastatic LM, remove gene status features with less than 15 ‘true’ samples.

  2. (ii)

    For the metastatic-only and primary-only LMs, remove gene status features with less than 10 ‘true’ samples.

  3. (iii)

    For the remaining boolean features, remove features with less than 5% ‘true’ samples.

  1. (2)

    Fit a LM using the lm() base R function to correlate log10(SV-type burden + 1) versus all features.

For each LM analysis, we used the following filtering criteria to identify the features that were correlated with increased SV-type burden:

  1. (1)

    Only keep LM analyses for which there was significant increase in SV-type burden for the respective cancer type (P < 0.01 as described in the previous section ‘Comparing SV burden between primary versus metastatic cancer’).

  2. (2)

    For primary + metastatic LM:

  1. (i)

    Regression P < 0.01

  2. (ii)

    Coefficient P < 0.01

  3. (iii)

    Coefficient more than 0

  1. (3)

    For metastatic-only LM or primary-only LM:

  1. (i)

    Coefficient P < 0.01

  2. (ii)

    Coefficient more than 0

Finally, to determine which features (of those correlated with increased SV-type burden) were enriched in metastatic tumours compared with primary tumours (and vice versa), we calculated Cliff’s delta for numeric features and Cramer’s V for boolean features. Cliff’s delta ranges from −1 to +1, with −1 representing complete enrichment in primary tumours, whereas +1 represents complete enrichment in metastatic tumours. Cramer’s V only ranges from 0 to 1 (with 1 representing enrichment in either primary or metastatic tumours), the sign of the log(odds ratio) was assigned as the sign of the Cramer’s V value so that it ranged from −1 to +1. Features with an effect size of more than 0 were considered as those that could explain the SV burden increase in metastatic cancer when compared with primary cancer.

Driver alterations

We relied on patient-specific cancer driver and fusion catalogues constructed by PURPLE5 and LINX55. Only drivers with a driver likelihood of more than 0.5 were retained. Fusion drivers were filtered for those that were previously reported in the literature. Similarly, we manually curated the list of drivers and removed SMAD3 hotspot mutations because of the high-burden mutations in low-mappability regions. The final driver catalogue contained a total of 453 driver genes and the final fusion catalogue contained 554 reported fusions.

To compare the number of drivers in primary and metastatic tumours, we then combined fusions with the LINX driver variants to calculate a patient-specific number of driver events. Drivers that concern the same driver gene but a different driver type were deemed to be single drivers (for example, TP53 mutation and TP53 deletion in the same sample were considered as one driver event). Cancer-type-specific Mann–Whitney test was performed to assess differences between primary and metastatic tumours. An adjusted q < 0.01 was deemed to be significant.

To assess the driver enrichment, a contingency matrix was constructed from the driver catalogue, containing the frequency of driver mutations per driver type (that is, deletion, amplification or mutations) and cancer type in each cohort (metastatic and primary). A second contingency matrix was constructed for the fusions. Partial amplifications were considered as amplifications, whereas homologous disruptions were considered as deletions. These contingency matrices were filtered for genes that show a minimum frequency of five mutated samples in either the primary or the metastatic cohorts. Then, a two-sided Fisher’s exact test for each gene, cancer type and mutation type was performed and the P value was adjusted for FDR per cancer type. Cramer’s V and the odds ratio were used as effect size measures. An adjusted P < 0.01 was deemed to be significant.

Therapeutic actionability of variants

To determine the amount of actionable variants observed in each sample, we compared our variants annotated by SnpEff (v5.1)56 to those derived from three different databases (OncoKB57, CIViC58 and CGI59) that were classified based on a common clinical evidence level (https://civic.readthedocs.io/en/latest/model/evidence/level.html) as previously described5. In our study we only considered A and B levels of evidence, which represent variants that have been FDA approved for treatment and are currently being evaluated in a late-stage clinical trial, respectively. A variant was determined to be ‘on-label’ when the cancer type matches the cancer type for which the treatment was approved for or is being investigated for, and ‘off-label’ otherwise. Only actionable variants of the sensitive category were considered (that is, tumours containing the variant are sensitive to a certain treatment). Sample-level actionable variants such as TMB high/low or MSI status were not evaluated, because of their tendency to overshadow the other variants, especially in the off-label category. Furthermore, wild-type actionable variants were not considered in this analysis for the same reason. Variants related to gene expression or methylation were not considered due to lack of available data. In addition, we found actionable variants derived from leukaemias to be very different from the solid tumours in our dataset, which is why we excluded them for this analysis. For the analysis of proportion of samples bearing therapeutically actionable variants, we considered that the highest evidence level was retained for each sample following the order A on/off-label to B on/off-label. To assess enrichment of actionable variants globally and at the A on-label level in metastatic tumours, a Fisher’s exact test was performed pan-cancer-wide and per cancer type. An adjusted P < 0.05 was deemed to be significant. Fold changes in frequency are only shown for cancer types with a global significant difference.

To determine which variants contribute the most to the observed significant frequency differences, individual actionable variants were tested for enrichment in metastatic tumours using a Fisher’s exact test per cancer type and tier level. P values were FDR adjusted per cancer type and q < 0.05 was deemed to be significant. In Extended Data Fig. 8, only actionable variants from cancer types with a global significant difference (see above) and that were found at a minimum frequency of 5% in either primary or metastatic cohort and a minimum frequency difference of 5% between them were shown. However, the differences across all screened variants are available as part of Supplementary Table 7.

TEDs

We aimed to pinpoint drivers that are potentially responsible for lack of response to certain cancer treatments in the metastatic cohort. Hence, we devised a test that identifies driver alterations that are enriched in groups of patients treated with a particular treatment type compared with the untreated group of patients from the same cancer type (see Extended Data Fig. 9a for illustration of the workflow).

Treatments were grouped according to their mechanism of action so that multiple drugs with a shared mechanism of action were grouped into the mechanistic treatment category (for example, cisplatin, oxaliplatin and carboplatin were grouped as platinum). We created 323 treatment and cancer-type groups by grouping patients with treatment annotation according to their treatment record before the biopsy. One patient might be involved in multiple groups if they have received multiple lines of therapy or a simultaneous combination of multiple drugs. Only 92 treatment and cancer-type groups with at least ten patients were further considered in the analysis.

Hence, for each cancer type (or subtype, in the case of breast and colorectal) and treatment group, we performed the following steps:

  1. (1)

    We first performed a driver discovery analysis in treatment and cancer-type (or subtype)-specific manner. We explored three types of somatic alterations: coding mutations, non-coding mutations and CN variants (see below for detailed description of each driver category). Driver elements from each alteration category were selected for further analysis.

  2. (2)

    For each driver alteration from (1), we compared the alteration frequency in the treated group to the untreated group of the same cancer type. Each driver category (coding and non-coding mutations and CN variants) were evaluated independently. We performed a Fisher’s exact test to assess the significance of the frequency differences. Similarly, we computed the odds ratio of the mutation frequencies for each driver alteration. The P values were adjusted with a multiple-testing correction using the Benjamini–Hochberg procedure (α = 0.05). An adjusted P value of 0.05 was used for coding mutations and CN variants. An adjusted P value of 0.1 was used for non-coding variants due to the overall low mutation frequency of the elements included in this category, which hampered the identification of significant differences.

  3. (3)

    We then annotated each driver element with information about the exclusivity in the treatment group. We labelled drivers as treatment exclusive if the mutation frequency in the untreated group was lower than 5% or we annotated as treatment enriched otherwise. In addition, we manually curated the identified drivers with literature references of their association with each treatment category.

  4. (4)

    Finally, the overlap of patients in multiple treatment groups (see above) in the same cancer type prompted us to prioritize the most significant treatment association for each driver gene in a particular cancer type. In other words, for each driver gene that was deemed as significantly associated with multiple treatment groups in the same cancer type, we selected the most significant treatment association, unless a driver-treatment annotation was clearly reported in the literature.

The full catalogue of TEDs and their mutation frequencies can be found in Supplementary Table 8.

Coding mutation drivers

We used dNdScv (v0.0.1)60 with default parameters to identify cancer driver genes from coding mutations. A global q < 0.1 was used as a threshold for significance. Mutation frequencies for each driver gene were extracted from the dNdScv output. We defined the mutation frequency as the number of samples bearing non-synonymous mutations.

Non-coding mutation drivers

We used ActiveDriverWGS61 (v1.1.2, default parameters) to identify non-coding driver elements in five regulatory regions of the genome including 3′ untranslated regions (UTRs), 5′ UTRs, long non-coding RNAs, proximal promoters and splice sites. For each element category, we extracted the genomic coordinates from Ensembl v101. Each regulatory region was independently tested. To select for significant hits, we filtered on adjusted P values (FDR < 0.1) and a minimum of three mutated samples. We defined the mutation frequency as the number of mutated samples for each significantly mutated element in the treatment group.

CN variant drivers

We ran GISTIC2 (ref. 62) (v2.0.23) on each of the 92 treatment and cancer-type groups using the following settings:

gistic2 -b <inputPath> -seg <inputSegmentation> -refgene hg19.UCSC.add_miR.140312.refgene.mat -genegistic 1 -gcm extreme -maxseg 4000 -broad 1 -brlen 0.98 -conf 0.95 -rx 0 -cap 3 -saveseg 0 -armpeel 1 -smallmem 0 -res 0.01 -ta 0.1 -td 0.1 -savedata 0 -savegene 1 -qvt 0.1.

The focal GISTIC peaks (q ≤ 0.1 and <1 Mb) were then annotated with functional elements using the coordinates from Ensembl v101. The frequency differences between treated and untreated cohorts on every gene was assessed with Fisher’s exact test as described above. For this, we first calculated the focal amplification and deep depletion status of every gene within each sample. A gene was amplified when the ploidy level of the gene was 2.5 ploidy levels higher than its genome-wide mean ploidy level (as measured by PURPLE), and deleted when the gene ploidy level was lower than 0.3 (that is, deep deletion). We observed that the majority of the peaks contained multiple significant gene candidates (after multiple correction q < 0.05) and therefore we retained the gene most closely positioned to the peak summit, which is the most significantly enriched region across the treated samples. Next, we also found recurrent peaks across multiple treatment groups per cancer type that are not, or less, present in the untreated control group because most of the Hartwig samples have received multiple treatment types. We therefore merged peaks with overlapping ranges to produce a single peak per genomic region per cancer type. For each collapsed peak, we selected the treatment type showing the lowest q value for the gene near the peak summit. Deletion and amplification peaks were processed separately.

Group-level aggregation of treatment resistance-associated variants

To estimate the contribution of TEDs to the total number of drivers per sample in the metastatic cohort, we excluded any TED from the catalogue of driver mutations (see the above section ‘Driver alterations’) in a cancer-type-specific, gene-specific and driver-type-specific manner.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.



Source link