Distinct Classes of Complex Structural Variation Uncovered across Thousands of Cancer Genome Graphs

Abstract

Cancer genomes often harbor hundreds of somatic DNA rearrangement junctions, many of which cannot be easily classified into simple (e.g., deletion) or complex (e.g., chromothripsis) structural variant classes. Applying a novel genome graph computational paradigm to analyze the topology of junction copy number (JCN) across 2,778 tumor whole-genome sequences, we uncovered three novel complex rearrangement phenomena: pyrgo, rigma, and tyfonas. Pyrgo are "towers" of low-JCN duplications associated with early-replicating regions, superenhancers, and breast or ovarian cancers. Rigma comprise "chasms" of low-JCN deletions enriched in late-replicating fragile sites and gastrointestinal carcinomas. Tyfonas are "typhoons" of high-JCN junctions and fold-back inversions associated with expressed protein-coding fusions, breakend hypermutation, and acral, but not cutaneous, melanomas. Clustering of tumors according to genome graph-derived features identified subgroups associated with DNA repair defects and poor prognosis.

Keywords: aneuploidy; cancer evolution; cancer genomics; chromothripsis; fragile sites; genome graphs; mutational processes; phasing; structural variation; superenhancers.

Figure 1.. Junction-balanced genome graphs enable complex structural variant characterization

(A) Schematic of elevated junction copy number (JCN) from the duplication of an allele harboring a DEL-like junction, resulting in a characteristic read depth and junction pattern from which a junction-balanced genome graph can be reconstructed. (B) Junction balance analysis (JaBbA) integrates binned WGS read depth data and junctions to estimate JCN and generates coherent models of rearranged genome structure. The topology of JCN can be used to define complex SV events as variant subgraphs. (C) The coverage of junction-spanning 10X linked-read barcodes in HCC1954 (y axis) correlates highly with the JaBbA JCN estimate, obtained through the analysis of HCC1954 short-read WGS (x axis). (D) Top, heatmap of the number of 10X linked-read barcodes shared between each pair of 1 kbp genomic bins. Bottom, JaBbA genome graphs within 400 kbp of the featured junctions. (E) Tumor type sample counts across 2,778 genome graphs. See Table S1 for abbreviations. * marks datasets with multiple samples. (F) Correlation of the purity and ploidy corrected read depth difference vs. JaBbA-fitted JCN. (G) Histogram of JCNs and their associated categories in the cohort. See also Fig. S1 and Tables S1–S3

Figure 2.. Rigma and pyrgo are novel patterns of clustered low copy rearrangements

(A) Quantile-quantile (Q-Q) plot of observed vs. expected −log₁₀(P) quantiles obtained from a gamma-Poisson model of low-JCN DUP-like junction density across genomic bins (see STAR Methods). Red dots indicate sample-specific model outliers (FDR<0.5). Top right, an example of a pyrgo-associated outlier window with a high density of low-JCN DUP-like rearrangements within a sample. Bottom right, a non-outlier window containing a simple duplication event. (B) Right, Q-Q plot similar to (A), but low-JCN DEL-like junctions. Top left, an example of a rigma-associated outlier window containing a high density of low-JCN DEL-like junctions within a sample. Bottom left, a non-outlier window containing a simple deletion event. (C-D) Fraction of samples within tumor types that harbor pyrgo and rigma, respectively. Significantly enriched tumor types (compared to all others) are marked by asterisks (Fisher’s exact test). Significance levels: *** (P < 1 × 10⁻³), ** (P < 0.01), * (P < 0.05). (E-F) Association of events with replication timing (ordinal logistic regression), superenhancers (Hnisz et al., 2013) (logistic regression), known fragile sites (Kumar et al., 2019) (logistic regression), and gene width (Wilcoxon rank-sum test, Bonferroni corrected P-values < 0.05). Error bars on bar plots represent 95% confidence intervals on the Bernoulli trial parameter. See also Fig. S2

Figure 3.. Rigma junctions evolve gradually and in trans

(A) Top, comparison of the total genomic territory covered by chromothripsis and rigma footprints. Bottom, the fraction of rearrangements that occur in cis (i.e. on the same predicted haplotype) when the longest possible contigs are inferred from the JaBbA graph. (B) Linked-read (LR) sequencing and local allelic deconvolution of a rigma found in a lung adenocarcinoma cell line (NCI-H838, see STAR Methods) demonstrates DEL-like junctions occur on separate haplotypes (i.e. in trans). (C) A comparison of the fraction of events that occur early (i.e. in multiple samples from the same patient) for simple deletions and rigma in the BE cohort, which has a median of 4 samples per patient across 80 patients. (D) Reconstruction of haplotypes across multiple samples from a single case in the BE WGS dataset. P-values obtained by Wilcoxon rank-sum test in (A) and Fisher’s exact test in (C). RD, read depth. Error bars on bar plots represent 95% confidence intervals on the Bernoulli trial parameter. See also Fig. S3

Figure 4.. Analysis of amplified subgraphs identifies tyfonas

