The genomic landscape of 2,023 colorectal cancers

Sample collection

The following steps were taken for sample collection. (1) Ethics approval was provided to the 100kGP by the HRA Committee East of England–Cambridge South research ethics committee (REC reference 14/EE/1112). Samples were obtained as part of the 100kGP cancer programme, an initiative for high-throughput tumour sequencing for NHS patients with cancer^49,50. (2) Thirteen Genomic Medicine Centres (GMCs) were established by the NHS and 100kGP, each with multiple affiliated hospitals across in the same region of the UK. (3) Patients undergoing resection for CRC were identified by specialist nurses and other staff. (4) All patients provided written informed consent, and blood samples were taken. (5) Tumour samples were assessed in histopathology cut-ups. Associated clinicopathological data were obtained from health records. (7) Frozen tumour sub-samples were taken and frozen. Haematoxylin and eosin sections were assessed for purity and other histological features of note. (8) Blood and tumour samples that passed quality control were sent for DNA extraction in regional genetics laboratories. (9) DNA was transferred to the 100kGP central national biorepositry. (10) WGS of paired tumour-constitutional (whole blood-derived) DNA was performed by Illumina. (11) Processed BAM files were transferred to Genomics England for additional processing, quality checking and data storage. (12) All sequencing and clinicopathological data were transferred to Colorectal Cancer Domain (GECIP) for further quality control and data analysis.

WGS and SV calling

Sequencing, mapping and variant calling were generally performed as previously described⁵¹, although we used a less stringent variant allele frequency (VAF) to enable analyses of subclonal mutations.

Sequencing and alignment

Samples were prepared using an Illumina TruSeq DNA PCR-free library preparation kit and sequenced on a HiSeq X, generating 150 bp paired-end reads. Tumour and constitutional DNAs were sequenced to average depths of 100× and 33×, respectively. Poor sequencing quality outliers were identified using principal component analysis and removed on the basis of the following quality metrics: percentage of mapped reads; percentage of chimeric DNA fragments; average insert size; AT/CG dropout; and unevenness of local coverage. Illumina’s North Star pipeline (v.2.6.53.23) was used for the primary WGS analysis. Sequence reads were aligned to the Homo sapiens GRCh38Decoy assembly using Isaac (v.03.16.02.19)⁵². Overall, PCR-free tumour and germline sequencing data for 2,492 fresh-frozen CRC samples were obtained from the 100kGP main program (v.8) release and used in our analysis.

Single-nucleotide variant and indel calling

Single-nucleotide variant and small indel calling was performed using Strelka (v2.4.7). In addition to the default Strelka filters, we applied the following exclusion filters:

• Variants with a germline allele frequency > 1% in the full Genomics England dataset.

• Variants with a population germline allele frequency > 1% in the gnomAD database⁵³.

• Somatic variants with frequency > 5% in the Genomics England cancer dataset. A 5% cut-off was chosen based on the frequency of recurrent non-synonymous variants in Cancer Gene Census genes⁵⁴.

• Variants overlapping simple repeats as defined by Tandem Repeats Finder⁵⁵.

• Indels in regions with high levels of sequencing noise where >10% of the base calls in a window extending 50 bp either side of the indel were filtered out by Strelka owing to the poor quality.

• Indels within 10 bp of 100kGP or gnomAD (v.3) germline indel with allele frequency > 1%.

• Variants in regions of poor mappability where the majority of overlapping 150 bp reads do not map uniquely to the variant position.

• SNVs resulting from systematic mapping and calling artefacts present in both tumour and control 100kGP sample sets. We tested whether the ratio of tumour allele depths at each somatic SNV site was significantly different to the ratio of allele depths at this site in a panel of control samples using Fisher’s exact test. The panel of control was composed of a cohort of 7,000 non-tumour genomes from the Genomics England dataset. At each genomic site, only individuals not carrying the relevant alternative allele were included in the count of allele depths. The mpileup function in bcftools (v.1.9) was used to count allele depths in the PoN. To replicate Strelka filters, duplicate reads were removed and quality thresholds set at mapping quality ≥ 5 and base quality ≥ 5. All somatic SNVs with a Fisher’s exact test phred score < 80 were filtered, with the threshold determined by optimizing precision and recall calculated from a TRACERx truth set⁵⁶.

Removing alignment bias introduced by soft clipping of semi-aligned reads

The Isaac –clip-semialigned parameter invokes the soft clipping of read ends until five consecutive bases are matched with the reference genome. This soft clipping therefore results in the loss of support for alternative alleles occurring within 5 bp of each read end, which leads to artefactually low VAFs. To address allelic bias introduced by this clipping, we introduced FixVAF to soft clip all reads by 5 bp at each end, regardless of whether any of the bases are variant sites or whether the reads support reference or alternate alleles⁵⁷. Reads containing small indels at variant positions were ignored (Supplementary Fig. 1).

Identifying MSI

