GoKawiilPatients, therapy data and whole-genome sequencing
Patient enrolment, sample and data collection, and whole-genome sequencing were done independently for each cohort3,4,5. For KiCS, samples were sequenced on Illumina HiSeq X with target depth of 30× for normal tissue and 30× or 60× for tumours. Samples from ZERO were also sequenced on Illumina HiSeq X with target depth of 30× for normal tissue and 60× or 90× for tumours. MSK samples were sequenced on Illumina NovaSeq 6000 with target depth of 50× for normal and 95× for tumours. For therapy data collection, all three programmes used the same template with detailed therapy-related fields to retrospectively extract and collect all the available treatment details from the patient charts (Supplementary Table 2). The deep panel sequencing (~1,000× depth, using a >800-gene panel) and processing was performed as described3.
Detection of somatic alterations
For all three cohorts, raw FASTQ files were aligned to the human reference genome (GRCh37d5) using BWA-MEM46. Somatic alterations were identified by comparing each tumour sample to its matched normal sample. For solid and central nervous system tumours, the normal sample was typically blood-derived. For haematologic malignancies, a skin biopsy was generally used as the normal. For KiCS, somatic single-nucleotide variants (SNVs) and small indels were identified using Mutect2 (GATK v4.1.3)47. CNAs in the genome were identified using PURPLE (v1.4.10)48 and SVs were called using gridss (v2.9; default parameters)49. For ZERO, somatic SNVs and small indels were identified using Strelka (v2.0.17)50 and CNAs were detected using PURPLE (v2.39)48, while SVs were called using gridss (v2.72)49. For MSK, SNVs were identified using Strelka2 (v2.9.1)50, Mutect2 (GATK v4.0.1.2)47 and CaVEMan (cgpCavemanWrapper v1.7.5)51. CNAs were detected using Battenberg (cgpBattenberg v1.4.0)52. Small indels were detected using Strelka2, Mutect2 and Pindel (cgpPindel v1.5.4)53, and filtered against a panel of 100 unmatched normals. For multiple callers, we used a consensus of minimum two out of three callers. Further, SVs were called using gridss (v2.2.2). All SNV and indel variants across the three cohorts were annotated using VEP (v3.4)54. We applied a multi-tier filtering strategy to retain high-confidence SV calls across all three cohorts. Variants were required to meet the following criteria: variant quality ≥500, mean mapping quality ≥20, assembly support from both the breakpoint and remote break-end, variant-supporting fragments ≥1, and either split reads ≥2 or read pairs ≥2. All variants required assembly-based support with precise breakpoint resolution. Additional filters included minimum SV length ≥50 bp (for non-break-end variants), variant allele frequency (VAF) ≥ 0.05, and exclusion of extreme strand bias (0.05 ≤ strand bias ≤ 0.95). For break-end variants representing interchromosomal translocations or complex rearrangements, both break-end pairs were required to independently satisfy all criteria to ensure robust structural variant calls. For more detail on whole-genome sequencing processing and variant filtering see refs. 3,4,5.
Mutational signature analysis
Plotting the distribution of samples mutation burden revealed two separate subsets of high and low-burden samples with separation line at 25,000 mutations (8.628 mut Mb−1; Extended Data Fig. 3). Therefore, to avoid stronger signals in high-burden hypermutators potentially masking the signatures in low-burden samples, we ran mutational signature de novo extraction separately on these two subsets. The low- and high-burden subsets consisted of 577 and 34 samples, respectively. A similar approach was adopted to extract DBS and ID signatures. Plotting the distribution of the samples’ DBS burden revealed two separate subsets of low and high-burden samples (444 and 166 samples, respectively). ID analysis revealed three subsets of low, intermediate and high burden levels (89, 487 and 35 samples, respectively; Extended Data Fig. 3).
SigProfilerMatrixGenerator (v1.1.31; default parameters)26 was used to generate SBS288, DBS78 and ID83 mutation matrices from our general cohort as well as clustered mutations. In addition to the classic 96 mutation channels (each single-base substitution and its immediate 5′ and 3′ flanking bases), the SBS288 mutational context includes the strand orientation for mutations occurring in the genic regions. Platinum-based drugs, as well as other chemotherapies, such as nitrogen mustards, thiopurines and various alkylating agents, are known to create DNA lesions preferentially repaired by transcription-coupled nucleotide excision repair (TC-NER) which targets lesions on the transcribed strand of active genes32,55,56. Because TC-NER is strand-specific, it introduces a transcriptional bias in mutation patterns. Thus, we used the SBS288 mutational context to enhance our ability to detect and interpret therapy-associated mutational signatures. All SBS, DBS and ID mutational profiles are described elsewhere26. SigProfilerExtractor (v1.1.3)57 was applied to the generated matrices to extract de novo mutational signatures from both the general cohort and the clustered mutation categories. This tool, which uses non-negative matrix factorization to decipher de novo mutational signatures and their relative activities in a cohort, was started using random initialization and was run for 100 replicates. The default values were used for all other parameters. Finally, SigProfilerExtractor decomposed the extracted mutational signatures into the COSMIC (v3.2) set of known signatures. In case of SBS288, the extracted signatures were first collapsed into SBS96 before decomposition into COSMIC signatures. SigProfilerAssignment (v0.1.1; default parameters)58 was used as the refitting approach for signature assignment.
Enrichment analysis
We devised a logistic regression model, with adjustment for age, sex and tumour purity, to infer associations between therapies and mutational signatures. The model (penalty = ‘l1’, solver = ‘liblinear’) estimated main effects for each drug (restricted to those administered to ≥5 tumours) as well as all pairwise interaction terms (that is, synergistic effects). To evaluate significance, we calculated the observed AUROC for each drug–signature association, followed by a permutation test. This test involved randomly shuffling the binary outcome labels (presence or absence of the signature) 5,000 times and recalculating the AUROC each time. We then computed a one-tailed empirical P value as the proportion of permutations in which the permuted AUROC was greater than or equal to the observed AUROC. Associations were considered significant if they met all three thresholds: P value < 0.05, AUROC ≥ 0.7 and regression coefficient ≥ 0.7 (OR ≈ 2).
RNA sequencing analysis
... continue reading