(A) Framework to identify features of complex amplicons through the analysis of amplified subgraphs. Each subgraph is annotated for three features: number of elevated (JCN>3) junctions, maximum JCN relative to the highest vertex CN in the subgraph, and total JCN of all fold back junctions relative to maximum JCN. (B) Hierarchical clustering of subgraph features reveals three stable clusters, representing distinct amplification SV classes: (C) double minutes, (D) BFBCs, and a new event pattern, (E) tyfonas. Top track in C-E is the JaBbA-estimated CN. Bottom, normalized read depth data. (F) Tumor types significantly enriched in each amplicon class. Significance determined by Fisher’s exact test comparing each tumor type against all others. All enrichments with FDR < 0.25 are shown. “Other” category denotes non-significant tumor types in the respective analysis. (G) Tyfonas enrichment in melanoma and sarcoma subtypes. (H) Enrichment of expressed protein-coding fusion transcripts by count (left) and density per Mbp of event territory (right) in tyfonas relative to other amplicon types and chromothripsis (Wilcoxon rank-sum test). Error bars on bar plots represent 95% confidence intervals on the Bernoulli trial parameter. See also Fig. S4

Figure 5.. Tyfonas are enriched in breakend hypermutation outside of SBS 2/13

(A) Breakend-centered coordinate system to analyze mutational patterns near junctions. Top, the cis (+ coordinates) sides of the breakend (i.e. attached to the junction) have been fused across the rearrangement junction (red-colored line), while the trans (− coordinates) sides are disconnected from this derivative allele. Bottom, relative SNV density is the CN corrected and 101 bp smoothed density of SNVs at every base pair on this axis normalized to the average SNV density between −5 and 0 kbp. (B) Relative SNV density of SBS 2/13 and all other SBS associated contexts near breakends demonstrates a peak in the first 1 kbp on the cis side of the breakend. (C) Normalized breakend density for SBS 2/13 (top) and other SBS (bottom) SNVs for each of amplified event types. Enrichment P-values obtained from Wald test by gamma-Poisson regression (see STAR Methods). Significance determined by Bonferroni correction at a threshold of < 0.05. See also Fig. S5

Figure 6.. Chromosomal integration of a MYB / MYCN -associated tyfonas in a small cell lung cancer cell line.

(A) Multi-platform profiling of a tyfonas in the small cell lung cancer cell line NCI-H526. (B) Top, heatmap of Hi-C in cell line NCI-H526 demonstrating MYCN and MYB co-amplification and contiguity with chromosome 1. Blue arrows on the Hi-C heatmap highlight the pixels supporting contiguity of tyfonas fragments with chromosome 1. Middle, JaBbA graph of the NCI-H526 tyfonas, where junction input was derived from short-read WGS and Bionano genomics optical mapping. The locations of MYCN, MYB, and the tyfonas BAC probes used for FISH experiments are shown. Bottom, Hi-C profiling of the normal cell line GM12878, obtained from (Rao et al., 2014). (C) Relative expression of MYCN and MYB across the Cancer Cell Line Encyclopedia (Ghandi et al., 2019b). (D) Metaphase FISH of NCI-H526 using a chromosome 2 probe targeting the center of the tyfonas near MYCN (red) and a second probe targeting the chromosome 2 centromere (green) (scale bars, 5 μm). (E) Candidate reconstruction of a linear allele spanning the tyfonas amplicon. The coordinates of the reconstructed allele (y axis) are shown using a nonlinear axis which begins at a junction adjacent to chromosome 1, and repeatedly threads between the MYCN and MYB genes, which are both amplified by the event. See also Fig. S6

Figure 7.. Genome graph-derived features define biologically distinct patient groups

(A) Heatmap of normalized junction burden across 2,487 unique patients yields 14 clusters named after their most prevalent event types (CT, chromothripsis; TYF, tyfonas; BFB, BFBC; BR, BFBC and rigma; SPRS, sparse; CP, chromoplexy; TIC, templated insertion chain; TRA, translocations; DM, double minute; PYR, pyrgo; DUP, various duplications; DDT, deletion, duplication, and TIC; INVD, inverted duplications). (B-C) Significantly enriched tumor types within and outside of selected clusters. Blue (vs. gray) bars indicate the fraction of cases within (vs. outside of) the cluster that have the given feature, e.g. tumor type (B) or mutation (C). Significance in (B) and (C) determined by Bayesian logistic regression (Wald Test), with significant results (FDR < 0.1) for 742 hypotheses in (B) and 4,774 hypotheses in (C) (see STAR Methods). (D) Kaplan-Meier curves across BR, PYR and TYF clusters. P-values obtained via log-rank test. (E) Cox proportional hazard analysis correcting for age, tumor type, sex, metastasis status, tumor mutational burden, SV junction burden, and TP53 status demonstrating confidence bounds on the adjusted hazard ratio for 13 clusters relative to the QUIET cluster. Significantly associated variables (FDR < 0.1, 13 hypotheses) colored in red. Error bars on bar plots represent 95% confidence intervals on the Bernoulli trial parameter. See also Fig. S7 and Tables S3, S4