Haplotype-specific assembly of shattered chromosomes in esophageal adenocarcinomas

Abstract

The epigenetic landscape of cancer is regulated by many factors, but primarily it derives from the underlying genome sequence. Chromothripsis is a catastrophic localized genome shattering event that drives, and often initiates, cancer evolution. We characterized five esophageal adenocarcinoma organoids with chromothripsis using long-read sequencing and transcriptome and epigenome profiling. Complex structural variation and subclonal variants meant that haplotype-aware de novo methods were required to generate contiguous cancer genome assemblies. Chromosomes were assembled separately and scaffolded using haplotype-resolved Hi-C reads, producing accurate assemblies even with up to 900 structural rearrangements. There were widespread differences between the chromothriptic and wild-type copies of chromosomes in topologically associated domains, chromatin accessibility, histone modifications, and gene expression. Differential epigenome peaks were most enriched within 10 kb of chromothriptic structural variants. Alterations in transcriptome and higher-order chromosome organization frequently occurred near differential epigenetic marks. Overall, chromothripsis reshapes gene regulation, causing coordinated changes in epigenetic landscape, transcription, and chromosome conformation.

Graphical abstract

Figure 1

Chromothripisis present and sequencing modalities used in this study (A–F) Rearrangement plots exhibiting chromothriptic regions in each sample. Black dots denote copy number. Vertical lines denote breakpoints: black represents translocations, green represents deletion orientation, teal represents tandem-duplication orientation, dark blue represents tail-to-tail inversion orientation, and brown represents head-to-head inversion orientation. HCM-SANG-0311-C15 has two chromothriptic chromosomes, and all other samples have one chromothriptic chromosome. (G) Graphical overview of study design. Patient-derived esophageal adenocarcinomas and matched normal blood were used to generate data. Genomic sequencing was used to generate assemblies upon which epigenetic datasets were layered.

Figure 2

Genome assemblies of complex rearrangements (A) Main steps in method for cancer genome assembly. (B and C) Left: dot plot alignments of the wild-type and chromothriptic assemblies, respectively, against the reference GRCh38 chromosome 6. A sequence identical to the reference GRCh38 genome would appear as a 45° diagonal line. Right: alignment plots produced using XMatchView for largest contigs from the wild-type and chromothriptic assemblies, respectively. Black line represents chromosome 6 from GRCh38, and the line below represents the contig. Blue regions are direct repeats, and pink regions are inverted repeats. An identity threshold was set at 90%. (D, G, and J) Rearrangements plots are as previously described. (E, F, H, I, K, and L) Dot plots alignments of each haplotype to the reference GRCh38 genome. Hap, haplotype; CT, chromothripsis; WT, wild-type.

Figure 3

Haplotype-resolved structural rearrangements, retrotransposons, and epigenetic modalities (A) Total retrotransposon classes in each sample. L1, LINE-1; SVAs, SINE-VNTR-Alus. (B) Total retrotransposon types in each sample. (C) On the chromosomes with evidence of chromothripsis, number of retrotransposon calls. (D) Total SV counts greater than 1 kb, filtered to remove germline SVs and SVs in repeat regions, which are likely to be erroneous calls. (E and F) SV size distribution for wild-type and chromothriptic chromosomes, respectively. (G) Proportion of reads assigned in different data types.

Figure 4

Differences in chromatin accessibility and histone modification between alleles (A and B) Log2 raw read counts of peaks that can be assigned in both samples. Marginal density plots are shown. (C) Heatmap of differentially bound and open chromatin sites on the wild-type and chromothriptic chromosome 6 in HCM-SANG-0300-C15. (D) Distance of peaks that are weaker on the chromothriptic chromosome to the nearest SV on the chromothriptic chromosome (left) and wild-type chromosome (right) relative to non-differential peaks. p values were calculated using the Wilcoxon rank-sum test. (E) Distance of peaks that are stronger on the chromothriptic chromosome to the nearest SV on the chromothriptic chromosome (left) and wild-type chromosome (right) relative to non-differential peaks. p values were calculated using the Wilcoxon rank-sum test. (F and G) Distance effect on fold change of peak strength split by mark profiled using ChIP-seq.

Figure 5

Differences in gene expression between alleles (A and B) Log2 raw expression of genes that can be assigned in both samples. Marginal density plots are shown. (C) Heatmap of 40 most differentially expressed genes on the wild-type and chromothriptic chromosome 6 in HCM-SANG-0300-C15. There is high concordance between biological repeats. (D) Distance of differential genes with higher gene expression on the chromothriptic allele to active peaks. The p value was calculated using the Wilcoxon rank-sum test. (E) Distance of differential genes with lower gene expression on the chromothriptic allele to active peaks. The p value was calculated using the Wilcoxon rank-sum test. (F) Distance of non-differential genes to active peaks. The p value was calculated using the Wilcoxon rank-sum test. (G) Distance of differential genes with higher gene expression on the chromothriptic allele to inactive peaks. The p value was calculated using the Wilcoxon rank-sum test. (H) Distance of differential genes with lower gene expression on the chromothriptic allele to inactive peaks. The p value was calculated using the Wilcoxon rank-sum test. (I) Distance of non-differential genes to inactive peaks. The p value was calculated using the Wilcoxon rank-sum test.

Figure 6

Differences in higher-order chromatin conformation between alleles (A and B) Example TAD boundary calls on wild-type and chromothriptic chromosomes, respectively. Top: Hi-C contact matrix with SV track showing where structural variants are present relative to the reference genome. Middle: TAD boundary calls using 150-kb bins. Signal values represent the difference in number of contacts between adjacent bins. When the differences are considered to represent TAD boundaries, orange dots are plotted. These correspond with small TAD structures. Bottom: TAD boundary calls using 1-Mb bins as described above. These represent larger TAD structures. (C) TAD sizes on wild-type versus chromothriptic chromosomes in HCM-SANG-0300-C15 called using different bin sizes. TADs are inferred as regions between boundary calls. p values were calculated using the Wilcoxon rank-sum test. Median values are stated on boxplots.

Figure 7

A 600-kb region of the custom assemblies Wild-type assembly on the left and chromothriptic assembly on the right. The regions shaded in gray are identical sequences in the two chromosomes, if ignoring indels, and SNPs and are the same as the sequence chr6:151,209,484–151,396,541 in the GRCh38 reference genome. This gray shaded region contains AKAP12 and ZBTB2. Colored blocks in the contiguous sequences track are contiguous sequences found in the reference genome. A block is contiguous but is not found in the reference GRCh38 genome adjacent to the next block. The arrows denote orientation of the regions relative to the reference GRCh38 genome. Each pixel in the Hi-C track is an interaction. Differential ChIP-seq and ATAC-seq peaks were called using DiffBind, and the scale measures relative peak score. The peak is shown on the chromosome where that mark is upregulated. There are multiple differentially active sites on the wild-type chromosomes and no differential sites on the chromothriptic chromosomes in this region. Histone tracks (H3KXXX) show raw haplotype-resolved histone read coverages. In Iso-seq tracks, black lines show splicing, and gray boxes represent exons.