3D genome of multiple myeloma reveals spatial genome disorganization associated with copy number variations

Abstract

The Hi-C method is widely used to study the functional roles of the three-dimensional (3D) architecture of genomes. Here, we integrate Hi-C, whole-genome sequencing (WGS) and RNA-seq to study the 3D genome architecture of multiple myeloma (MM) and how it associates with genomic variation and gene expression. Our results show that Hi-C interaction matrices are biased by copy number variations (CNVs) and can be used to detect CNVs. Also, combining Hi-C and WGS data can improve the detection of translocations. We find that CNV breakpoints significantly overlap with topologically associating domain (TAD) boundaries. Compared to normal B cells, the numbers of TADs increases by 25% in MM, the average size of TADs is smaller, and about 20% of genomic regions switch their chromatin A/B compartment types. In summary, we report a 3D genome interaction map of aneuploid MM cells and reveal the relationship among CNVs, translocations, 3D genome reorganization, and gene expression regulation.

Fig. 1

Correcting copy number variation in cancer Hi-C data. a Whole genome Hi-C interaction matrix of RPMI-8226 cells. Black rectangles indicate inter-chromosomal translocations. b CNVs of RPMI-8226 cells estimated from different sequencing data sets. Methods for filtered Hi-C data are described in the Material and Methods. c Scatterplot between the copy numbers estimated from WGS (x-axis) and filtered Hi-C data (y-axis) shows high similarity. d Median raw Hi-C interaction counts of the diagonal bins in each CNV block vs. the copy number of CNV blocks. A CNV block is defined as a continuous chromosome region with the same estimated copy number from WGS data. e, f Same as d, but the y-axis data are HiCNorm e and ICE-normalized f interaction values. The p-values of figures d–f were determined by F-test

Fig. 2

Hi-C data reveal translocation events. a Fifty-six inter-chromosomal translocation events identified in RPMI-8226 by WGS. From outer circle to inner circle: panel I: average RNA-seq count for every 200 kb chromosome bin; panel II: CNV data from WGS (black dots: copy number equal to 2, red dots: copy numbers larger than 2, green dots: copy numbers less than 2); panel III: orange lines: translocation events, blue lines: translocation events that are among the top 100 highest inter-chromosomal Hi-C interactions in b. b The 100 highest inter-chromosomal Hi-C interactions in RPMI-8226 cells (5 M resolution): from outer circle to inner circle: panel I and II: the same as in a for RNA-seq counts and CNV data; panel III: orange lines: top 100 highest interactions, blue lines: the common events with a. c The chimeric chromosome formed between chr16 and chr22 identified by spectral karyotyping (SKY). The SKY image of RPMI-8226 is from the NCBI SKY database. d The translocation event of t(16;22)(q23.2;q11.22) in RPMI-8226 cells plotted by IGV tool, showing only WGS paired reads supporting this translocation event. The reads depth and alignment panels are shown. Read alignments are sorted by read location and colored by read strands. e Top panel: the intra and inter-chromosome interaction maps of chromosome 16 and 22 in RPMI-8226. Bottom panel: copy number variations of chromosome 16 and 22. f Enlarged inter-chromosomal interactions corresponding to the translocation t(16;22)(q23.2;q11.22)

Fig. 3

CNV breakpoints associated with TAD boundaries. a Hi-C interaction matrix of a region (chr15: 80 Mb-92 Mb) in RPMI-8226 shows the association between TAD boundaries and CNV breakpoints. Top: Hi-C interaction matrix, middle: TAD boundaries (vertical bars) and insulations scores, bottom: copy number variations. The dashed triangles indicate two TADs coinciding with CNV blocks. b The distribution plots of the distance between all CNV breakpoints (red line) or random sites (blue line) to their nearest TAD boundaries in RPMI-8226 cells. c The averaged insulation scores at TAD boundaries (black line) or CNV breakpoints (red line) in RPMI-8226 are lower compared to surrounding regions, in contrast to random sites (blue line). d The averaged insulation scores around different sets of genomic locations in B cells (GM12878). Blue line: randomly selected sites; black line: TAD boundaries in GM12878; red line: CNV breakpoints from the CNVD database

Fig. 4

Reorganization of 3D chromosome is associated with gene expression changes. a The A/B compartments in RPMI-8226 and U266 compared with normal B cells (GM12878) in chromosome 1. b Genome-wide proportions of A/B compartment changes among RPMI-8226, U266 and GM12878 cells. c Boxplots of expression changes of genes grouped by their A/B compartment changes (t-test). d Example of conserved and changed TADs in a region (chr14: 24 Mb-30 Mb) comparing RPMI-8226, U266 and GM12878 cells. e The average size of TADs in multiple myeloma cells compared with GM12878. f The number of conserved and changed TADs in multiple myeloma cells compared with GM12878

Fig. 5

The 3D structure and expression changes of multiple myeloma-related genes. a Gene pathway enrichment analysis for the switched on/off genes. Switched on genes are defined when their genomic compartment changed from B to A in both MM cell lines (compared with GM12878). Switched off genes are defined when their genomic compartment changed from A to B in both MM cell lines (compared with GM12878). Genes with mutations that affect protein coding were discarded. b Example of cytokine receptor gene cluster including IL1R1, IL1R2, IL18R1. The IL1R2 gene is located at chr2:101991844-102028544, and the plotted locus is from 1 Mb upstream to 1 Mb down-stream region surrounding the IL1R2 gene. Red color gene names indicate the genes that are enriched in the cytokine-cytokine receptor interaction pathway. In GM12878, IL1R1, IL1R2, IL18R1 and MAP4K4 were highly expressed. Insulation score, TAD boundaries (vertical red lines), copy numbers, ChIP-seq peaks of H3K27ac, H3K4me1, H3K4me3, H3K36me3, H3K27me3, H3K9me3, RNA-seq expressions and RefSeq gene tracks are shown below the Hi-C heatmaps. H3K4me3 is annotated as promoters, H3K4me1/H3K27ac as active enhancers, H3K36me3 as transcribed gene body and H3K9me3 as heterochromatin. c The same genomic region as in b using the data of U266 cells. The expression of IL1R1, IL1R2, IL18R1 and MAP4K4 are downregulated and the TAD boundaries are altered, accompanied by reduced epigenetic marks such as H3K27ac, H3Kme1 and H3K4me3