Pan-cancer landscape of loose ends

We next applied JaBbA v1 to 1,330 high-purity tumor and matched normal SRS profiles previously analyzed in Hadi et al.⁴ (see Methods for details), identifying 154,322 (clonal and somatic) junctions (median of 63 per tumor sample) and 48,835 somatic loose ends (median of 21 per tumor sample). The somatic loose end burden per sample varied across a 200-fold range and was correlated (Spearman R²=0.68) with the junction burden (Fig. 1c,d).

Junction breakends may be reciprocal, meaning that they are near (within 10kb) of another breakend with opposite orientation. Reciprocal breakends are usually copy-neutral (Fig. ​(Fig.1e,1e, top left) which makes them difficult to detect through classic CN analyses. JaBbA’s bookkeeping of mass balance across segments and junctions enables sensitive detection of reciprocal and nonreciprocal SVs at both copy-neutral and copy-altered genomic regions (Extended Data Fig. 4a–e). Across cancer, we found that most (85%) cancer junctions were both nonreciprocal and copy-altered (Fig. ​(Fig.1e,1e, bottom left). Such junctions can arise from inherently nonreciprocal SVs, such as simple deletions, or begin as reciprocal translocations that undergo subsequent loss or gain of one of the derivative alleles (Extended Data Fig. ​Fig.4f).4f). Like somatic junction breakends, somatic loose ends were predominantly (92%) copy-altered (Fig. ​(Fig.1e,1e, bottom right), although copy-neutral loose ends were also identified (Fig. ​(Fig.1e,1e, top right). Taken together, these results suggest that loose ends arise by breakage and repair mutational processes similar to those generating junction breakends.

Open in a separate window

Extended Data Fig. 4

Examples of copy-neutral and copy-altered SVs detected by JaBbA.

(a) Fraction of breakends recovered by JaBbA in simulated genomes, by breakend class (described in Extended Data Fig. 3). Box plot: line (median), body (IQR), whiskers (1.5 times IQR) (b-c) Examples of simple (inversion, translocation) and complex (chromoplexy) copy-neutral reciprocal SVs. (d,e) Examples of complex non-reciprocal SVs associated with CN alterations. (f) Example of a copy-altered reciprocal translocation. In each example in panels b-f, the top graph shows the JaBbA genome graph and the bottom graph shows the binned purity / ploidy transformed tumor read depth. Genome graph legend for plots b-f is shown in panel b.

Source data

Loose ends harbor repetitive and foreign sequences

To study the sequence context around loose ends, we defined a canonical axis originating at the loose end with coordinates increasing along the DNA strand whose 3′ terminus matches the side of a segment on which a loose end is found, which we refer to as the loose end’s ‘forward’ strand (Fig. ​(Fig.2a).2a). We next asked whether loose ends occurred preferentially at reference sequence repeats. Indeed, we found that unmappable bases were enriched near loose ends, most frequently LINE elements (Fig. ​(Fig.2b2b and Extended Data Fig. ​Fig.5a).5a). We next reasoned that some loose ends would result from the somatic fusion of mappable bases to unalignable sequences. Confirming this, we found a tumor-specific enrichment of repetitive and foreign sequences, including satellite and viral sequences, mated to reads on the forward (but not reverse) strand of somatic loose ends (Fig. ​(Fig.2c2c and Extended Data Fig. ​Fig.5b5b).

Open in a separate window

Fig. 2

Loose ends pinpoint missing cancer SVs.

a, Loose end coordinates are centered at each loose end and increase in the 5′ to 3′ direction along the forward strand. For a loose end arising from the right side of its associated reference genomic segment (that is, the side with larger reference genomic coordinates), the forward strand is the positive reference strand, that is, the strand with increasing reference coordinates along its 5’ to 3’ direction. Conversely, for a loose end arising from the left side of its associated reference genomic segment, the forward strand is the negative reference strand. b, Density of unmappable bases around loose ends (Methods). c, Density of uniquely mapping reads with unmapped (i.e. non-uniquely aligning) mates around loose ends. d, Subclassification of loose ends based on local assembly and consensus alignment (Methods). RC, reverse complement. e, Alluvial plot showing each loose end class (bottom row) and the mappability tier of the distal (middle row) and proximal (top row) ends of breakend sequences obtained through local assembly or consensus alignment. f, Alluvial plot comparing SRS and LRS breakend calls. The fraction of breakends identified by LRS only, SRS only and both platforms (LRS and SRS) is shown (left), stratified by whether the breakend was reciprocal to another breakend in the same sample (right). LRS breakends were taken from tumor-specific junctions found by at least two of four LRS SV callers (SVIM⁴⁹, cuteSV⁵⁰, Sniffles2 (ref. ⁵¹) and SAVANA⁵²). SRS breakends comprise junction breakends and loose ends in the JaBbA v1 genome graph. g, Stacked barplots showing the fraction of complex SVs called from genome graphs with versus without the addition of LRS junctions. DEL, deletion; DUP, duplication; TIC, templated insertion chain. h, Examples of a rigma (left) and pyrgo (right) identified by LRS and missed by SRS. Tracks from top to bottom show tumor haplotype reconstructions, a genome graph with SRS and LRS junctions, a genome graph with only SRS junctions, and purity- and ploidy-transformed read depth. Pyrgo and rigma boundaries are delineated by a purple line at the top of each plot.

Source data

Open in a separate window

Extended Data Fig. 5

Genomic context of loose ends.

(a) Density of unmappable bases overlapping RepeatMasker-annotated LINE elements (top) and repeat classes in standard loose end coordinates (b) Density of reads mated to viral sequences (top) and satellite repeats (bottom) in the vicinity of loose ends. See Fig. ​Fig.2a2a for definition of forward and reverse strands. (c) Enrichment of large (>1 Mbp) CN-unmappable segment boundaries near partially mapped breakends (red) in comparison to other loose ends (pink). (d) Fraction of breakends in CN-mappable regions of the genome not fully resolved by SRS per sample, grouped by tumor type. Box plot: line (median), body (IQR), whiskers (1.5 times IQR) (n = 48 AML, Acute Myeloid Leukemia; 20 BLCA, Bladder Urothelial Carcinoma; 127 BRCA, Breast Invasive Carcinoma; 15 CESC, Cervical Squamous Cell Carcinoma and Endocervical Adenocarcinoma; 57 COAD, Colon Adenocarcinoma; 186 ESAD, Esophageal Adenocarcinoma; 43 GBM, Glioblastoma Multiforme; 30 HNSC, Head and Neck Squamous Cell Carcinoma; 40 KICH, Kidney Chromophobe Carcinoma; 24 KIRC, Kidney Renal Clear Cell Carcinoma; 19 KIRP, Kidney Renal Papillary Cell Carcinoma; 41 LGG, Brain Lower Grade Glioma; 28 LIHC, Liver Hepatocellular Carcinoma; 46 LUAD, Lung Adenocarcinoma; 34 LUSC, Lung Squamous Cell Carcinoma; 86 MALY, Malignant Lymphoma; 170 MELA, Melanoma; 48 OV, Ovarian Adenocarcinoma; 70 PRAD, Prostate Adenocarcinoma; 37 SARC, Sarcoma; 25 STAD, Stomach Adenocarcinoma; 12 THCA, Thyroid Carcinoma; 34 UCEC, Uterine Corpus Endometrial Carcinoma). (e) Fraction of rearrangement breakends in tumor suppressor genes and oncogenes that are not fully resolved by SRS.

Source data

To identify distinct classes of repetitive SVs missing from SRS whole-genome profiles, we systematically classified tumor-specific sequences fused to each somatic loose end through assembly or consensus alignment (Fig. ​(Fig.2d2d and Methods). Overall, 55% of somatic loose ends showed evidence of tumor-specific fusion to a distal sequence. For over half of these (33% of somatic loose ends), the distal sequence aligned uniquely, indicating that these were fully mapped breakends missed by the initial junction caller (Fig. ​(Fig.2e)2e) (SvAbA¹⁸). In 23% of somatic loose ends (3% of detected breakends), the distal sequence was repetitive or foreign and could not be unambiguously placed on any reference (ambiguously mapped breakends; Fig. ​Fig.2e).2e). Finally, 45% of somatic loose ends (6% of detected breakends) did not map to any distal location (partially mapped breakends; Fig. ​Fig.2e).2e). Notably, partially mapped breakends were enriched in boundaries of large (>1-Mb) CN-unmappable regions (odds ratio (OR)=3.8; P<2×10⁻¹⁶) (Extended Data Fig. ​Fig.5c),5c), indicating that some represented CN changes shifted away from a CN-unmappable SV breakend (for example, centromeric breakends causing arm-level chromosomal changes).

Combining fully mapped breakends across both loose ends and junctions indicated that 91% of JaBbA v1 breakends could be uniquely mapped. Notably, the fraction of partially or ambiguously mapped breakends did not vary substantially across cancer types (Extended Data Fig. ​Fig.5d;5d; range of 5–33%) or established cancer drivers (Extended Data Fig. ​Fig.5e;5e; range of 0–38%), although we observed tumor types (for example, acute myeloid leukemia) and cancer genes (SMARCB1, TSC2 and FGFR3) with higher (>25%) fractional burdens. Given the estimated recall of JaBbA v1 (~96%), these results suggest that 87% of cancer SVs in the 87% of the genome that is CN-mappable can be fully resolved by SRS.

Long-molecule validation

To orthogonally assess these SRS-derived estimates of missing somatic SVs, we profiled the whole genomes of 11 melanoma (n=10) and breast cancer (n=1) tumor samples and their matched normal tissues with both SRS and Oxford Nanopore Technologies long-read sequencing (LRS; median read N50 of 11kb; median coverage of 73× and 32× for tumor and normal samples, respectively). After calling large (>10-kb) somatic SVs in CN-mappable regions (Methods), we found a strong overlap (87%, 7,258 breakends) between LRS and SRS breakends, including 77% overlap with fully mapped SRS breakends (Fig. ​(Fig.2f).2f). The majority of junction calls identified by either platform had local read depth changes that were consistent with breakend topology; reciprocal breakends were copy-neutral, whereas nonreciprocal breakends showed a CN drop along their forward strand (Extended Data Fig. ​Fig.6a).6a). This analysis along with manual inspection of long and short read support at inidivudal junctions (Extended Data Fig. ​Fig.6b)6b) suggested that both SRS-only and LRS-only junctions comprise largely true positives; combining SRS and LRS breakend counts suggests that SRS missed ~12% of breakends. This result is consistent with our simulation-based estimate of recall (Fig. ​(Fig.1b1b and Extended Data Fig. ​Fig.3f).3f). Notably, we found a similar proportion of reciprocal and non-reciprocal breakends among those detected and missed by SRS (Fig. ​(Fig.2f),2f), indicating that reciprocal and copy-neutral breakends do not comprise the bulk of missed structural variation in cancer genomes. These results confirm our SRS findings that most cancer SVs are nonreciprocal and copy-altered (Fig. ​(Fig.1e1e).

Open in a separate window

Extended Data Fig. 6

Comparison of breakends detected by SRS and either LRS or sLRS.

(a) CN changes around breakends (JaBbA v1 loose ends and junction breakends) across reciprocal (left) and non-reciprocal (right) breakends detected by SRS and/or LRS (see Fig. ​Fig.22 legend and Methods for more details). Bins are plotted in breakend coordinates (similar to loose end coordinates, oriented to the forward strand of the loose end or junction breakend) and normalized so that the relative CN of the segment harboring the breakend is 0. (b) Example of junctions called by SRS and LRS. Tracks from top to bottom show shared qnames per pair of genomic bins, SRS read pairs supporting the variant junction, the SRS-only genome graph, and the tumor read depth. The first example shows a call made by only short reads, the next shows a call made by both LRS and SRS, the rightmost panel shows calls made by only LRS. Of note, the rightmost deletion is missed by JaBbA as it spans a CN-unmappable region (red). (c) Alluvial plot showing high-confidence breakends called by either sLRS, SRS, or by both platforms, and whether the breakend is reciprocal or non-reciprocal. sLRS breakends are taken from tumor-specific junctions found by at least two of the three sLRS SV callers (LinkedSV⁵⁷, GROC-SV⁵⁸, and NAIBR⁵⁹). SRS breakends comprise junction breakends and loose ends in the JaBbA v1 genome graph. (d) Fraction of complex and simple SV calls made using sLRS and SRS junctions or SRS junctions only. (e) CN changes around reciprocal (left) and non-reciprocal (right) breakends detected by SRS and/or sLRS, similar to (a).

Source data

We next asked whether LRS improved SV event detection, which relies on the recognition of high-order patterns across multiple junctions^3,4. Although LRS did not help identify many additional simple or complex events relative to SRS (Fig. ​(Fig.2g),2g), LRS junctions also resolved breakends at complex SVs found by SRS, including for chromothripsis, pyrgo, rigma and templated insertion chains^3,4. The incorporation of LRS junctions enabled more complete haplotype reconstruction at loci where SRS found loose ends (Fig. ​(Fig.2h2h).

As additional validation of our results, we analyzed 27 high-purity (purity of >0.5) breast cancer and matched normal samples with both SRS and synthetic LRS (sLRS) whole-genome profiles (10x Genomics linked reads, median N50 molecule length of 23kb, median coverage of 173× and 98× in tumor and normal samples, respectively; Methods)¹⁹. Similar to LRS, most sLRS SV calls (Methods) overlapped with SRS breakends, showed concordant patterns of reciprocality and CN change, and yielded similar complex SV calls in sLRS junction-augmented genome graphs (Extended Data Fig. 6c–e). These breast cancer and melanoma LRS and sLRS results are consistent with our pan-cancer finding that SRS captures most large cancer SVs in CN-mappable regions.

Loose ends reveal neotelomeres

We next sought to investigate specific mutational processes engendering loose ends. We observed that a fraction (4.8%) of ambiguously mapped loose ends (0.01% of all breakends) were fused to telomere repeats, as evidenced by telomere repeat-positive sequences mated to reads on the positive loose end strand (Fig. ​(Fig.3a).3a). We refer to these breakends as telomere repeat-positive loose ends and surmised that they might represent neotelomeres, telomere-stabilized chromosome ends at previously interstitial genomic loci.

Open in a separate window

Fig. 3

Loose ends reveal neotelomeres.

a, Density of reads mated to telomere repeats near loose ends. TR, telomere repeat. b, Density of reads mated to GRTRs and CRTRs on the forward and reverse strands of GRTR⁺, CRTR⁺ and telomere repeat-negative loose ends. c, Potential etiologies of telomere repeats fused to loose ends. d, Density of sLRS barcodes harboring telomere repeat-positive read pairs near GRTR⁺ loose ends relative to telomere repeat-negative loose ends. Telomere repeat-positive read pairs are defined as a read pair in which one mate is spanned entirely by G-rich telomere repeats and the other by C-rich telomere repeats. e, Fraction of loose ends fused to a unique distal interstitial location via sLRS (Methods). Data are shown as mean±95% confidence interval (CI) (n=71 GRTR⁺ loose ends and 28 CRTR⁺ loose ends from 14 tumor samples). f, Estimated telomere length at GRTR⁺ loose ends compared to telomere repeat-negative loose ends (n=71 GRTR⁺ loose ends from 14 tumor samples). In box plots, the line represents the median, the body represents the IQR and whiskers extend to 1.5 times the IQR. g, Schematic of the neotelomere detection assay. h, Southern blot showing a diffuse ~23-kb band in U2OS cells but not in a control cell line (Saos-2). Both cell lines show the 4.4-kb HindIII control band. i, The U2OS-specific band disappears after exonuclease (Bal-31) digestion of genomic DNA before HindIII digestion (left) without altering the overall size distribution of DNA fragments (right). EtBr, ethidium bromide. Time refers to length of Bal-31 exposure. For h,i, the experiment was repeated four times with similar results. Panels show uncropped images of the entire gel lanes. j, Fraction of tumors with a neotelomere (that is, GRTR⁺ loose end) across the given categories. ‘ATRX low’ corresponds to ATRX RPKM of <500 and ‘TERT low’ corresponds to TERT RPKM of 0. Error bars are the 95% CIs on the binomial proportion. LOF, loss of function. k, Tumor type enrichment of GRTR⁺ loose ends. Tumor types with a statistically significant association with GRTR⁺ loose end burden are highlighted in red. See the Fig. ​Fig.1d1d legend for definitions of the abbreviations. In j,k, P values were calculated by two-sided Wald’s test on the coefficients of a negative binomial generalized linear model (Methods). In k, values with log(OR)>log(1.5) and false discovery rate (FDR)<0.1 after Benjamini–Hochberg correction are highlighted in red.

Source data

Telomere repeat-positive mates were found on the forward strand of telomere repeat-positive loose ends, but not on the reverse strand or in matched normal samples (Fig. ​(Fig.3a),3a), indicating that these were neither telomere insertions^20,21 nor constitutional neotelomeres^22,23. Deeper analysis of telomere repeats at loose ends revealed strong strand bias, with loose ends harboring either G-rich (GRTR) or C-rich (CRTR) repeats but not both (Fig. ​(Fig.3b).3b). The GRTR pattern is consistent with a neotelomere, whereas the CRTR pattern is consistent with the fusion of an interstitial sequence to a native chromosome end (Fig. ​(Fig.3c,3c, right). The predominance of the GRTR pattern among telomere repeat-positive loose ends, in combination with the tumor specificity and forward strand bias, suggested that somatic neotelomeres are frequent in cancer.

To better assess sequences fused to GRTR⁺ loose ends, we profiled three cancer cell lines (U2OS, NCI-H526 and NCI-H838) with sLRS (Methods). We found telomere repeat-positive linked reads within 5kb of 26 of 31 GRTR⁺ loose ends (83.8%) (Methods). Telomere repeat-positive linked reads were found up to 50kb upstream of each GRTR⁺ loose end, indicating power to map distal fusion partners at these loci (Fig. ​(Fig.3d).3d). In contrast to sLRS junctions and telomere repeat-negative loose ends, linked reads at GRTR⁺ loose ends rarely (<1.5%) mapped to distant chromosomal locations, consistent with new chromosome ends (Fig. ​(Fig.3e).3e). Quantitative analysis of repeat counts at linked reads mapping to these loci (Methods) revealed 2.4±1.3 (s.d.) kb of telomere repeats per GRTR⁺ locus, in line with previous estimates of native cancer telomere lengths²⁰ (Fig. ​(Fig.3f3f).

To confirm that GRTR⁺ loose ends were indeed chromosome ends, we performed Southern blot analysis on restriction-digested U2OS and control (Saos-2) genomic DNA using radiolabeled probes against two U2OS GRTR⁺ loose ends. At each locus (Fig. ​(Fig.3g3g and Extended Data Fig. ​Fig.7a),7a), we found a small (<5-kb) band consistent with an unaltered reference allele and a longer U2OS-specific diffuse band consistent with a neotelomere (Fig. ​(Fig.3h3h and Extended Data Fig. ​Fig.7b).7b). To further investigate the nature of these nonreference bands, we subjected intact genomic DNA to exonuclease (Bal-31) digestion²⁴. The U2OS-specific (but not wild-type) bands disappeared with prolonged exonuclease exposure (Fig. ​(Fig.3i3i and Extended Data Fig. ​Fig.7c),7c), consistent with their origin at a chromosome end. These results establish these two U2OS GRTR⁺ loose ends as bona fide neotelomeres.

Open in a separate window

Extended Data Fig. 7

U2OS neotelomere validation.

(a) Schematic of neotelomere detection assay. Genomic DNA was subjected to EcoRI digestion followed by Southern blotting with a radiolabeled probe to the site of a chr 10 GRTR+ loose end (bottom track). Alleles that lack a neotelomere at this location will show a focal≈5 kbp control band. A neotelomere will yield a diffuse band that is sensitive to digestion by the exonuclease Bal-31 prior to EcoRI digestion. (b) A Southern blot shows a diffuse high molecular weight band in U2OS but not in a control cell line (SAOS-2), consistent with a neotelomere. Both cell lines show the≈5 kbp EcoRI control band. (c) The U2OS-specific band disappears when double stranded DNA ends of intact genomic DNA are digested by Bal-31 for over 5 minutes prior to EcoRI digestion and probe hybridization, confirming a new chromosome end (left). Bal-31 digestion does not substantially alter the overall size distribution of ethidium bromide (EthBr) labeled DNA fragments in the EcoRI digest (right). Time refers to length of Bal-31 exposure. Experiment repeated three times with similar results. Panels show uncropped images of entire gel lanes.

We next hypothesized that telomerase-mediated healing of double-stranded DNA breaks might give rise to neotelomeres (Fig. ​(Fig.3c,3c, left)²⁵. However, neotelomeres were not found more frequently in tumors that amplified TERT or expressed it at high levels (CN>2 ploidy, expression z score>2). Instead, neotelomeres were enriched in samples with low or negligible TERT expression (reads per kilobase per million mapped reads (RPKM)=0) (Fig. ​(Fig.3j).3j). Tumors that lack telomerase may activate the alternative lengthening of telomeres (ALT) pathway, a break-induced replication (BIR) process (Fig. ​(Fig.3c,3c, middle) suppressed by ATRX²⁶. Indeed, we found that neotelomeres were significantly more common in tumors harboring truncating mutations in ATRX than in ATRX-wild-type cancers (Fig. ​(Fig.3j).3j). Furthermore, we found that several ALT-associated cancers, including sarcomas (18%; OR=6.47; P=1.95×10⁻⁵) and low-grade gliomas (12.3%; OR=3.92; P=4.1×10⁻³), had the highest rate of GRTR⁺ loose ends relative to other tumor types (Fig. ​(Fig.3k).3k). These results indicate that GRTR⁺ loose ends and neotelomeres may be a new hallmark of the ALT phenotype.

Loose ends link viral integration to amplicon formation

Surveying additional mutational processes engendering loose ends, we found ambiguously mapped somatic breakends fused to viral sequences, indicating junctional viral integration at large SVs (Extended Data Fig. ​Fig.8a).8a). While the integration of viral sequences into otherwise unrearranged loci (Extended Data Fig. ​Fig.8a,8a, left) has been widely studied in cancer^27,28, the role of viruses in causing chromosomal-scale SVs (Extended Data Fig. ​Fig.8a,8a, right) has been a topic of only recent interest^29–31. Somatic loose ends harboring tumor-specific viral sequence (viral loose ends) were rare overall (~1% of cancers), although enriched in cancer types with viral etiology in our dataset⁴: cervical squamous cell carcinoma (CESC; 32%), liver hepatocellular carcinoma (LIHC; 13%) and head and neck squamous cell carcinoma (HNSC; 7%) (Extended Data Fig. ​Fig.8b).8b). Consistent with previously characterized viral integration patterns, we found viral loose ends fused to oncogenic HPV sequences in CESC and HNSC and hepatitis B virus (HBV) sequences in LIHC²⁷.

Open in a separate window

Extended Data Fig. 8

Viral loose ends are amplified in virus-driven cancers.

(a) Schematic contrasting viral loose ends and insertions. Virus integrations causing SVs will give rise a tumor-specific viral contig associated with a pair of distant loose ends (right) while virus insertions arise from viral integration into unrearranged DNA (left). (b) Bar plot showing fraction of samples containing at least one viral loose end by tumor type. Error bars show mean fraction±95% confidence interval. Only tumor types with ≥10 samples are shown (n = 48 AML, 20 BLCA, 127 BRCA, 15 CESC, 57 COAD, 186 ESAD, 43 GBM, 30 HNSC, 40 KICH, 24 KIRC, 19 KIRP, 41 LGG, 28 LIHC, 46 LUAD, 34 LUSC, 86 MALY, 170 MELA, 48 OV, 70 PRAD, 37 SARC, 25 STAD, 12 THCA, 34 UCEC). See Extended Data Fig. 5d caption for expanded tumor type abbreviations. (c) Association between viral loose ends and high-CN amplicons. Left: Viral loose ends are more likely to have high CN (≥ 7) compared to all other loose ends (n = 32 viral loose ends, 68354 other loose ends). Error bars show mean fraction±95% confidence interval. (d) CN of viral loose ends separated by virus type. Box plots: line (median), body (IQR), whiskers (1.5 times IQR). P-values calculated with two-sided Wilcoxon rank sum test. (e-f) Example viral loose ends from TCGA HNSC sample 4077 (e) and CESC sample A0TN (f). From top, “walk" representation of proposed circular allele; alignments of read pairs with discordant alignments between human and viral genome; JaBbA v1 genome graph showing two loose ends corresponding to the discordant read alignments; binned purity / ploidy transformed read depth. Genome graph legend in panel e applies to both e and f.

Source data

Breakends initiating complex amplifications are themselves likely to be amplified⁴. Viral loose ends were frequently amplified (CN>7) relative to nonviral loose ends (P=1.7×10⁻⁴; OR=8.66) (Extended Data Fig. ​Fig.8c),8c), and HPV-16 loose ends had higher mean CN than either HPV-18 or HBV loose ends (P=8.2×10⁻³ and P=2.2×10⁻⁵, respectively, Extended Data Fig. ​Fig.8d).8d). Among these was an HNSC tumor (TCGA-4077) locus where two high-copy viral loose ends on chromosome 14 flanking an intronic region of the RAD51B gene were fused to opposite ends of the HPV-16 genome (Extended Data Fig. ​Fig.8e).8e). This locus is consistent with an ecDNA where HPV-16 is fused between two ends of a long-range duplication junction. This and other similar amplicon structures with high-copy viral loose ends (Extended Data Fig. 8e,f) point to HPV-16 integration as an initiating event in SV evolution, rather than a viral insertion into an existing ecDNA.

Crossover between parental homologs is rare in cancer

We next asked whether loose ends could be used to assess the contribution of aberrant homologous recombination to cancer rearrangements. Homologous recombination-driven crossover between parental homologs (allelic homologous recombination, or AHR) is a hallmark of meiosis³². Although AHR has been observed in somatic cells³³, its contribution to cancer structural variation is unclear. AHR crossovers lead to segmental uniparental disomy (UPD) in approximately half of segregants (Fig. ​(Fig.4a,4a, left). In balanced allelic graphs, AHR crossovers manifest as reciprocal pairs of partially mapped and copy-neutral loose ends on distinct parental homologs (Fig. ​(Fig.4b,4b, left, and Methods). Notably, this form of UPD (AHR-UPD) is mechanistically distinct from UPD arising through progressive acquisition of nonhomologous recombination (for example, end joining)-driven rearrangements and/or chromosomal missegregation (progressive UPD, or P-UPD; Fig. 4a,b, right).

Open in a separate window

Fig. 4

AHR rarely drives CN-mappable breakends.

a, Schematic showing mechanistic differences between AHR and P-UPD, two mechanisms that give rise to segmental UPD. b, Examples of AHR-UPD (left) and P-UPD (right). The allelic graph (top subpanel) shows parental homolog-specific CN, which matches purity- and ploidy-transformed allelic SNP read counts (scatterplot, second subpanel from top) (Supplementary Information). The AHR-UPD locus shows no breakends in the total CN JaBbA v1 graph (third subpanel) but a pair of loose ends in the allelic graph. By contrast, the P-UPD locus does not harbor a pair of allelic graph loose ends, but rather contains a copy-altered breakend in both the allelic and total CN JaBbA v1 graphs. Het, heterozygous. c, Width distribution of segments produced by AHR-UPD, P-UPD and all other LOH (n=545 AHR-UPD ranges, 39,877 P-UPD ranges and 61,469 other LOH ranges from 1,330 tumors). In box plots, the line represents the median, the body represents the IQR and whiskers extend to 1.5 times the IQR. d, Fractional contribution of P-UPD, AHR-UPD and other forms of LOH to the total number of LOH segments. e, Schematic of NAHR. f, Number of estimated (y axis) versus true (x axis) NAHR-mediated breakends per simulated sample (n=500 simulated genomes). The blue line shows the line of best fit, with Pearson’s correlation coefficient provided on the graph; error bands show the standard error of the prediction. The P value was calculated from the t distribution of Pearson’s correlation coefficient test statistic. g, Fraction of somatic junctions, somatic loose ends and germline loose ends consistent with NAHR rearrangements in the SRS pan-cancer whole-genome cohort (n=1,330 samples). Error bars represent the 95% CIs on the binomial proportion. Germ, germline; Som, somatic. h, Fraction of germline and somatic LRS junctions, SRS junctions and SRS loose ends consistent with NAHR in a separate melanoma and breast cancer cohort with paired SRS and LRS whole-genome profiles (n=11 samples). Error bars represent 95% CIs on the binomial proportion.

Source data

In our simulations (Extended Data Fig. ​Fig.3a3a and Methods), JaBbA v1 distinguished AHR-UPD from P-UPD with both high precision (84.4%) and high recall (87.4%), substantially outperforming previous allelic CN algorithms (with precision ranging from 11–44%) (Extended Data Fig. 9a,b). Analysis of segment width distributions showed that AHR-UPD was distinct from P-UPD, whose distribution closely mirrored that of other forms of loss of heterozygosity (LOH; Fig. ​Fig.4c).4c). Likewise, AHR-UPD events were large (median width of 19.8Mb), unlike P-UPD events (median width of 0.69Mb) and other forms of LOH (median width of 0.62Mb), which were focal (Fig. ​(Fig.4c4c).

Open in a separate window

Extended Data Fig. 9

Benchmarking AHR detection.

(a) Confusion matrix comparing AHR-UPD and P-UPD labels assigned using the allele-specific CN profile produced by various CN inference algorithms on a simulated dataset of 200 samples. AHR-UPD calls were defined as segments with major allele CN = 2 and minor allele CN = 0, immediately adjacent to a heterozygous segment with total CN = 2. The remaining segments with major allele CN = 2 and minor allele CN = 0 are annotated as P-UPD. The calls made by each CN inference algorithm are shown as rows, while the ground truth call is shown as columns. The number in each box represents the Gbp of each call across all simulated samples, while the color scale is normalized by row (representing precision). Only segments with width>10 kbp were considered in this analysis. (b) Overall precision, recall, and F1 score of AHR calls made by each algorithm. (c) Prevalence of AHR in each cancer type. Only tumor types with ≥10 cases are included. (n = 48 AML, 20 BLCA, 127 BRCA, 15 CESC, 57 COAD, 186 ESAD, 43 GBM, 30 HNSC, 40 KICH, 24 KIRC, 19 KIRP, 41 LGG, 28 LIHC, 46 LUAD, 34 LUSC, 86 MALY, 170 MELA, 48 OV, 70 PRAD, 37 SARC, 25 STAD, 12 THCA, 34 UCEC). See Extended Data Fig. 5d caption for expanded tumor type abbreviations. Error bars show mean fraction per tumor type±95% confidence interval.

Source data

Although AHR was found in many cancers (24% of all tumors) and specific tumor types (for example, 55% of cases of malignant lymphoma) (Extended Data Fig. ​Fig.9c),9c), it contributed to a minority of UPD events, most of which were progressive (31% P-UPD versus 1% AHR-UPD by total width) (Fig. ​(Fig.4d).4d). Overall, a small minority of detected cancer breakends (<1%) arose by AHR (including non-UPD LOH). On the basis of an approximate rate of 0.5 AHR events per tumor and 100 cell divisions in the average ancestral cancer clone, and barring effects of selection, we estimate a rate of 10⁻¹² AHR events per base pair per cell division. This is four orders of magnitude lower than the rate of meiotic recombination in human gametes, suggesting that AHR events are infrequent in somatic evolution³⁴.

Germline but not somatic loose ends are consistent with NAHR

A second mechanism by which aberrant homologous recombination can cause large SVs is through non-AHR (NAHR), or crossover between long (>500-bp) stretches of nearly identical genomic sequences at distant haploid coordinates^32,35,36 (Fig. ​(Fig.4e).4e). We reasoned that such SVs would engender pairs of loose ends with substantial (>500-bp) strand-specific sequence homology in their vicinity (Extended Data Fig. 10a and Methods)³⁶. Indeed, the burden of homologous loose end pairs accurately reflected the true NAHR burden across a compendium of simulated SRS tumor whole-genome profiles (Extended Data Fig. ​Fig.3a)3a) harboring a wide range of NAHR SV fractions (1–10%) (Fig. ​(Fig.4f4f).

Open in a separate window

Extended Data Fig. 10

Detecting NAHR in SRS, LRS, and sLRS whole genome profiles.

(a) Schematic of method for detecting NAHR at SRS junction or loose end associated breakends. Locus 1 and 2 are sites of breakends (blue lines) either belonging to a junction or a pair of loose ends in a tumor sample. A loose end was designated putative NAHR if there was at least one other loose end within the sample sharing>500 bp of strand-specific sequence homology (>96% identity) within 10 kbp of each breakend (see Methods). Similarly, a junction was annotated as putative NAHR if>500 bp homologous sequences were found within 10 kbp of each breakend belonging to the junction. (b) Fraction of sLRS junctions (n = 3352 somatic, 1978 germline) and SRS junctions (n = 4740 germline, 5694 somatic) and loose ends (n = 282 germline, 2397 somatic) consistent with NAHR, separated by whether the rearrangement was germline or somatic (n = 27 tumors). Error bars show mean fraction±95% confidence interval. Germ, germline; Som, somatic. (c) Example of a somatic LRS junction that was designated putative NAHR. Tracks from top to bottom show the number of shared reads with split alignments between pairs of loci in the tumor and matched normal LRS, purity / ploidy transformed read depth in the tumor and matched normal, the JaBbA v1 graph, and fraction of bases having exact homology. (d) Example of two germline loose ends identified as putative NAHR and validated as a germline NAHR junction in the LRS data (top heatmaps). Genome graph legend as in panel c. (e) Example of two somatic loose ends identified as putative NAHR in the sLRS data. Top two tracks show number of shared sLRS barcodes in the tumor and matched normal. Genome graph legend as in panel c. (f) Mean per-megabase density of somatic breakends on each autosome in CN-mappable regions across n = 1330 tumor genomes. Box plot elements: line (median), body (IQR), whiskers (1.5 IQR). (g) Pearson correlation between eligible homologous position pairs per chromosome and total estimated somatic NAHR breakend count per chromosome. Error bands show standard error of the prediction and P-value was calculated from t-distribution of Pearson’s correlation coefficient test statistic.

Source data

Analyzing breakend pairs within each tumor, we found that approximately 20% of germline loose ends (Methods) were consistent with NAHR in contrast to only about 0.5% of somatic loose ends (and 0.06% of all somatic SV breakends) (Fig. ​(Fig.4g).4g). These findings are consistent with prior observations about the substantial role of NAHR in germline variation^8,37. The somatic NAHR burden did not vary by tumor type nor was it lower in tumors harboring biallelic pathogenic mutations in DNA repair genes, including frequently mutated homologous recombination pathway mediators (BRCA1, BRCA2, PALB2 and RAD51C). In summary, given a mean of 0.16 somatic NAHR events per tumor occurring across an estimated eligible territory of 2.8×10⁸ homologous position pairs, we estimate a somatic NAHR density of 6×10⁻¹⁰ events per cancer genome bp² (Methods).

To validate these SRS findings in long-molecule whole-genome profiles, we analyzed 38 melanoma and breast cancer cases profiled with SRS and either LRS or sLRS. Both LRS and sLRS data confirmed our SRS findings that somatic NAHR SVs were rare (<1% of LRS junction calls) while germline NAHR SV events were common (Fig. ​(Fig.4h4h and Extended Data Fig. 10b–e). Notably, we did not identify any reciprocal somatic NAHR rearrangements, a class of SVs that may potentially be missed through analysis of SRS loose ends.

JaBbA v1 algorithm

The JaBbA v1 algorithm builds on the previous JaBbA (v0.1) algorithm introduced in Hadi et al.⁴ with two key modifications: (1) updating the JaBbA statistical model to a Laplace distribution, which improved performance and convergence, and (2) balancing allelic genome graphs across parental SNP homologs to enable breakend phasing and identification of allelic loose ends. These and other pipeline changes are visually summarized in Extended Data Figs. ​Figs.1a1a and 2a–d. The updated algorithm is described in further detail below.

JaBbA v1 pipeline

In addition to algorithmic improvements, the JaBbA v1 pipeline includes additional data processing improvements compared to the previous version⁴: (1) correction of sample-specific bias in tumor read depth and (2) rigorous definition of CN-mappable regions. Unlike previously⁴, 1-kb binned read depth x is obtained via dryclean⁵³, a robust principal-components analysis-based algorithm to remove systematic low-rank biases in binned read depth using a panel of normal samples (Supplementary Note 2). In addition, we mask bins that occur in CN-unmappable regions of the genome (see below for details) and use purity and ploidy estimates to transform read depth into CN units (Supplementary Note 2). The JaBbA v1 pipeline then applies CBS⁵⁴ to x and takes the union of the resulting segment endpoints with SvAbA²² junction breakends to construct a preliminary genome graph.

In practice, we use three iterations of total CN MIP optimization (equation (3)) followed by allelic mass balance. After each total CN iteration, the results are processed to yield a simplified graph where reference adjacent segments are merged if a loose end or variant junction with CN>0 does not exist at their interface. The first MIP iteration takes as input only large (>10-kb) and high-confidence (FILTER = PASS) SvAbA junctions and CBS segment breakends. Clusters of high-confidence reciprocal SVs are constrained into the model (that is, by including in the set E_F in equation (3)). The second MIP iteration augments the graph from the first iteration with low-confidence SvAbA junctions located within 10kb of fitted loose ends. The final MIP iteration refits the graph from the second iteration but adds a noninteger chromosome-specific offset that prevents hypersegmentation from small inaccuracies in purity and ploidy estimation or subclonal CN changes. Allelic mass balance is then run on the balanced genome graph output of the final MIP iteration. To optimize AHR detection, before allelic mass balance, large (>1-Mb) segments of the balanced genome graph are further split by running CBS on minor allelic count vectors (see Supplementary Note 2 for details of chromosomal bias correction and AHR detection).

Mappability analysis

We performed exhaustive self-alignment of all 101-mers in the GRCh37 reference to identify base-unmappable positions, that is, those whose corresponding 101-mer gave rise to at least one full-length supplementary alignment or harbored at least one masked (N) base. A position was then called CN-unmappable if more than 90% of the bases in a 1-kb window around that position were base-unmappable; otherwise, it was called CN-mappable. An analogous approach was used to determine GRCh38 and 150-bp mappability (Extended Data Fig. ​Fig.2d).2d). Additional details are provided in Supplementary Note 2.

Short-read whole-genome sequencing

SRS whole-genome profiles for 1,330 high-purity (>0.5) tumors (inferred sex: 586 female, 744 male; age: 4–90 years, 164 unknown) and 326 cell lines (provided sex: 139 female, 170 male, 17 unknown; age: 0.25–74.05 years, 46 unknown) were obtained from a previous study published by Hadi et al.⁴. Additional SRS whole-genome libraries were prepared using the Illumina TruSeq DNA PCR-free Library Preparation Kit and profiled on an Illumina NovaSeq 6000 sequencer with 2×150-bp cycles. Following GRCh37 alignment, data processing and standard whole-genome variant calling, we ran the JaBbA v1 pipeline (see above) and previously published somatic CN callers (JaBbA v0.1 (ref. ⁴), ASCAT v2.5.2 (ref. ¹⁴), FACETS v0.6.2 (ref. ¹⁷), Sequenza v3.0.0 (ref. ¹⁶) and TITAN v1.28 (ref. ¹⁵)). Additional details regarding SRS library preparation, data processing and variant calling are provided in Supplementary Note 3.

Long-read whole-genome sequencing

LRS profiles were generated for ten melanomas and one triple-negative breast cancer collected at MSKCC (collection details above; inferred sex: six female, five male; age: >17 years). Following high-molecular-weight (HMW) DNA extraction, LRS was performed on the Oxford Nanopore Technologies PromethION sequencer using R10 chemistry with two flow cells per tumor and one flow cell per normal sample. Following GRCh37 long-read alignment (minimap2 v2.17), LRS SV junction calls were identified from the two-way consensus of four callers: cuteSV (release v2.0.2; ref. ⁵⁰), SAVANA (release 0.2.3; ref. ⁵²), SVIM (release 2.0.0; ref. ⁴⁹) and Sniffles2 (release v2.0.7; ref. ⁵¹). Callers were run on tumor and normal samples separately (cuteSV, SVIM, Sniffles2) or in paired mode (SAVANA). Tumor and normal junction calls with identical orientation and 1-kb padded overlap were merged across algorithms. Additional details regarding LRS library preparation and data processing are provided in Supplementary Note 4.

Synthetic long-read whole-genome sequencing

sLRS whole-genome profiling was performed on 25 breast cancer tumor–normal pairs (collection details above; inferred sex: 25 female, 0 male; age: >17 years) and a panel of 8 ATCC cell lines (provided sex: 5 male, 3 female; age: unknown) previously profiled with SRS by the Cancer Cell Line Encyclopedia⁵⁵. In brief, HMW DNA was subjected to 10x Genomics Chromium Genome library preparation and sequenced on an Illumina NovaSeq 6000 sequencing system to approximately 30× base and 170× physical coverage. 10x Genomics linked reads were aligned to GRCh37 using the EMerAld aligner (v0.6.2)⁵⁶. To nominate SV junctions, we applied a consensus of three algorithms (LinkedSV, https://github.com/WGLab/LinkedSV (commit 1b77a14)⁵⁷; GROC-SV, https://github.com/grocsvs/grocsvs (v0.2.6)⁵⁸; NAIBR, https://github.com/raphael-group/NAIBR (commit 15eba96)⁵⁹) run on tumor and normal sLRS alignments. Tumor and normal junction calls with identical orientation and 1-kb padded overlap were merged across algorithms. Somatic SVs were then called as junctions found in tumors by two or more algorithms and undetected in the normal sample. Additional details regarding sLRS profiling and data processing are provided in Supplementary Note 4.

Short- versus long-read platform comparisons

We used 1-kb strand-specific overlap of SRS breakends (junctions and loose ends) and LRS/sLRS junction breakends to assess concordance between SV calling platforms. To assess the ability of LRS or sLRS junctions to resolve SRS loose ends, we applied an additional iteration of junction balance to SRS-derived balanced genome graphs, including additional LRS or sLRS junctions as input. We then overlapped loose ends in the original SRS genome graph with junctions incorporated into the LRS/sLRS genome graph. If a loose end was within 1kb of an LRS breakend or within 10kb of an sLRS junction breakend on the same strand, we considered that loose end to have been resolved by LRS/sLRS. We applied gGnome⁴ to annotate and compare complex SV events across SRS, LRS and sLRS JaBbA graphs and used genomic overlaps to identify shared versus platform-specific calls.

Loose end classification

To identify candidate distal mappings for loose ends, we used Fermi⁶⁰ local assembly (https://github.com/mskilab-org/RSeqLib) and realignment of loose end-associated reads and mates. Fermi contigs were assessed for tumor and normal read support through BWA realignment of reads to contigs and the reference (https://github.com/mskilab-org/readsupport) to uncover tumor-specific contigs with distal alignments. To find additional distal mappings, we also analyzed consensus distal alignments of loose end-associated reads. We then labeled loose ends as ‘fully mapped’ or ‘ambiguously mapped’, respectively, if they had a unique or ambiguous tumor-specific distal mapping and ‘partially mapped’ otherwise. See Supplementary Note 5 for full details of loose end classification.

Neotelomere analysis and validation

We identified telomere repeat-positive sequences as those matching one of a panel of G-rich and C-rich telomere repeat trimers and their six cyclic permutations. A loose end was considered telomere repeat-positive if a tumor-specific telomere repeat-positive contig (see above) was found at the loose end. Given an sLRS telomere repeat-positive loose end, we counted read pairs comprising exclusively telomere repeats across all sLRS barcodes associated with the locus and multiplied the maximum value by the median intramolecular distance between reads across all molecules in that sLRS library to estimate the neotelomere length. To validate neotelomere candidates, genomic DNA was isolated from U2OS and Saos-2 cells and, where indicated, treated with Bal-31 exonuclease²⁴. Bal-31-digested DNA was isolated by phenol extraction and ethanol precipitation and was then digested with the appropriate restriction enzyme. Gel electrophoresis of the DNA, Southern blotting and hybridization with Klenow-labeled radioactive probes were performed. See Supplementary Note 5 for additional details of neotelomere analysis and validation.

Nominating NAHR junctions

We identified all pairs of sequences with ≥500bp of homology (96% sequence identity) through exhaustive BWA-mem self-alignment of all 101-mers on both strands of GRCh37. Loose ends b₁ and b₂ were annotated as a putative NAHR junction if a sequence of ≥500bp within 10kb of b₁ on its forward strand was found to be homologous to a sequence within 10kb of b₂ on its negative strand (Extended Data Fig. 10a). Similarly, junctions were annotated as NAHR if their breakends b₁ and b₂ demonstrated the above property.

Distinguishing mechanisms of UPD

After allelic CN inference using JaBbA or other tools, UPD segments (total CN=2 and minor allele CN=0) were identified. UPD segments reference adjacent to another segment of CN=2 without LOH (minor allele CN=1) were called AHR-UPD; otherwise, segments were called P-UPD.

Simulating tumor and normal SV profiles

We simulated 500 SRS whole-genome SV profiles on GRCh37 by rearranging the fully phased NA12878 Platinum genome⁶¹. To simulate phased rearrangement junctions, we randomly sampled and shifted pan-cancer SvAbA junctions⁴ and assigned each a random NA12878 haplotype. We also simulated NAHR junctions at a prevalence of 0.1–10% by linking pairs of homologous positions in the genome. Junction breakends were used to define allelic segments, and both were assigned a phased integer CN (‘balance’ function, gGnome⁴) with a target ploidy. We sampled junctions to simulate imperfect sensitivity for junction detection (accounting for sampling effects due to finite read depth or stromal admixture) at a rate proportional to tumor purity. Realistic purity and ploidy values were sampled from pan-cancer distributions⁴. To simulate read depth at bins or SNPs, the normalized purity-adjusted CN of each 1-kb bin was multiplied by a coverage factor to achieve a target genome-wide per-base read depth (80× in tumors and 40× in normal samples) and a bias factor was computed from normalized read counts for that bin in a random normal diploid sample. This product was used as the mean parameter for a Poisson distribution, which was sampled to obtain the final total (or allelic) read depth. See Supplementary Note 6 for additional simulation details.

Benchmarking

To benchmark breakend detection, we compared the endpoints of simulated ‘ground-truth’ breakends to CN calls from JaBbA v1 (this paper), JaBbA v0.1 (ref. ⁴), ASCAT (v2.5.2)¹⁴, FACETS (v0.6.2)¹⁷, Sequenza (v3.0)¹⁶ and TITAN (v1.28)¹⁵. True-positive breakends were defined as those found within 10kb of a ground-truth breakend on the same strand. We applied a similar approach to assess true-positive rates among JaBbA v1 versus v0.1 loose ends. To assess the accuracy of total CN inference, we computed the root mean square error between estimated and ground-truth total (or allelic) CN values across 10-kb genomic bins. See Supplementary Note 6 for additional benchmarking details.

Identifying NAHR-eligible sites in T2T CHM13 v2

We sampled 1 million 500-bp substrings from T2T CHM13 and realigned them to T2T CHM13 using BWA-mem, identifying all alignments with cigar 500M. This yielded nearly 9 million position pairs (p₁,p₂), where p₁ and p₂ represent the starting coordinate of the query and alignment, respectively. We then divided the self-alignments into three categories (CNU×CNU, CNM×CNU and CNM×CNM) on the basis of the overlap of p₁ and p₂ with a CN-mappable region lifted from GRCh37 to T2T CHM13. See Supplementary Note 7 for details of T2T CHM13 NAHR analysis.

Estimating the genome-wide unmappable breakend fraction

To extrapolate SRS findings from the CN-mappable to the CN-unmappable genome, we applied two principles. First, NAHR rearrangements occur in proportion to the number of homologous position pairs in the genome. Second, non-NAHR rearrangements (including AHR and non-HR SVs including end joining) occur in proportion to the number of bases in the genome. Our CN-mappable NAHR analysis found 216 somatic NAHR events across 2.8×10⁸ NAHR-eligible position pairs in 1,330 genomes, giving a somatic NAHR density of 6×10⁻¹⁰ events per bp². Given the 2.7×10¹⁰ NAHR positions in the 13% of the genome that is CN-unmappable, we estimate an approximate NAHR burden of 16.2 breakends per tumor genome. Given the 681 AHR and 357,000 non-HR CN-mappable breakends in the 1,330 tumor samples, we estimate 0.6 AHR and 310 non-HR events per genome. Putting these numbers together, we estimate that ~17% of large SV breakends occur in CN-unmappable regions, and, given JaBbA’s 96% recall in CN-mappable regions, 80% of large SV breakends will be detected and 73% will be fully resolved by SRS. Given an estimated HR burden of 16.8 SVs per tumor genome, the fractional HR SV burden is approximately 5%. See Supplementary Note 7 for additional details of these calculations.

Statistics and reproducibility

All statistical analysis was performed as stated in the figure legends using the R programming language (v4.0.2). Statistical methods were not used to predetermine sample size. The study design did not involve blinding or randomization. See Supplementary Note 8 for additional details on statistics and reproducibility as well as loose end association analyses.

M.I., J.M.B., X.Y, A.D., J.R. and H.T. were supported by M.I.’s Burroughs Wellcome Fund Career Award for Medical Scientists, Weill Cornell Medicine Department of Pathology and Laboratory Medicine startup funds and/or NYU Perlmutter Cancer Center startup funds. Z.-N.C. was supported by an F30 predoctoral fellowship from the NIH/NCI (F30CA268747) and a Medical Scientist Training Program grant from the National Institute of General Medical Sciences of the NIH under award number T32GM007739 to the Weill Cornell/Rockefeller/Sloan Kettering Tri-institutional MD PhD Program. K.H. was supported by an NIH/NCI F31 graduate research fellowship (F31CA232465). T.d.L. was supported by R35CA210036. J.M.B., M.I. and T.d.L. were additionally supported by Starr Cancer Consortium award I13-0019. M.I., S.N.P. and J.R.-F. were additionally supported by Starr Cancer Consortium award I11-0051. We thank the members of the Imieliński laboratory for help with manuscript proofreading. We thank C. Black and the New York Genome Center ResComp team for high-performance computing support. We thank S. Dider and T. Dey in the Imielinski lab for software engineering support. The funders had no role in study design, data collection and analysis, decision to publish or preparation of the manuscript.