Tumours with MSI were identified using MSINGS⁵⁸ following the previously described procedure for background model generation (https://github.com/sheenamt/msings/blob/master/Recommendations_for_custom_assays). A set of 132 tumours with known MSI status (106 MSS, 26 MSI) was randomized into test and training sets of 53 MSS and 13 MSI cases (that is, 2 sets of 66 cases). Microsatellite sites were generated using MISA⁵⁹. Only sites overlapping regions of good mappability were considered. Sites measured as unstable in >5 MSS test tumours and sites not unstable in >1 test MSI tumours were removed. The background model produced using the training set was able to perfectly distinguish between MSI and MSS samples in the test set using default MSINGs settings and was then applied to the full CRC cohort.

Identifying pathogenic POL variants

Tumours with pathogenic somatic or germline variants in POLE or POLD1 were identified considering the 22 known pathogenic variants a previously reported⁶⁰. In total, 18 tumours (17 MSS, 1 MSI) had a pathogenic germline (n = 1) or somatic (n = 17) POLE variant and these were considered as a separate POL group in all subsequent analyses. All of the highest mutational burden tumours were either MSI or had a known pathogenic POLE variant, which indicated that no pathogenic polymerase proofreading domain mutations were missed. Tumours with pathogenic POLE variants also exhibited high SBS10a and SBS10b activity, which are established indicators of POLE exonuclease domain mutations¹¹.

CNA calling

Somatic CNAs were called using a framework implemented in the R package CleanCNA (Supplementary Fig. 2). Genome-wide subclonal CNAs were first called using Battenberg (v.2.2.8)⁶¹. To check the quality of these CNA calls, we applied DPClust⁶¹ and CNAqc⁶² to the CNA profiles and SNV VAFs. DPClust clusters variants by their cancer cell fraction (CCF), whereas CNAqc compares observed and expected peaks in SNV VAF distributions to assess CNA calling accuracy. A sample was classified as ‘pass’ if it met both of the following criteria, and ‘fail’ otherwise as follows:

1. 1.

   A clonal cluster of SNVs (0.95 ≤ CCF ≤ 1.05) was identified by DPClust. This clonal cluster was required to have either the highest CCF of all SNV clusters or contain the largest number of SNVs. SNV clusters containing <1% of all sample SNVs were removed before assessment.

2. 2.

   The difference in purity estimates from Battenberg and CNAqc was <5%. CNAqc estimates sample purity considering peaks in SNV VAF distributions in genome regions with one of five copy number states (1:0, 1:1, 2:0, 2:1, 2:2).

CNAs were profiled a maximum of four times per sample and the procedure was stopped if both criteria were met. After a failure, CNA were re-called using Battenberg with re-estimated sample purity and tumour ploidy. After the first fail, purity and ploidy were re-estimated using information from DPClust, where CCF_top is the CCF of the SNV cluster with the greatest CCF:

After the second fail, purity and ploidy were re-estimated using Ccube⁶³, and after the third and fourth fails, purity and ploidy were re-estimated using CNAqc. If a sample failed after four re-runs, then it was removed from downstream analyses reliant on CNAs. Pass CNA profiles were produced for 1,765 out of 2,023 samples.

SV calling

SVs (also referred to as chromosomal rearrangements) represent two reference positions (referred to as rearrangement breakpoints) that are non-adjacent in the reference genome and juxtaposed in a specific orientation. We identified somatic rearrangements using a graph-based consensus approach comprising Delly⁶⁴, Lumpy⁶⁵ and Manta⁶⁶ while also considering support from CNAs (Supplementary Fig. 3). Rearrangements were first called using the three individual callers with default parameters. Delly was run with post-filtering of somatic SVs using all normal samples, as described in the Delly documentation. Rearrangements from the three individual callers were further filtered if any reads supporting the variant were identified in the matched normal, if <2% of tumour reads supported the variant or if either variant breakpoint was in a telomeric or centromeric region or on a non-standard reference contig (that is, not chromosomes 1–22, X or Y). Remaining rearrangements were merged with a modified version of PCAWG Merge SV, which uses a graph-based approach to identify and merge rearrangements identified by multiple callers, allowing a maximum 400 bp difference in breakpoint position to account for variant calling ambiguity¹⁶. Rearrangements were included in the final dataset if they were identified by at least two callers, or by a single caller but with a breakpoint within 3 kb of a CNA segment boundary. SVs were only called in the 1,765 out of 2,023 samples with CNA profiles passing quality control criteria.

Retrotransposition events are mechanistically distinct from other SV-generating events. We searched for retrotransposition events using xTea for LINE-1 elements^67,68,69, as other retrotransposition categories (Alu elements, SINE-VNTR-Alu elements and processed pseudogenes, among others) collectively constitute ≤3% of retrotransposition events across human cancers⁶⁶. We subsequently decided to exclude retrotranspositions from our current SV analysis report, to await later separate publication.

Putative kinase gene fusions were identified considering the following genes: ALK, BRAF, EGFR, ERBB2, ERBB4, FGFR1, FGFR2, FGFR3, KIT, MET, NTRK1, NTRK2, NTRK3, ROS1 and RET²². Fusions were required to involve the kinase domain of the 3′ gene and to have correct strand orientation.

Clinical data

Clinical data were obtained from the GMCs, NHS Digital (NHSD) and Public Health England’s National Cancer Registration and Analysis Service (PHE-NCRAS) through the Genomics England Research Environment as part the 100kGP main program v.10 release. Survival data were obtained from the 100kGP main program v.13 release. Tumour samples sequenced by Genomics England were matched to their respective PHE-NCRAS records using the date of tumour sampling reported by Genomics England and dates of biopsy or treatment reported by PHE-NCRAS, allowing a maximum discrepancy of 7 days.

Clinical data included sex, age at tumour sampling, date of cancer diagnosis, date of last reported follow-up and date of death, tumour histology, tumour type (primary, recurrence of primary or metastases), anatomical site sampled, anatomical site of primary tumour, Dukes stage, and tumour grade (differentiation). For some variables, data were obtained from multiple sources (GMC, NHSD, PHE-NCRAS), and any conflicts between these sources were resolved by individual inspection. If Dukes staging was not available, it was inferred from TNM staging if reported. Anatomical site of primary tumour was reported at different resolutions by the different data sources (for example, one source may report site as proximal colon, whereas another may report it as caecum). To resolve and standardize the site, we therefore constructed an anatomical ontology based on ICD-10-CM codes and assigned sample terms to this ontology. This enabled us to consider anatomical site at two main levels of resolution: less specific (proximal colon, distal colon and rectum) and more specific (caecum, ascending colon, hepatic flexure, transverse colon, splenic flexure, descending colon, sigmoid colon, rectosigmoid colon and rectum). Certain analyses were also performed on the basis of a combined analysis of proximal and distal colon (colon). The proximal colon comprised the caecum, ascending colon, hepatic flexure and transverse colon, whereas the distal colon comprised the splenic flexure, descending colon and sigmoid colon. The rectosigmoid junction was considered part of the rectum. All associations between clinical and molecular data, and between different molecular data, are reported based on tests unless otherwise stated.

Germline mutations in the Mendelian CRC predisposition genes (APC, MSH2, MLH1, MSH6, MUTYH, SMAD4, BMPR1A, GREM1, STK11, NTHL1, MBD4, POLE and POLD1) were explored in the sequenced constitutional DNA. Disease-causing changes were identified based on ClinVar annotation as ‘pathogenic’ or ‘likely pathogenic’, with the exception of POLE and POLD1, which used the method described in the section ‘Identifying pathogenic POL variants’. Evidence of pathogenic biallelic changes was required to diagnose the recessive conditions (MUTYH, NTHL1 and MBD4) and no such cases were found. Twenty patients (aged 30–79 years) were identified as having a previously unreported CRC predisposition caused by germline mutations in Lynch syndrome or polymerase proofreading polyposis genes (seven MSH2, five MLH1, six MSH6, one POLE, one POLD1).

Based on principal component analysis of germline genotypes, 90.2% (n = 1,819) patients were of European ancestry, with 2.6% (n = 52) African, 0.7 (n = 15) East Asian, 3.2% (n = 64) South Asian and 3.3% (n = 67) mixed ancestry (Supplementary Fig. 4). There was strong agreement between 16 self-reported ancestry groups and principal component analysis classification.

Sample selection

Because tumour sample purity and sequencing data quality affect the sensitivity and precision of variant calling⁷⁰, we excluded samples using the following quality control procedures (Supplementary Table 2).

• Tumour samples were excluded if cross-contamination of the tumour sample was >1%, as estimated by VerifyBamID⁷¹.

• Tumour samples were excluded if cross-contamination of the matched germline sample was >1%, as estimated by VerifyBamID.

• Estimating tumour sample purity is particularly difficult when purity is low. We therefore used the distribution of single-nucleotide variant VAFs to identify low purity samples, as a low average SNV VAF can be indicative of low sample purity⁷². Tumour samples with a median SNV VAF < 0.1 were excluded, with this threshold chosen based on the smaller numbers of potential driver variants observed in MSS CRC samples when compared with all MSS CRC samples (Supplementary Fig. 5). Here driver mutations were defined as any potentially pathogenic coding variant called in 63 driver genes previously identified in MSS CRC^3,4,7,8.

• Tumour samples were excluded if <100 SNVs were called, as this number is below the smallest number of SNVs previously reported in CRC whole genomes^{2,3,4,5,6,7,8,9} and therefore suggestive of low sample purity or sequencing data quality.

• Tumour samples were excluded if many mutations were associated with a probable artefactual mutational signature¹⁵.

In total 286 out of 2,492 (11.5%) tumour samples were excluded based on the above criteria.

Tumour samples were also excluded if essential clinical data were missing or there were unresolvable conflicts between the sources from which clinical data were obtained (GMCs, NHSD, PHE-NCRAS) (Supplementary Table 2). In total, 183 out of 2,206 (8.3%) of tumour samples that passed tumour sample purity and sequencing data quality control were excluded based on clinical data, using the following criteria:

• GMC, NHSD and PHE-NCRAS reported conflicting years of birth.

• Sex reported by GMC, NHSD and/or PHE-NCRAS did not match the sex inferred from sequencing data.

• GMC, NHSD and PHE-NCRAS did not report tumour histology or reported conflicting histology.

• Tumour was not classified as a colorectal adenocarcinoma.

• Missing or conflicting data meant it was unclear whether the primary tumour or a metastasis was sampled.

• If multiple primary tumours or multiple metastases from a single individual were sequenced, the primary tumour or metastasis sample with the highest purity was included, and all other primary tumour or metastasis samples were excluded. This procedure was completed after all other exclusion criteria had been applied. Primary tumours and metastases were considered separately for this procedure.

Based on these criteria, 2,023 colorectal adenocarcinoma samples were suitable for analysis (Supplementary Table 2). This cohort comprised 1,898 primary tumours, 122 metastases and 3 recurrences of primary tumours from 2,017 patients. Six patients (all MSS) had both a primary tumour and a metastasis sample sequenced and each tumour was included. One hundred and nineteen metastases were MSS, the other three comprising two MSI and one POL cancer. Some subsequent analyses excluded the MSI and POL metastases (details in Supplementary Tables). The three recurrences were MSS (n = 1) and MSI (n = 2), and these were included in the appropriate primary cancer group for further analyses. A single cancer was POL and MSI, and this was included in the POL group for further analyses. Clinical data completeness is detailed in Supplementary Table 31.

Single-nucleotide variant and indel drivers

Mutation annotation

Somatic mutations were annotated to Ensembl (v.101, GRCh38) using Variant Effect Predictor (VEP)⁷³. The following parameters were used: vep -i –assembly GRCh38 –no_stats –cache –offline –symbol –protein -o

–vcf –canonical –dir –hgvs –hgvsg –fasta –plugin CADD, –plugin UTRannotator,.

The CADD score file was obtained using CADD (v.1.6)^74,75,76, with scores attained for all SNV and indel mutations using the CADD software available from GitHub (https://github.com/kircherlab/CADD-scripts) before being utilized by the VEP CADD plugin.

Protein-coding driver identification

Protein-coding driver genes were identified using the IntOGen pipeline (v.2020, downloaded February 2021)¹⁸. Identification was performed separately in MSS primary, MSI (all primary), POL (all primary) and MSS metastasis sample sets, with the aim of optimizing correction for varying background mutation rates and spectra among these four groups. Subsequent analyses restricted discovery to specific anatomical locations or cluster groups in MSS primary tumours.

Pre-processing of input mutations

Somatic mutations passing the filtering criteria described above were subject to initial sample and mutation pre-processing. In the case of multiple tumours from the same patient, the primary tumour was used. Within each cohort (that is, MSS primary, primary MSI, primary POL, MSS metastasis), tumours were flagged for exclusion from downstream driver gene identification if they contained >10,000 mutations and had an outlier mutation count, defined as upper quartile + (1.5 × interquartile range). Mutations present in a Hartwig Consortium panel of control set were also excluded⁷⁷. Unless otherwise specified, mutations were mapped to canonical protein-coding transcripts from Ensembl (v.101, GRCh38).

Driver identification methods

Seven driver gene identification methods were run through the IntOGen pipeline (Supplementary Fig. 6):

1. 1.

   dNdSCV (v.0.1.0)⁶ is designed to detect genes under positive selection that show an excess of non-synonymous (missense, nonsense, essential splice) mutations after correction for local trinucleotide context. In the primary POL cohort the parameter ‘max_coding_muts_per_sample = Inf’ was used because of the high proportion of hypermutated tumours.

2. 2.

   OncodriveFML (v.2.4.0)⁷⁸ aims to detect driver genes that show an enrichment of mutations with high functional impact. CADD scores were used as measure of functional impact^74,75,76.

3. 3.

   OncodriveCLUSTL (v.1.1.3)⁷⁹ is a method designed to detect driver genes that are enriched for linear mutation clusters. In the primary POL cohort, pentamer signatures were used rather than trinucleotide signatures because of the improved performance of the pentanucleotide-based background models compared with that of trinucleotides in these tumours.

4. 4.

   cBaSE (v.1.1.3)^18,80 aims to detect driver genes under positive selection that exhibit a significant mutation count bias after correction by trinucleotide context.

5. 5.

   MutPanning (v.2)⁸¹ is designed to detect driver genes that exhibit enrichment of mutations with unusual nucleotide contexts compared with a background model.

6. 6.

   HotMaps3D (v.1.1.3)^18,82 detects driver genes containing missense mutations that are spatially clustered together in the three-dimensional structure of the protein. Protein structures were downloaded from The Protein Data Bank⁸³ in March 2020.

7. 7.

   smRegions (v.1)⁸⁴ detects genes containing an enrichment of non-synonymous mutations in regions of interest, such as protein domains, after correcting for trinucleotide context. This analysis utilized information from protein family (Pfam) domains that were mapped to Ensembl (v.101) canonical transcripts.

Combination of driver identification methods

The results of the seven driver identification methods were combined in similar manner as previously described¹⁸. In brief, the driver combination procedure considered the top 100 ranked genes and their associated P and Q values in each of the seven driver identification methods. Somatically mutated genes assigned as tier 1 or tier 2 in the COSMIC Cancer Gene Census (CGC; v.92)⁵⁴ were designated as the truth set of known drivers. Through comparison of the relative enrichment of CGC genes in the top ranked gene lists, a per-method weighting was obtained. Per-method ranked lists were combined using Schulze’s voting method to generate a consensus ranking, with combined P values estimated using a weighted Stouffer Z score method.

Driver candidates were then classified into the following tiers:

• Tier 1: candidates for which the consensus ranking was higher than the ranking of the first gene with Stouffer Q ≤ 0.05. These represent high-confidence drivers.

• Tier 2: candidates not meeting the criteria for tier 1, but which are CGC genes and showed a combined Stouffer Q_CGC < 0.25. These represent a set of ‘rescued’ known cancer drivers.

• Tier 3: candidates not meeting the criteria for tier 1 or tier 2 but with Stouffer Q < 0.05. These represent lower confidence drivers.

• Tier 4: candidates not meeting criteria for tier 1 or tier 2 and with Stouffer Q > 0.05. These genes are not likely to be drivers.

Post-processing of candidate drivers

Candidate driver genes were filtered based on the following annotations:

1. 1.

   Automatic fail: a candidate driver gene would be excluded from further consideration if annotated with at least one of the following:

   1. a.

      Tier 4: categorized as tier 4 by the combination procedure.

   2. b.

      Single method: only significant (Q < 0.1) in one of the seven methods (non-CGC genes).

   3. c.

      Expression: gene has very low or no expression in a relevant tumour type based on data from The Cancer Genome Atlas (TCGA).

   4. d.

      Olfactory receptor: gene is in list of olfactory receptor genes.

   5. e.

      Known artefact: gene is in a list of known artefacts or long genes (for example, TTN).

2. 2.

   Manual review: if a gene is not excluded based on any automatic fail filters, it is retained as a candidate driver:

   1. a.

      Germline: non-tier 1-CGC gene has ≥1 mutations per sample and oe_syn/ms/lof > 1.5 based on gnomAD (v.2.1) constraint metric estimates.

   2. b.

      Sample 3 Muts: non-CGC gene for which there are ≥3 mutations in ≥1 tumour.

   3. c.

      Literature: non-CGC gene for which there are no literature annotations according to CancerMine⁸⁵.

3. 3.

   Automatic pass: is not flagged by any automatic fail or manual review filters.

Candidate driver roles were assigned on the basis of dN/dS ratios for missense (wmis) and nonsense (wnon) mutations for the given gene derived from dNdSCV (https://bitbucket.org/intogen/intogen-plus/src/master/core/intogen_core/postprocess/drivers/role.py):

• A distance metric was calculated by distance = ((wmis – wnon))/√2

• Candidate drivers with distance >0.1 represent those with an excess of missense to nonsense mutations and are therefore considered oncogenes.

• Candidate drivers with distance <0.1 represent those with an excess of nonsense to missense mutations and are therefore considered TSGs.

• Otherwise, the role of the candidate driver is unclear and considered ambiguous.

In the case of multiple cohorts being run representing subsets of a given tumour type, a consensus role was designated comparing between each subtype role:

Gene candidates were annotated by their overlap with any IntOGen cohorts from a previous IntOGen pan-cancer analysis (1 February 2020) as well as from a pan-cancer TCGA analysis².

Noncoding driver identification

Defining sets of noncoding regions

Regions from candidate noncoding elements overlapping coding sequence (CDS) or exon regions from canonical protein-coding transcripts were removed using bedops (v.2.4.39)⁸⁶.

The following sets of noncoding regions were defined:

1. 1.

   Core promoters (n = 19,283). Defined based on the transcription start site (TSS) of canonical protein-coding transcripts: 200 bp < TSS < 50 bp. CDS regions were removed.

2. 2.

   Distal promoters (n = 19,296). Defined based on the TSS of canonical protein-coding transcripts: 2 kb < TSS. CDS regions were removed.

3. 3.

   5′ untranslated regions (UTRs; n = 18,613). Defined based on canonical protein-coding transcripts. CDS regions were removed.

4. 4.

   3′ UTRs (n = 18,806). Defined based on canonical protein-coding transcripts. CDS regions were removed.

5. 5.

   lincRNAs (n = 16,510). Based on exon regions from transcripts annotated as lincRNAs in Ensembl (v.101). Exon regions from canonical protein-coding transcripts were removed.

6. 6.

   miRNAs (n = 1,793). Based on regions from transcripts annotated as miRNAs in Ensembl (v.101). Exon regions from canonical protein-coding transcripts were removed.

7. 7.

   Non-canonical splice regions (n = 18,163). Defined from regions extending 30 bp into the intron from essential splice donor or acceptor sites in canonical protein-coding transcripts. Exon regions from canonical protein-coding transcripts were removed.

8. 8.

   Enhancers (n = 130,996). Defined from Ensembl (v.101) regulatory elements annotated as ‘enhancer’. Exon regions from canonical protein-coding transcripts were removed.

9. 9.

   Open chromatin regions (n = 95,344). Defined from Ensembl (v.101) regulatory elements annotated as ‘open chromatin’. Exon regions from canonical protein-coding transcripts were removed.

10. 10.

    CTCF sites (n = 173,711). Defined from Ensembl (v.101) regulatory elements annotated as ‘CTCF sites’. Exon regions from canonical protein-coding transcripts were removed.

11. 11.

    Transcription factor-binding sites (n = 29,259). Defined from Ensembl (v.101) regulatory elements annotated as ‘TF binding sites’. Exon regions from canonical protein-coding transcripts were removed.

Detecting noncoding drivers

Potential noncoding driver mutations were identified in non-hypermutated MSS primary tumours (n = 1,442). OncodriveFML (v.2.4.0) was run on sets of noncoding regions according to the following amended parameters from the protein-coding analysis: indel-max indels are treated as a set of substitutions, with the functional impact of the indel mutation being the maximum of all the substitutions, and the background simulated as substitutions. A Q < 0.01 threshold was considered as significant (Supplementary Fig. 7).

SNV mutations exhibiting extreme strand bias

SNV mutations that otherwise passed filtering criteria as previously detailed were further scrutinized for excessive strand bias (Strelka INFO field SNVSB > 10). This highlighted many missense mutations that cause a recurrent missense change in CACNA1E (p.Ile95Leu); these exhibited excessive strand bias and were therefore deemed false calls.

Driver mutation annotation

Non-synonymous mutations in the 682 gene transcripts considered by OncoKB (v.3.3) were annotated using the OncoKB API⁸⁷. In the first instance, the HGVSg identifier was used; in the rare instances that this failed, a combination of gene symbol, consequence and HGVSp were used to map mutations to OncoKB annotations.

Annotation of oncogenic mutations

Non-synonymous mutations in candidate driver genes were annotated as pathogenic if any of the following criteria were met:

1. 1.

   The mutation is annotated by OncoKB as ‘oncogenic, ‘likely oncogenic’ or ‘predicted oncogenic’.

2. 2.

   The driver is classified as an oncogene, the mutation consequence is missense, and the mutation is recurrent (seen in ≥3 tumours in cohort).

3. 3.

   The driver is classified as a TSG or ambiguous and either:

   1. a.

      Consequence is protein-truncating (splice acceptor, splice donor, frameshift, stop lost, stop gained or start lost).

   2. b.

      Consequence is missense and mutation is recurrent (seen in ≥3 tumours in cohort).

For POLE, oncogenic annotations were restricted to missense mutations in the exonuclease domain (amino acid residues 268–471).

Non-synonymous mutations not meeting these criteria were considered as variants of uncertain significance.

Lollipop plots of driver gene mutations

Lollipop plots of driver gene mutations (Supplementary Result 2) were generated using the Rpackage trackViewer⁷⁹. Pfam protein domains mapping to the Ensembl (v.101) canonical transcripts were plotted. The protein position was taken from the first position in the HGVSp annotation, apart from splice donor and acceptor mutations, for which the codon nearest to the HGVSc transcript position was assigned as the protein position.

Timing driver mutations

The relative evolutionary timings of candidate driver mutations were obtained using MutationTimeR³¹. Copy number input for MutationTimeR was prepared from Battenberg segmentation files, with the clonal frequency of each segment taken as the tumour purity. In the case of subclonal calls, the clonal frequency was calculated by multiplying the tumour purity by the clonal fraction. The clusters input for MutationTimeR was prepared from DPClust cluster estimates. The VAF proportion was calculated by multiplying the estimated cluster CCF by the tumour purity. Superclonal clusters (CCF > 1.1) were removed. VCF input for MutationTimeR was obtained from the small somatic SNV/indel variant VCFs, which had been filtered as previously described. For SNVs, alt and ref depths were obtained using FixVAF. For indels, ref and alt depths were obtained from tier 2 Strelka TAR and TIR fields, respectively. Only mutations within Battenberg copy-number segments were retained (note that for male XY tumours with only 1 copy of the X chromosome, copy number information is restricted to the pseudoautosomal region and Battenberg was not run on the Y chromosome).

MutationTimeR was run with 1,000 bootstraps. For tumours previously defined as having undergone WGD, the parameter isWgd was set to true. Mutations were then classified into estimated simple clonal states (as per figure 1a of ref. ³¹): clonal (early), mutation on ≥ 2 copies per cell; clonal (late), mutation on 1 copy per cell, no retained allele; clonal (NA), mutation on 1 copy per cell, either on amplified or retained allele; subclonal, mutation on <1 copy per cell.

Mutational signature attribution

SeqInfo VCFs produced as part of SigProfilerMatrixGenerator¹⁷ were used to map somatic mutations from input VCFs to their SBS96, DBS78 or ID83 contexts and then to the final SigProfilerExtractor COSMIC (v.3.2) decomposed signature probabilities. For different purposes, mutational signatures were variously measured as follows: presence–absence, for example, when assessing shared aetiology; proportional activity (essentially proportion of mutations fitted to any signature in that tumour), useful for comparing between signatures in the same sample; and number of mutations ascribed, estimated as (activity × burden of mutations of SBS, DBS or ID type fitted to any signature), approximating to burden of mutations from that signature in that tumour.

Annotation of DBS mutations

Per-tumour VCFs containing DBS mutations, either directly called originally by Strelka or originally called by Strelka as two adjacent SNVs and reconstructed as DBS mutations, were created and mutation consequences were re-calculated using VEP as above.

Patterns of somatic CNA

WGD classification

Tumours were classified as WGD considering the average genome copy number state (ψ_ave) as follows:

Where S is the number of copy number genome segments, \({C}_{{i}_{{\rm{Maj}}}}\) and \({C}_{{i}_{{\rm{Min}}}}\) are the major and minor allele copy numbers, respectively, for genome segment i, and L_i is the base pair length of genome segment i. If there was evidence of subclonal alteration, then the copy number states corresponding to the largest tumour cell fraction were considered. Tumours were classified as WGD if 2.9–2H < ψ_ave and non-WGD otherwise, where H is the fraction of the genome with a minor allele copy number of 0 (ref. ³²).

Classification of CNAs

Individual CNAs were grouped into six categories: homozygous deletion (HD), LOH, including copy-neutral LOH, other loss (OLOSS), no change (NOC), gain (Gain) and amplification (AMP). The classification considers whether a tumour has undergone WGD (Supplementary Table 38).

For cases in which subclonal CNAs existed, the copy number state corresponding to the largest cell fraction was used. Classification into one of the six categories overlaps significantly between non-WGD and WGD tumours, with differences relating to total copy number. Differences include the following:

• In non-WGD tumours, segments were classified as LOH if 1 allele had a copy number state of 0 and the total copy number (t_CN) ≤2. In WGD tumours, segments were classified as LOH if 1 allele had a copy number state of 0 and t_CN ≤4.

• Non-WGD tumours do not have an OLOSS category.

• NOC was defined as 1+1 in non-WGD tumours and 2+2 in WGD tumours.

• In non-WGD tumours, segments were classified as Gain if 2 < t_CN ≤ 5. In WGD tumours, segments were classified as Gain if 4 < t_CN ≤10.

• In non-WGD tumours, segments were classified as AMP if t_CN > 5. In WGD tumours, segments were classified as AMP if t_CN > 10.

Positional enrichment of CNAs

Preparing GISTIC input

Recurrent arm-level copy number events, as well as focal amplifications and deletions, were identified using GISTIC (v2.0.2.3)³⁴. For all samples with CNA profiles passing quality criteria, a copy number segmentation file suitable for GISTIC input was generated using Battenberg output. Chromosomal coordinates and major (n_Maj) and minor (n_Min) copy number states were obtained for each copy number segment identified by Battenberg. In the case of subclonal copy number segments, n_Maj and n_Min values corresponding to the largest tumour cell fraction were considered.

Per-segment normalized copy number (SegCN) values were calculated differently for tumours with WGD (for which ploidy was assumed to be four) and without WGD (for which ploidy was assumed to be two). SegCN was thresholded to a minimum of –2 and maximum of 2.

For non-WGD tumours, SegCN was calculated as follows:

For non-WGD tumours from males, X chromosome SegCN was calculated as follows:

For WGD tumours, SegCN was calculated as follows:

For WGD tumours from males, X chromosome SegCN was calculated as follows:

Running GISTIC

GISTIC was run using the following parameters: -conf 0.99 -broad 1 -qvt 0.25 -genegistic 1 -gcm extreme -brlen 0.5 -rx 0 -twoside 1 -scent median -armpeel 1 -arb 1 -refgene hg38.UCSC.add_miR.160920.refgene.mat.

Prioritizing probable gene targets of focal amplifications and deletions

Candidate target genes at focal amplifications and deletions were annotated using the following criteria:

1. 1.

   Overlap with genes at focal amplifications and deletions reported in a previous pan-cancer study that used GISTIC⁸⁸. Comparisons were made both with the overall pan-cancer GISTIC analysis, and GISTIC analysis was restricted to the given tumour type. Special consideration was given to genes specifically highlighted by the previous study⁸⁸ as being candidates.

2. 2.

   Overlap with Cosmic Cancer Gene Census genes and whether their annotated role (oncogene (OG), TSG or ambiguous) is consistent with the copy number change (OG with amplifications and TSG with deletions)²².

3. 3.

   Overlap with driver genes identified in this study and whether their probable role (OG, TSG or ambiguous) is consistent with the copy number change (OG with amplifications and TSG with deletions).

Based on the above criteria, consensus driver genes were manually assigned to peaks. Comparisons were made with all potential gene synonyms as available from the HUGO gene nomenclature name committee (https://www.genenames.org/).

Defining copy number segments overlapping recurrent CNAs

Alterations from the broad analysis with Q < 0.05 were taken to indicate recurrent arm-level events. Copy number segments constituting greater than half of the total chromosome arm size were taken to indicate arm-level events.

In the case of focal events identified by GISTIC, the ‘wide region’ was used to compare potential extent of overlap with copy number segments. Segments were defined as overlapping focal events if either the segment interval constituted greater than half of the focal region, or vice versa, using pybedtools and bedtools (v.2.3.0)^89,90.

Tumours were considered to have specific arm-level or focal deletions if an overlapping copy number segment was annotated as HD or LOH (as described above). Similarly, tumours were considered to have specific arm-level or focal amplifications if an overlapping copy number segment was annotated as Gain or AMP. In the case of subclonal CNAs, n_Maj and n_Min values corresponding to the largest cell fractions were considered.

ecDNA detection

With the caveat that that there is no definitive way to distinguish ecDNA and intrachromosomal amplification in heterogeneously staining regions, potential ecDNA molecules were detected from tumour bam files using AmpliconArchitect (v.1.2)²⁹. In brief, per-tumour seed regions were prepared from Battenberg copy number segmentation output if a segment was >100 kb and the total copy number was >5. AmpliconArchitect was then run using these seed regions to extract overlapping sequence reads from the tumour BAM file and to construct candidate amplicons.

Candidate amplicons were classified using AmpliconClassifier (v.0.4.6) into the following categories: (1) cyclic (truly circularized ecDNA); (2) complex non-cyclic; (3) linear amplification; and (4) no amplification or invalid. Amplicons were highlighted if containing a known highly amplified oncogene (MDM2, MYC, EGFR, CDK4, ERBB2, SOX2, TERT, CCND1, E2F3, CCNE1, CDK6, MDM4, NEDD9, MCL1, AKT3, BCL2L1, ZNF217, KRAS, PDGFRA, AKT1, MYCL, NKX2-1, IGF1R and PAX8, as previously reported³⁰).

Estimation of telomere content

Telomere content was estimated from tumour and germline BAM files using TelomereHunter (v.1.1.0)⁹¹ and Telomerecat (v.3.3.0)⁹² with default parameters.

Telomere content was normalized by log₂(tumour content/normal content).

Patterns of somatic structural variation

Classification of simple and complex SVs

Rearrangements identified by the graph-based consensus approach were grouped into footprints and clusters based on their proximity within the genome, the overall number of events in the genome and the size of these events using ClusterSV⁹³. Rearrangement footprints represent sets of rearrangement breakpoints that are positionally associated, whereas rearrangement clusters represent sets of rearrangements that are mechanistically associated. Rearrangement footprints were described using the string approach as previously proposed⁹³. Simple and complex events were defined as clusters comprising ≤2 or ≥3 individual rearrangements, respectively. Simple events were classified as deletions, tandem duplications, balanced inversions, balanced translocations or unbalanced translocations, whereas complex events were classified as chromoplexy or chromothripsis (detailed further below). Simple and complex events that did not meet the criteria of any of these classifications were described as simple unclassified or complex unclassified, respectively.

Chromothripsis events were inferred using established criteria^93,94. A rearrangement cluster was defined as chromothripsis if it met all the following criteria:

• A contiguous series of four genome segments oscillating between two copy number states, or five genome segments oscillating between three copy number states.

• At least six interleaved intrachromosomal rearrangements, as per a previous study⁹⁴.

• No evidence (FDR > 0.2) that the distribution of intrachromosomal fragment join orientations diverge from a multinomial distribution with equal probabilities for each of the four orientation categories (duplication-like, deletion-like, head-to-head inversion and tail-to-tail inversion).

A rearrangement cluster was defined as chromoplexy if it met all the following criteria:

• Contains a chain of rearrangements spanning at least three chromosomes⁹⁵. SV chains were identified using a graph-based approach, in which nodes represent breakpoints, and are connected by an edge if they are not involved in the same rearrangement and fall within 1 Mb of each other. The graph-based approach was implemented using the R package igraph⁹⁶.

• At least 50% of rearrangement footprints in the cluster represent balanced translocations, either with no observed copy number change or a deletion bridge between the break ends.

• Consists of between 3 and 30 rearrangements.

Identification of simple structural variation hotspots

Rates of somatic structural variation differ throughout the genome and are influenced by local genomic features⁹⁷. Genome regions enriched for simple SVs (Supplementary Table 10) were therefore identified using a permutation-based approach considering genomic features associated with structural variation occurrence. Deletions, tandem duplications, balanced inversions, balanced interchromosomal translocations and unclassified simple SVs were considered separately. Individual rearrangements forming parts of complex SVs were excluded from this analysis. MSS primary and MSI tumours were also analysed separately, whereas primary POL tumours and metastases were not considered owing to low sample numbers.

Evaluating relationships between genomic features and SV frequencies

Negative binomial regression was used to test associations between genomic features and numbers of SVs of each simple class⁹⁷. The following features were included in the models: average total copy number across the bin in the CRC sample set, GC content, the presence of genes highly or lowly expressed in CRC, ALU repeats, other genomic repeats, segmental duplications, fragile sites, replication timing, and DNase, H3K36me3 and H3K9me3 peaks. Highly and lowly expressed genes were defined as those with mean RSEM value in the top 25% and bottom 75% of protein-coding genes in TCGA CRC samples with RNA sequencing⁷. ALU and other genomic repeats were obtained from the UCSC Genome Browser⁹⁸. Segmental duplications were obtained for GRCh38 from the Segmental Duplication Database⁹⁹. Fragile sites were obtained from a previous study⁶⁶. Replication timing data from CRC epithelial cells (HCT116) were obtained from ReplicationDomain¹⁰⁰. DNase-seq data (ENCFF443KCU) and ChIP–seq data for histones H3K36me3 (ENCFF553QXG) and H3K9me3 (ENCFF482DLD) were obtained for the large intestine from ENCODE¹⁰¹.

Permuting SVs

SVs were simulated to test whether the number of SVs observed in a region was greater than expected by chance given the local genomic features¹⁰². SVs were simulated for each simple SV class, preserving the number and length (distance between intrachromosomal SV break ends) of SVs observed in the CRC sample sets. To simulate SVs, the genome was divided into non-overlapping 1 Mb bins and the genomic features (listed above) of each bin summarized. All genomic features were normalized to a mean of 0 and standard deviation of 1 to aid comparisons. The number of break ends expected in each bin was then estimated using the effect estimates from the previously generated negative binomial regression model. For each observed SV, a SV was simulated by sampling a bin under probabilities proportional to the expected numbers of break ends in each bin. For intrachromosomal SVs, a partner break end was then simulated by selecting the position either upstream or downstream (with equal probability) equal in distance to the distance between the two break ends in the observed SV. For interchromosomal SVs, a partner break end was simulated by sampling a bin under probabilities proportional to the expected numbers of break ends in each bin, excluding bins on the same chromosome. SVs were re-simulated if either break end fell within an uncallable region (a telomere or centromere). SVs were simulated 1,000 times to generate a null distribution of expected SV numbers for the 1-Mb bins.

Identifying SV hotspots

Piece-wise constant fitting (PCF) was used to identify regions of the genome containing greater numbers of SV break ends than expected¹⁰². SV break ends were first sorted by position and the distance between successive break ends calculated. PCF was then applied to the log₁₀ of these inter-mutational distances (IMDs). SV hotspots were identified by first computing the observed (d^obs_i) and expected (d^exp_i) number of breakends per base pair for each PCF segment (i):

Where a_i is the number of break ends in the segment, s_i is the length of the segment in base pairs, n is the number of bins overlapping the segment, b_j is the expected number of SVs in bin j, and s^bin is the bin size (1 Mb). A simple SV enrichment factor β_i^simple is then computed for each PCF segment as follows:

The PCF algorithm requires parameters γ (that controls the smoothness of the segmentation) and k_min (the minimum number of mutations in a segment). FDRs at each β^simple value were estimated by applying PCF to both the observed and simulated SV sets and dividing the mean number of segments with a β^simple value at least as great in the simulated SV sets by the number of segments with a β^simple value at least as great in the observed SV set. A maximum FDR of one was set and FDR values equal to zero were changed to the lowest non-zero FDR value observed. Optimal γ and k_min values were chosen by repeating this process for values of γ between 1 and 20, and values of k_min between 2 and 20, and selecting values that maximized the number of hotspots identified while minimizing the FDR. In the final analysis, γ = 10 was used throughout, whereas k_min = 2 was used for translocations in MSS primary samples, k_min = 4 was used for unclassified simple variants in primary MSI samples, and k_min = 10 was used otherwise. SV hotspots for which no SVs were supported by CNAs were considered potential artefacts and removed. Overlapping SV hotspots identified in the same sample sets were collapsed.

Classification of SV hotspots as fragile sites

SV hotspots were classified as fragile sites if they satisfied at least three of the following six criteria (this threshold was chosen by assessing the co-occurrence of these criteria):

• Were late replicating¹⁰³. Replication timing data from CRC epithelial cells (HCT116) were obtained from ReplicationDomain. Late-replicating regions were defined as those with mean Repli-Seq values ≤ 0.

• Had low gene density¹⁰⁴. A threshold of five genes per megabase was used.

• Overlapped a gene greater than 300 kb in size. This threshold was chosen as fragile sites generally occur in chromosome regions containing genes at least 300 kb in size¹⁰⁵.

• The overlapping gene of greatest size was the focus of the SV enrichment. This was assessed by computing the ratio between SV break point densities in the overlapping gene of greatest size and intergenic regions flanking 1 Mb upstream and downstream. A threshold of five was used.

• Overlapped a fragile site as previously reported⁶⁶. These fragile sites were originally obtained from either the NCBI or literature curation and were mapped from NCRI36 to GRCh38 co-ordinates using LiftOver⁹⁸.

• Overlapped a fragile site identified in a pan-cancer analysis of whole-genome-sequenced tumours¹⁰² and mapped from GRCh37 to GRCh38 co-ordinates using LiftOver⁹⁸.

SV hotspots were not considered as potential fragile sites if they contained an identified CRC driver gene. SV hotspots at potential fragile sites likely occur for mechanistic rather than selective reasons and were therefore not considered further⁶⁶.

Identification of candidate gene targets of recurrent SVs

Genes were reported as candidate targets of recurrent SVs if they had been identified as targets in previously analyses^4,7,8,102, were known CRC driver genes overlapping an SV hotspot or were the sole expressed gene in the hotspot region. Numbers of samples with a focal change at a candidate gene were computed considering SVs <3 Mb in size¹⁰⁶ at least partially overlapping the gene coding sequence.

Enrichment of complex structural variation

Genome regions enriched for complex SVs were identified using a permutation-based approach, considering chromothripsis, chromoplexy and unclassified complex SVs separately. MSS primary and MSI tumours were also analysed separately, whereas POL tumours and were not considered owing to low sample numbers. The genome was first split into non-overlapping 1 Mb bins and the observed number of tumour samples with complex SV footprints (g^pbs_i) overlapping each bin (j) counted. Complex SV footprint positions were next permuted 100,000 times by randomly sampling genome regions equal in size to the footprints. The expected number of tumour samples with complex SV footprints (g^exp_i) overlapping each 100-kb bin was then estimated as the mean number of tumour samples with SV footprints overlapping the bin across all permutations. A complex SV enrichment factor β_i^complex was calculated for bin (j) as follows:

FDRs at each β^complex value were estimated by computing β^complex for each bin in both the observed and permuted SV sets and dividing the mean number of bins with a β^complex value at least as great in the permuted SV sets by the number of bins with a β^complex value at least as great in the observed SV set. A maximum FDR of 1 was set and FDR values equal to zero were changed to the lowest non-zero FDR value.

Mutational processes

Characterizing SBS, DBS and indel signatures

SBS, DBS and indel signatures were extracted de novo and related to known COSMIC signatures (v.3.2) using SigProfilerExtractor¹¹. SBS, DBS and indel signatures were extracted using random initialization, 500 NMF replicates, and between 10,000 and 1,000,000 NMF iterations. We assumed the presence of between 1 and 30 SBS signatures (minimum signatures and maximum signatures parameters, respectively), 1 and 15 DBS signatures, and 1 and 10 indel signatures. Default settings were used for all other parameters. Investigation of the new DBS-A signature (Supplementary Table 3) hinted towards the signature being a technical artefact of the high number of short indels at homopolymer regions occurring in MSI samples.

Characterizing SV signatures

SV signatures were extracted considering only simple SVs, specifically deletions, tandem duplications, balanced and unbalanced inversions, and balanced and unbalanced interchromosomal translocations. Deletion and tandem duplication size distributions are multimodal, and we therefore classified these variants as <10 kb, 10 kb to 1 Mb, and > 1Mb. Variant site replication timing is also multimodal and we therefore classified variants as late, mid or early replicating considering mean Repli-Seq thresholds of <–2, –2 to 2, and >2 using CRC epithelial cell data from ReplicationDomain (Supplementary Fig. 8).

Mechanisms of fragile site instability differ from other SVs, and deletions and tandem duplications at fragile sites were therefore considered separately⁹¹. Signatures were extracted using a hierarchical Dirichlet process (HDP) implemented in the R package hdp (v.0.1.5)⁹¹. The hierarchical Dirichlet process structure was initialized with one common grandparent node, a parent node for each of the MSS, MSI and POL tumour subtypes, and a child node for each of the 1,765 tumour samples in which SVs were called. Four separate Markov chain Monte Carlo posterior sampling chains were run with 5,000 burn-in iterations, extracting 12 SV signatures. Extraction stability was assessed by splitting the cohort into halves, maintaining proportions of MSS, MSI and POL tumours, and re-extracting signatures from each half. Nine signatures extracted from the cohort halves showed high similarity between halves (cosine similarity > 0.9) and high similarity with signatures extracted from the full cohort. These nine signatures were named SV1–SV9 and considered in subsequent analyses.

To investigate DNA repair mechanism perturbation, we correlated driver gene mutation with SV signature activity. A gene was considered mutated if it harboured a likely pathogenic germline SNP or indel (variants annotated as ‘pathogenic’ or ‘likely pathogenic’ in ClinVar¹⁰⁷), a likely oncogenic somatic SNV or indel, or a homozygous deletion at a gene exon. Pairwise associations between gene mutation and SV signature activity were tested for using multiple linear regression, including gene mutation status, age at sampling, primary tumour site and tumour sample purity as independent variables. Genes were considered if they were mutated in at least 1% of tumours. TP53 mutation is associated with increased CIN, and TP53 was therefore included in all models. The Yeo–Johnson extension to the Box–Cox transformation was applied to mutation numbers to reduce heteroscedacity and to ensure distributions were approximately normal¹⁰⁸. Samples with missing independent variable values were excluded. Owing to mutational burden heterogeneity, only MSS primary tumours were considered in this analysis. P values were adjusted for multiple testing using Bonferroni correction and a threshold of P = 0.05 considered significant.

Characterizing copy number signatures

SigProfilerExtractor was used to extract copy number (CN) signatures in the 1,765 tumours with profiled CNAs¹². Where Battenberg identified a subclonal CNA, the copy number states corresponding to the largest tumour cell fraction were used, as SigProfilerExtractor cannot consider subclonal copy number states. Each copy number segment was assigned to 1 of 48 categories using SigProfilerMatrixGenerator, considering heterozygous or homozygous state, total copy number and segment length^12,17. Combinations of 1–30 de novo signatures were extracted and the recommended solution was accepted, balancing cosine distance with average stability (Supplementary Fig. 9), with the selection plot showing the mean sample cosine difference and average stability for de novo extraction of 1–30 CN signature. The accepted solution contained four de novo signatures.

De novo CN signatures were then deconvolved into their matching component COSMIC CN signatures from COSMIC (v.3) to identify six contributing COSMIC signatures, as shown below (CN1 (near-diploid state); CN2 (genome doubling); CN6 (chromothripsis/amplification with WGD); CN9 (CIN without WGD); CN17 (chromosomal-scale LOH); and CN20 (unknown aetiology)). CNV48A is a heterogeneous signature, dominated by heterozygous segments of 3–8 copies. It is decomposed into three COSMIC signatures: CN17, associated with HRD and TD (42.18%); CN6, associated with chromothripsis (29/72%); and CN20, which has a currently unexplained aetiology (28.1%). CNV48B is comprised primarily of heterozygous segments of 3-4 copies with a length of >40 Mb it is deconvoluted into a single cosmic signature CN2, associated with tetraploidy. CNV48C is dominated by heterozygous segments with a copy number of 2 and is decomposed to CN1, indicative of a diploid state. CNV48D is dominated by LOH segments with a copy number of 1 and heterozygous segments with a copy number of 2 and to a lesser extent 3–4, it deconvoluted into CN9, which has previously been associated with chromosomally unstable diploid tumours.

Each CN signature was assigned as being active or inactive in each sample. Associations with MSI status, ploidy and HRD status were calculated using Fisher’s exact test, comparing samples with and without the phenotype with those that had or did not have the active signature.

Predicting HRD

Evidence of HRD was assessed using HRDetect¹⁰⁹. HRDetect considers six genomic features predictive of HRD: (1) proportion of deletions with microhomology, (2) SBS3 contribution, (3) SBS8 contribution, (4) rearrangement signature RS3 contribution, (5) rearrangement signature RS5 contribution and (6) HRD index. HRDetect requires CNA data and was therefore run only on the 1,765 out of 2,023 tumours passing CNA calling. SBS3 and SBS8 contribution estimates were obtained from SigProfiler. Rearrangement signatures RS3 and RS5 were computed using HRDetect, using a previously reported rearrangement signature⁷³. Although HRDetect was trained on breast cancers, it has demonstrated high efficacy when applied to other cancer types¹⁰⁹. It was not possible to retrain HRDetect using our CRC samples, as few tumours exhibited a pathogenic germline BRCA1 or BRCA2 variant with somatic loss of heterozygosity of the wild-type allele.

Pathway analysis

Analysis of disrupted pathways

Altered pathways were identified by integrating coding and noncoding mutations using ActivePathways¹¹⁰. MSS, MSI and POL cancers were considered separately. Six mutation features were used: coding driver P values from IntOGen¹⁸ and 3′ UTR, 5′ UTR, core promoter, distal promoter and non-canonical splice site P values from OncodriveFML⁷⁸. We tested Reactome pathways obtained from MSigDB¹¹¹. All protein-coding genes included in at least one Reactome pathway were considered as the background gene set.

Driver mutation co-occurrence

Simple methods such as Fisher’s exact test and multiple regression were used to assess pairwise co-occurrence of driver genes. As these methods assume a null in which the probability of a gene alteration is independent of another gene, we also investigated use of the DISCOVER algorithm¹¹², which accounts for mutational heterogeneity at both the gene and tumour level. In practice, we reported simple association statistics, as we wished to include positively or negatively co-occurring driver genes or mutations, irrespective of a shared aetiology (for example, both genes containing short repeats prone to small indels in MSI tumours).

Cluster analysis

To search for groups of tumours with similar features, we used consensus clustering^113,114. We clustered 1,471 primary, treatment-naive tumours with CNA data using the following 304 clinical and molecular features: SNV, indel, SV and CNA burdens; all SBS, DBS, ID, SV and CN signature burdens; binary presence of mutations in 196 driver genes; ploidy; WGD status; fraction of the genome with LOH; mean ploidy across each chromosome arm divided by total ploidy (excluding the short arms of acrocentric chromosomes); age at sampling; sex; and subtype.

The features were ranked and normalized such that the resulting values were between zero and one. Hierarchical agglomerative clustering was run on these features using the diceR R package with the following distance metrics and linkage criteria:

• Distance metrics: Euclidean, Manhattan, cosine, correlation, Jaccard, eJaccard and fJaccard (from the R package proxy).

• Linkage criteria: average, complete, median, mcquitty, ward.D and ward.D2 (from R’s hclust function).

Each combination of distance metric and linkage criterion was run 10 times on random samples of 80% of the tumours. The number of clusters was varied from two to ten. We looked for robust clustering using the following criteria:

• The clustering must closely recapitulate the MSS, MSI, and POL subtypes

• High average clustering consensus¹¹⁴

• Absence of tiny clusters (<5 samples)

The ward.D2 linkage^115,116 consistently performed better than the other linkage criteria. With this linkage, Euclidean and Manhattan distances gave good clustering, but we chose Euclidean because the Manhattan distance failed to reproduce the POL subtype when the number of clusters was greater than six.

To increase the robustness of the clustering, we removed tumours that had an item consensus <0.7 and re-clustered using the resulting consensus matrix. This step removed 471 tumours that were difficult to cluster consistently and led to an increase in mean cluster consensus from 0.77 to 0.91. Following these steps, all samples had their subtype correctly classified, except for two MSI samples misclassified as MSS.

Immune profiling

HLA haplotyping

HLA typing of blood-derived normal samples was conducted using HLATyper, which is part of the Illumina Whole Genome Sequencing Service Informatic pipeline. The highest-ranking allele pair prediction for each type-I HLA allele (A, B and C) was taken to define a six-allele HLA set for each case.

Immune-escape prediction

We predicted three separate mechanisms of immune escape: (1) HLA gene mutation; (2) HLA gene LOH; and (3) mutation and LOH of other APGs.

Somatic mutations in the HLA locus were predicted using POLYSOLVER¹¹⁷. First, alleles were converted to a POLYSOLVER-compatible format (lower case, digits separated by underscore) and outputted into a patient-specific winners.hla.txt file. Next, the POLYSOLVER mutation-detection script (shell_call_hla_mutations_from_type) was run on matched tumour–normal pairs to call tumour-specific alterations in HLA-aligned sequencing reads using MuTect¹¹⁸. Strelka (v.2.9.9)⁷⁰ was also run to detect short insertions and deletions in HLA-aligned reads, as it offers increased sensitivity over POLYSOLVER’s default caller. Finally, both SNVs and indels passing quality control were annotated with POLYSOLVER’s annotation script (shell_annotate_hla_mutations).

LOH at the HLA locus was predicted using LOHHLA¹¹⁹. The same winners.hla.txt files were used as input, with POLYSOLVER’s comprehensive deduplicated FASTA of HLA haplotype sequences as reference. A type-I allele of a patient was annotated as allelic imbalance (AI) if the P value corresponding to the difference in evidence for the two alleles was <0.01. Alleles with AI were further labelled as LOH if the following criteria held: (1) the predicted copy number of the lost allele was <0.50 with CI < 0.70; (2) the copy number of the kept allele was >0.75; and (3) the number of mismatched sites between alleles was >10.

We also evaluated somatic mutations and copy number status of the following APGs¹²⁰: B2M, CALR, CANX, CIITA, ERAP1, ERAP2, HSPBP1, IRF1, PDIA3, PSMA7, PSME1, PSME2, PSME3, TAP1 and TAP2. First, somatic mutations were annotated using ANNOVAR¹²¹. An APG was deemed mutated if it contained any non-synonymous, frameshift, stop-loss, or stop-gain mutation in its exons. The copy number status of each gene was evaluated using Battenberg output.

A sample was defined as immune escaped if it showed at least one of the following: (1) HLA mutation; (2) HLA LOH; or (3) APG mutation. HLA AI was not considered to provide immune escape as AI can arise from multiple sources (including subclonal LOH and unequal focal gains of the locus) and therefore the effect of AI on antigen presentation is uncertain. For cases when HLA alterations could not be fully evaluated (see ‘Sample subsetting and statistical analysis’ below), but no HLA or APG alteration was detected, the immune escape status was considered unknown as we could not eliminate the possibility of immune escape.

Neoantigen prediction

We predicted neoantigens using NeoPredPipe, a Python-based pipeline combining ANNOVAR and netMHCpan (v.4.0)^122,123,124. In brief, all somatic SNVs and indels were annotated using ANNOVAR and for all non-synonymous exonic mutations the mutated peptide sequence was predicted. We took any 9- and 10-mer spanning the mutated amino acid (or acids), resulting in either (1) a 19-amino acid window for SNVs or (2) a peptide until the next predicted stop codon for frameshift mutations. These peptides were evaluated according to their novelty and predicted binding strength to the patient’s six-allele HLA set comprised of the HLA-A, HLA-B and HLA-C genes. Peptides that were new compared with the healthy human proteome with binding rank of two or below (among the best 2% of binders compared with a set of random peptides) were reported as neoantigens. All patient-specific HLA alleles were used for neoantigen prediction, regardless of mutation or LOH status of the HLA locus.

We considered a mutation neoantigenic if at least one of its downstream mutated peptides was a neoantigen with respect to any of the patient’s six HLA alleles. We defined neoantigen burden as the total number of neoantigenic mutations in the sample. We also evaluated the following alternative measures: (1) number of peptide–HLA binding pairs; (2) number of strong binder (best 0.5% of peptides) peptide–HLA binding pairs; (3) number of neoantigenic mutations in genes expressed in CRC (expression ≥10 TPM in ≥10% of TCGA CRCs)⁷. We found that all these measures were highly correlated with our definition of neoantigen burden: (1) R = 0.993, (2) R = 0.989, (3) 0.983; P < 10^–16 for all (Supplementary Fig. 10).

Sample subsetting and statistical analysis

Eighty-five samples were excluded from neoantigen calling because netMHCpan was unable to predict at least one of their HLA haplotypes. Overall, 217 samples had 1 or more haplotypes incompatible with POLYSOLVER, for which HLA mutation and LOH calling was restricted to the compatible haplotypes (1, 2 and 3 haplotypes were excluded in 171, 37 and 9 samples, respectively). In addition, LOH was not considered for 15 patients because they were homozygous for all type-I HLA genes. In total, 1,744 out of 2,023 samples had complete neoantigen and HLA alteration information available.

As CRC subtypes (MSS, MSI and POL) have substantially different mutation and immune properties, all analyses were completed separately for each subtype. Pairwise comparisons were conducted using Wilcoxon tests. Analysis of immune differences associated with tumour site was restricted to MSS primary samples, with samples that lacked specific information (site information missing or only specified as ‘colon’) excluded, leaving n = 1,100 samples.

Multivariate regression between immune escape types and neoantigen burden was performed using the lm function against the logarithm of neoantigen burden and therefore defined the fold change in burden associated with each escape type. Multivariate regression, including clinical characteristics, was carried out similarly, using the logarithm of total mutation burden as an additional independent variable. The number of POL samples was insufficient for statistical analysis and the regression analyses were therefore only conducted for MSS and MSI tumours.

PHBR analysis

We computed the immunogenicity of a given mutation in a given patient using PPHBR⁴⁰, which takes into account all novel peptides produced by that mutation and all HLA alleles present in the patient. Low PHBR values correspond to mutations that are likely to be presented on the cell surface and hence with a high immunogenic potential, whereas high PHBR mutations are less immunogenic. The overall immunogenic potential of a mutation within a cohort is defined as the median of PHBR values within that cohort. For each mutation and HLA haplotype pair considered, we generated all 8–11-mers overlapping the mutation and evaluated their binding affinity to the HLA allele using the ‘all-predictions’ mode of NeoPredPipe. The best (lowest) rank was recorded. For a given patient, PHBR were computed as the harmonic mean of six best rank values corresponding to the patient’s six HLA haplotypes (homozygous alleles were counted twice). We computed PHBR values for all single nucleotide mutations located in driver genes that were present in at least four cancers in the cohort. The 85 samples with incompatible HLA alleles were excluded.

To evaluate the effect of HLA alterations on PHBR values, we repeated the same analysis for affected patients with a reduced set (<6) of HLA alleles that were unaltered. To measure the level of patient- (HLA-) dependent selection on driver genes, we compared PHBR values for mutations in these genes between patients that did not carry the mutation and patients that did. Negative values indicate that mutations of the gene are enriched in patients for whom they have lower immunogenic potential. PHBR values between patients with no mutations and patients with mutations were compared using Wilcoxon rank-test, and P values were adjusted for multiple testing using Benjamini–Hochberg correction.

For comparisons such as those shown in Extended Data Fig. 7, the immunogenic potential of individual mutations was quantified using the median of PHBR values associated with that single nucleotide change for each patient belonging to a specific cohort or subcohort. Immunogenicity of groups of driver genes (for example, metastasis-specific drivers) was evaluated by considering all mutations observed in the genes and median PHBR computed across the entire cohort of CRCs or MSS primary CRCs, as indicated. Values for an individual mutation across different cohorts were compared using paired Wilcoxon rank-tests.

Mitochondrial genome characterization

Calling mitochondrial somatic SNVs and indels

Somatic mitochondrial SNVs and indels were called using Mutect2 (v.4.1.4.1)¹²⁵, with the light strand as reference based on the human mtDNA revised Cambridge reference sequence (rCRS). Somatic mitochondrial variants were excluded if they had the following:

• Low mapping quality score (<20).

• Low base quality score (<20).

• An alternative allele frequency <1%.

• Missing alternative reads in any stand direction.

• Location within hypermutated regions (302–316, 514–525 or 3106–3109).

Mutational distributions of SNVs, categorized by the six possible pyrimidine substitution classes, were constructed to analyse mutational processes. Distributions of substitutions on the D-loop, including and excluding variants between the two origins of replication (\({\text{O}}_{H}\) and Ori-b, between sites 16,197 and 191) were also analysed by substitution class⁶⁵. Pathogenic variants were identified using ClinVar¹⁰⁷, considering annotations where at least one submitter provided an ‘interpretation with assertion criteria and evidence’.

Mitochondrial copy number estimation

Autosomal and mitochondrial genome coverage was computed using fastMitoCalc¹²⁶. Using estimated sample purity (\(\rho \)), tumour ploidy (\(\varphi \)) and mean coverage depth, tumour sample mitochondrial DNA copy number was estimated as previously described³¹:

Mitochondrial copy number was estimated for only the 1,765 out of 2,023 tumours that passed CNA calling and therefore had purity and tumour ploidy estimates.

Linear regression was used to correlate mtDNA copy number with age at sampling, tumour stage, site of primary tumour, sex and tumour purity. The Yeo–Johnson extension of the Box–Cox transformation was applied to mtDNA copy number. Linear regression was applied considering all tumours and segregating MSS and MSI tumours. Regression results were adjusted for multiple testing using the Benjamini–Hochberg procedure.

Selection of mitochondrial mutation and POLG correlation

For the 13 mitochondrial protein-coding genes, selective pressure was quantified by calculating the respective dN/dS values using the R package dNdScv, with non-mtDNA chromosomes removed from the reference genome⁶. A global mitochondrial dN/dS value was also estimated, excluding MT-ND6 due to a suspected replication bias. Results were adjusted for multiple testing using the Benjamini–Hochberg procedure. In addition, it was investigated whether POLG mutations resulted in altered mitochondria mutational burden compared to other tumours. Only the primary MSI cohort was analysed for this trait, as other subcohorts had too few tumours with non-synonymous POLG mutations.

Genomic impact of previous treatments

Whether individuals had received systemic treatment or colorectum-targeting radiotherapy before sampling was based on data from NHSD and PHE-NCRAS. For NHSD, records related to systematic treatment were obtained from the Admitted Patient Care and Outpatients tables using associated Office of Population Censuses and Surveys (OPCS)-4 codes. For PHE-NCRAS, records related to systemic treatment were obtained from the AV_TREATMENT table using the event description codes, and from the Systemic Anti-Cancer Therapy (SACT) table. For PHE-NCRAS, records related to radiotherapy were obtained from the AV_TREATMENT table using the event description codes, and from the National Radiotherapy Dataset (RTDS) table considering records associated with a CRC diagnosis.

In total, 315 participants received systemic treatment or radiotherapy before tumour sampling. A total of 278 participants received systemic therapy before CRC sampling for sequencing, and information on the drugs administered was available for 182 of these participants. For 253 participants, the systemic treatment was used to treat CRC, whereas for 25 participants, it was used previously to treat another cancer. Overall, 94 participants received capecitabine, 23 received cetuximab, 93 received fluorouracil, 39 received irinotecan, 109 received oxaliplatin, 46 received steroids and 28 received other drugs. In total, 118 participants received colorectum-targeted radiotherapy before tumour sampling.

Associations between systemic treatment and colorectum-targeting radiotherapy before sampling with mutational signature activity were tested using multiple logistic regression. Previous treatment with radiotherapy, capecitabine, cetuximab, fluorouracil, irinotecan, oxaliplatin and steroids was included in the models as binary independent variants. Other treatments administered before sampling occurred in fewer than five individuals and were therefore not included in the models. One model was created for each of the identified SBS, ID, DBS and SV signatures, with signature presence encoded as a binary dependent variable based on whether any evidence of the signature was identified in each sample. In total, 96 samples that received treatment before sampling, but for which the specific administered drugs were unknown, were not included. Both primary tumours and metastases were considered in these analyses. Treatment coefficient P values were adjusted for multiple testing using Bonferroni correction and a threshold of P = 0.05, considered significant. Treatment duration was measured as the time between the first and last treatment administration.

Metastasis-specific analyses

Tumours were split between primary (n = 1,354) and metastatic (n = 105) MSS samples. Only MSS samples were included as there was just one MSI metastasis and no POL metastasis. Five primary tumours were matched to metastasis samples in this cohort, but for the purposes of the analysis all samples were treated as unmatched. To determine mutational burden, VCF files were filtered for PASS variants and the number of SNVs and indels summed. These were then divided by the total genome length (3,088.27 Mb). For the binned copy number analysis, the genome was first partitioned into 2,766 1 Mb windows. For each sample, the absolute allele-specific copy number within each bin was recorded. If two copy number segments overlapped a bin, the copy number of the segment with the larger overlap was recorded. Copy numbers were then classified according to the section ‘Classification of CNAs’. For each aberration type (gain or deletion/LOH) the proportion of primary tumours with that aberration was compared to the proportion of metastatic samples with two-sided Fisher’s exact tests. The difference between the proportions was then plotted as a trace along the genome with stars indicating significantly different bins. P values were corrected for multiple testing (FDR < 0.05). Absolute copy number calls were divided by mean integer ploidy to account for differences in ploidy between the two groups. The adjusted copy numbers for each bin were then compared between primaries and metastases using Wilcoxon signed-rank tests while correcting for multiple testing (FDR < 0.05). The difference in the mean (ploidy-adjusted) copy number was then plotted as a trace along the genome, with stars indicating significant bins.

Microbiome

Microbial identification

Microbial sequences¹²⁶ were identified using GATK PathSeq¹²⁷ aligned against the default PathSeq microbial genome bundles. A minimum clipped read length of 60 bp was used with all other parameters set to their defaults. Unambiguously assigned reads were used for the decontamination steps. Thereafter the adjusted score output was used, sharing ambiguous reads between species. Score output for each sample was converted to microbial cells per human cell for each taxon by adjusting for microbial and human average genome size (average human genome calculated from copy number and tumour cell percentage data).

This analysis showed that metastases had extremely low microbial content and therefore subsequent steps included only primary tumours unless otherwise stated. Reads passing PathSeq filters were realigned against the E. coli colibactin gene cluster¹²⁸ using bwa¹⁰², and matching reads counted.

Contaminants

Potential contaminant species were identified using methods developed by The Cancer Microbiome Atlas¹²⁹. In brief, the prevalence of species found in primary tumours and matched blood was compared (Extended Data Fig. 8a). Samples were called as positive for a species if two or more unambiguously aligned reads from the species was found. Species were deemed as probable tumour sample origin if a Fisher one-sided exact test found them to be more prevalent in the tumour sample than the blood sample (FDR < 0.05) and blood sample prevalence was <20% of samples. Genus level scores were recalculated from species scores by only including the species scores that survived this decontamination step. To mitigate the effects of species with mixed biological and contaminant components¹³⁰, downstream steps were adjusted for NHS Hospital Trust where possible (see below) as the processing laboratory was a plausible source of contamination.

Identifying taxa associated with CRC

CRC-associated taxa were identified by pooling all species level read numbers from eight published stool metagenomic studies^{2,131,132,133,134}. Application of LEfSe to these data identified 73 species and 37 genera associated with CRC¹³⁵. Bacterial species were classified as oral microbes if they were identified as ‘oral taxon’ or ‘oral species’ by PathSeq or if they were present in the expanded Human Oral Microbe Database¹³⁶.

Comparing microbiome and clinicopathological data

Microbial relative abundances were compared to clinicopathological data using decontaminated PathSeq output. Only tumours with complete data for the relevant categories were included in each comparison. Genus and species level alpha diversity was measured using the Shannon index and beta diversity using Bray–Curtis dissimilarity of relative abundance. Differences in beta diversity were measured by PERMANOVA using the adonis function¹³⁷ in Vegan using default settings, with permutations confined to within NHS Trusts using the ‘strata’ setting to minimize cross-site contamination differences. Taxa differing between clinicopathological categories were measured using MaAsLin2¹³⁸, with minimum abundance of 0, minimum prevalence of 0.1, and NHS Trust added as a random effect to minimize cross-site contamination differences.

Statistical analysis and clinicopathological correlates

Statistical tests were two-sided and unpaired unless otherwise stated. Fisher’s exact and χ² tests were used for categorical variables. Wilcoxon (rank-sum) tests, t-tests and Kruskal–Wallis tests were used for quantitative variables. Multivariable analyses are described below.

Correlating variables

Multiple linear regression was used to investigate the relationship between clinicopathological features and numbers of SNVs, indels, CNAs and SVs, and numbers of mutations attributed to SBS, ID, DBS and SV signatures. Number of CNAs was defined as the number of genome segments for which the clonal or subclonal copy number state was not 1:1 in non-WGD tumours or was not 2:2 in WGD tumours.

Multiple logistic regression was used to investigate the relationship between the presence or absence of clinicopathological features and driver gene mutation, recurrent arm-level CNAs, recurrent focal CNAs, WGD and evidence of CN signatures. Unlike SBS, ID, DBS and SV signatures, the activities of CN signatures do not represent numbers of mutations attributed to the signature¹². We primarily considered the presence or absence of CN signatures, but also assessed measures of activity or burdens where stated.

MSS primary and primary MSI tumours were considered separately. Signatures were tested if they were identified in at least 1% of the tumour set, driver genes were considered if they were mutated in at least 5% of the tumour set, and arm-level and focal copy number alterations were considered if identified as recurrent by GISTIC. TP53 mutation is associated with increased CIN, and TP53 somatic mutation status was therefore included in mutation number models. Considering multiple variables together in a single model is essential given that many of these variables are correlated, including age, primary tumour site and stage. The Yeo–Johnson extension to the Box–Cox transformation was applied to mutation numbers to reduce heteroscedacity and to ensure distributions were approximately normal¹⁰⁸. Samples with missing independent variable values were excluded. Primary tumour site and tumour stage were considered as ordinal variables. Primary tumour site was encoded as a single ordinal variable with the following values: caecum = 1; ascending colon = 2; hepatic flexure = 3; transverse colon = 4; splenic flexure = 5; descending colon = 6; sigmoid colon = 7, rectosigmoid junction = 8; rectum = 9. Exploratory analyses with location as a binary variable (proximal versus distal colorectum) or ternary variable (proximal colon, distal colon, rectum) were also performed in some cases (Extended Data Fig. 9b). Tumour stage was also encoded a single ordinal variable with values corresponding to the four Dukes stages. Unless otherwise stated, for each individual variable, P values were adjusted for multiple testing using Bonferroni correction and a threshold of P = 0.05 considered significant.

Survival analysis

Correlation of clinicopathological and genomic variables with all-cause mortality (overall survival) was assessed using Cox proportional hazards models. Follow-up time was measured from the date that the tumour was sampled (as a proxy for date of presentation or diagnosis) to the corresponding patient’s most recent time of contact. The median follow-up time was 1,075 days. Only individuals for whom the primary tumour was sequenced were included. To avoid proportional hazards assumption violation, individuals with MSS and MSI tumours were considered separately. Individuals were excluded if tumour sampling occurred before 1 January 2015 or the time between CRC diagnosis and tumour sampling was greater than 1 year. Hazard ratios were adjusted for sex, patient age at sampling, primary tumour location and Dukes stage. Owing to small numbers of deaths, Dukes stages A and B were combined. Analyses were performed regarding location as a binary variable (proximal versus distal colorectum) and as an ordinal variable (locations 1–9 from caecum to rectum).

After excluding individuals with missing covariate data, the MSS and MSI cohorts comprised 836 (144 deaths) and 272 (48 deaths) individuals. The following variables were analysed:

• Total mutational burden (SNVs and indels).

• SBS, DBS and ID mutational signature activity as binary indicators. Signatures were analysed if they were identified in <50% of tumours in the respective cohort.

• Immune escape status.

For analyses that required CNA profiles, smaller MSS and MSI cohorts comprising 810 (141 deaths) and 222 (40 deaths) individuals were used. The following variables were analysed using these smaller cohorts:

• Driver gene mutation status. Driver genes were considered mutated in a tumour if: (1) they contained an oncogenic mutation as defined by OncoKB and dNdScv annotation, (2) were homozygously deleted, or (3) were affected by a large copy number gain (total copy number state >5 for non-WGD tumours and total copy number state >10 copies for WGD tumours).

• WGD status.

• Chromosome-arm-level gains and deletions.

• Total SV number.

For each cohort, variables were only tested if at least 5% of deaths were present in each category. A variable was considered correlated with survival if it improved model fit using ANOVA and the z-test provided association evidence. The Benjamini–Hochberg procedure was used to determine FDR to adjust for multiple testing. Proportional hazards assumption violations were analysed for each test. In multiple Cox regression analysis, P = 0.05 was considered significant.

Normal colorectal epithelial cell signatures

Numbers and proportions of SNVs associated with each SBS signature were obtained from a previous study⁴⁴. For cases in which multiple crypts from the same colon region were sampled in a single individual, the median number and proportion of SNVs associated with each SBS signature was computed across these samples. For cases in which multiple crypts from the same colon region had been sampled in a single participant, the median number and proportion of variants attributed to each signature was considered. Supplementary Fig. 11 shows data from ref. ⁴⁴. IDA closely resembles ID18. P values were computed using Wilcoxon tests.

Software used

Supplementary Table 1 lists software versions used in this study and their URLs.

Reporting summary

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