Balanced chromosomal rearrangements offer insights into coding and noncoding genomic features associated with developmental disorders

International aggregation, sequencing, and genome-wide analyses of BCRs in 406 DD cases and 304 controls

We have previously shown that chromosomal rearrangements that appear balanced at cytogenetic resolution can involve extensive complexity ranging from multiple cryptic breakpoints to balanced chromosomal shattering, or chromothripsis, at sequence resolution.^11,15,34 Here, we focused analyses on the most interpretable classes of BCRs by aggregating a cohort of 710 unrelated individuals harboring a “simple” BCR (i.e., breakpoints at two genomic positions without significant imbalance or additional complexity) initially identified by cytogenetic methods and subsequently resolved using either short-insert or long-insert WGS (Fig. 1A; Supplementary Fig. 1; Supplementary Methods). This cohort included 406 cases diagnosed with a DD or congenital anomaly in which a BCR was confirmed to have arisen de novo or segregate with phenotype and 304 unaffected control adults with no early-onset pediatric phenotype (see Fig. 1B and Supplementary Table 1 for complete descriptions). Over 55% of the BCR breakpoints have not been previously published. As a comparison for these 710 empirically-identified BCRs (n = 1,420 breakpoints; Supplementary Table 2), we also generated a set of 30,400 simulated BCRs under the null hypothesis that breakpoints should be randomly distributed throughout the genome. These simulated BCRs were randomly sampled in silico from the genome while matching properties of the 304 BCRs empirically identified in controls, including structural variant (SV) type, inversion size, as well as excluding N-masked regions known to be inaccessible to short read alignments (Supplementary Methods).^35,36

• Download figure
• Open in new tab

Fig. 1 The properties of BCRs in the healthy human germline and DD cases

(A) We mapped the breakpoints of 710 simple (i.e., two-breakpoint) BCRs that were originally detected with cytogenetic methods. Here, we provide the genome-wide BCR breakpoint density in 10Mb windows per chromosome. (B) These 710 BCRs were identified in the genomes of 406 individuals affected by DDs and 304 unaffected controls. For purposes of comparison, we also generated 30,400 synthetic BCRs in silico by resampling the distribution of control BCRs 100 times randomly from the genome. (C) BCR breakpoints were distributed across the chromosomes as expected between affected, control, and simulated subsets except for chromosome 14, which exhibited a significant enrichment of breakpoints in affected samples, and chromosome X, which was depleted of breakpoints in controls. All comparisons were Bonferroni-adjusted for 72 total tests. (D) Balanced translocations involving at least one sex chromosome were significantly depleted in control samples compared to either affected samples or simulated null expectations. (E) We conducted association tests per cytoband for BCR breakpoints in affected samples vs. simulated null expectations. After correcting for 862 independent tests, just one cytoband was significant at the Bonferroni-adjusted threshold of P ≤ 5.8×10⁻⁵): 17q24.3, which corresponds to a locus with well-described LRPEs in DDs involving SOX9 and KCNJ2.^24,25 (F) The distribution of autosomal BCR breakpoints in controls across a “meta-chromosome” arm (i.e., chromosome size-normalized position) was significantly different from those of affected samples or simulated null expectations, with controls exhibiting depletions near centromeres and enrichments very close to telomeres. (G) Loci corresponding to empirically identified autosome-autosome translocation breakpoints (i.e., those sequenced in affected or unaffected genomes) contacted each other in 3D within the nucleus of human embryonic stem cells 1.08-fold more frequently than expected vs. simulated null expectations.

We first sought to understand the global patterns of BCRs throughout the genome by comparing the rates of breakpoints per chromosome between cases, controls, and simulations (Supplementary Fig. 2). We found that BCR frequency was approximately proportional to chromosome length with two exceptions: translocations were enriched on chromosome 14 in DD cases (P = 5.8×10⁻⁶ for cases vs. random simulations) and were depleted on chromosome X in controls (P = 1.4×10⁻⁶ for controls vs. random simulations) (Fig. 1C and Supplementary Fig. 2). The distribution of breakpoints across chromosome 14 did not appear to cluster in any particular location (Supplementary Fig. 3) and we did not find any features (i.e., compartment state, replication timing, recombination frequency, or gene disruption) that could account for the enrichment. In contrast, when we subset chromosome X, we observed that both cases and controls were 2.7-fold to 5.3-fold depleted for breakpoints on the q-arm (Supplementary Fig. 2 and 3). The Xq depletion in controls may be partly explained by the exclusion of males with oligo/ azoospermia and females with premature ovarian failure given their known association with X-autosome translocations.^37,38 However, most of our DD cases were too young to be assessed for infertility, thus the Xq depletion in cases is unlikely to be related to a similar ascertainment bias. Overall, we identified that most (65%; 20/31) translocation breakpoints on chromosome X localized to the p-arm (length=58.6 Mb), which exhibited a 6.9-fold enrichment of case vs. control breakpoints (P = 0.002). Interestingly, 90.0% (18/20) of the Xp translocations identified in DD cases were found in females. Finally, controls were depleted for translocations involving either sex chromosome: just 2.1% of control translocations involved either chromosome X or Y, which was significantly less than translocations in DD cases (8.9%; P = 9.2×10⁻⁴) or randomly simulated BCR carriers (12.1%; P = 5.86×10⁻⁷) (Fig. 1D).

Our assessments of BCR breakpoint distributions per chromosome led to two additional discoveries. First, a single cytoband on chromosome 17 was significantly enriched for BCRs in cases (P = 1.2×10⁻⁵ vs. random simulations) and surpassed a genome-wide significance threshold adjusted for all 862 cytobands tested across all chromosomes (Fig. 1E; Supplementary Fig. 4).³⁹ This cytoband, 17q24.3, matches the location of a well-described pathogenic LRPE in DDs and congenital anomalies caused by SVs altering the local TAD organization and dysregulating SOX9 and KCNJ2.^24,25 Second, after transforming the position of each breakpoint into a percentile relative to the length of its corresponding chromosome arm (i.e., meta-chromosome), we found that translocation breakpoints in controls were biased towards the most distal ends of chromosomes (Kolmogorov-Smirnov test; P = 0.002 for control vs. simulation and P = 0.021 for control vs. cases; Fig. 1F; Supplementary Fig. 5). For example, control BCR breakpoints were roughly three-fold enriched within the terminal 2% of each chromosome arm: 2.7-fold vs. DD cases and 3.3-fold vs. random simulations. This might suggest that translocations occurring near telomeres–which do not rearrange most of the affected chromosome–are more likely to be tolerated in the general population without leading to severe disease.

We also sought to identify genomic features that predispose to BCR formation by annotating all BCR breakpoints with features relating to chromosome maintenance (e.g., recombination rate, replication timing), chromatin accessibility, sequence context (e.g., repetitive elements, sequence homology), and three-dimensional (3D) nuclear organization (e.g., Hi-C contact frequency, nuclear compartment state).⁴⁰ Most features showed no significant differences from expectations after correcting for multiple testing. One feature of note was that our empirically-observed translocations (i.e., those sequenced in DD cases and controls) were slightly more likely to form between pairs of chromosomes in close proximity to each other in 3D within the nucleus than predicted from simulated breakpoints that did not account for this biological organization (1.08-fold increase; P = 0.002; Fig. 1G). This result was true in a fetal lung fibroblast cell line (IMR90) and replicated in a second dataset derived from embryonic stem cells.⁴⁰ These findings might suggest a weak influence on the formation of BCRs between chromosomal regions that co-localize within the nucleus, as suggested by analyses of tumor genomes and cytogenetic data from germline BCR carriers.⁴¹

BCRs in DD cases are strongly enriched for direct disruption of established disease genes

Previous studies have demonstrated that BCRs confer substantial risk for DDs through direct disruption of haploinsufficient, developmentally critical genes.^8,11,15 However, the absence of matched cohorts of unaffected control BCR carriers has historically hindered the quantification of disease risk contributed by gene-disruptive BCRs. Here, we annotated all BCR breakpoints for direct gene disruptions using Gencode v19 and compared the frequency of gene-disrupting autosomal BCRs between DD cases, controls, and random simulations.⁴² Most BCRs disrupted at least one protein-coding gene and there was no difference between cases and controls (68.1% of DD cases and 67.6% of controls) or expectations from random simulations (69.0% expected; Fig. 2). We further subdivided protein-coding genes into four tiers based on the evidence for their association with disease (Supplementary Table 3). Briefly, these included genes associated with dominant DDs (i.e., Tier 1, n = 812), genes associated with all other diseases (Tier 2, n = 3,129), mutationally constrained genes with no prior disease association (Tier 3; n = 1,257), and all remaining protein-coding genes (Tier 4; n = 15,188). We found a strong enrichment of cases with BCR breakpoints directly disrupting Tier 1 genes compared to controls (21.3% of cases versus 3.4% of controls; odds ratio [OR] = 7.45; 95% confidence interval [CI] = 3.74-16.50; Fisher’s exact test; P = 1.60×10⁻¹²) and compared to simulations (21.3% of cases vs. 6.6% of simulations; OR = 3.67; 95% CI = 2.80-4.76; P = 2.06×10⁻¹⁸), but not for Tiers 2-4.

• Download figure
• Open in new tab

Fig. 2 Risk conferred for DDs by BCRs disrupting genes and other transcribed loci

We annotated all BCR breakpoints for predicted overlap with protein-coding and noncoding genes present in Gencode v19.¹⁴ Here, we further subset these BCRs based on the properties of the gene(s) disrupted at either breakpoint, including: any gene present in Gencode; protein-coding genes; “Tier 1” genes including those known to be associated with dominant DDs; “Tier 2” genes including all remaining disease-associated genes; “Tier 3” genes including all genes in the top decile of loss-of-function constraint⁴⁷ but with no existing disease association; “Tier 4” genes including all remaining genes not captured in the preceding tiers; lincRNAs; all noncoding RNAs other than lincRNAs present in Gencode. For each subset of genes, we computed the rate of BCRs disrupting at least one qualifying gene between cases, controls, and simulated BCRs, and further computed the odds ratios of cases vs. controls and cases vs. simulated null expectations. Only BCRs disrupting Tier 1 genes were significantly enriched in cases vs. controls and cases vs. simulated BCRs after correcting for multiple comparisons.

Motivated by the strong association between disruption of dominant DD genes and BCRs in cases, we systematically searched for genes that were recurrently disrupted by BCRs in cases beyond expectations by conducting association tests for each autosomal gene (Supplementary Table 4). These analyses identified four protein-coding genes disrupted in at least three independent DD cases and none in controls (Table 1). Among these, just one gene, TCF4,^11,15 surpassed a strict exome-wide significance threshold (disrupted in 1.6% of DD cases vs. 0.01% of simulated BCRs; P = 4.3×10-10; OR = 149.5). Haploinsufficiency of TCF4 is the dominant genetic cause of Pitt-Hopkins Syndrome and has been associated with autism spectrum disorder (ASD) and broadly defined neurodevelopmental disorders (NDDs).^43,44 The three remaining protein-coding genes did not reach exome-wide significance but had suggestive (P ≤ 0.005) evidence of association with DDs in our analyses. These included two established DD genes (AUTS2, MBD5)^45,46 and one candidate DD gene, CDK6, which was disrupted in three cases that presented with developmental delay (n = 3), speech delay (n = 2), microcephaly (n = 2), and cardiac defects (n = 1). CDK6 is highly constrained against damaging point mutations,⁴⁷ is ubiquitously expressed across tissues,^48,49 and encodes a cyclin-dependent kinase with major roles in skin, blood, and breast cancers.⁵⁰ CDK6 has also been associated with a recessive form of primary microcephaly,⁵¹ but has not been previously associated with dominant germline disease.

Pathogenic positional effects from disruption of three-dimensional chromatin structures

Our analyses identified an unambiguous, strong association between DDs and BCRs disrupting dominant DD genes; however, the majority of autosomal BCRs in DD cases did not disrupt a known disease gene (n = 289; 78.7%), and one-third (n = 113; 30.8%) did not disrupt any annotated protein-coding gene. We therefore considered three other models by which BCRs might confer DD risk through noncoding mechanisms based on: (i) their disruption of noncoding genes (e.g., lincRNAs), (ii) their linear distance from known disease genes, and (iii) the disruption of TADs containing known disease genes. We excluded all 78 (21.3%) cases and 10 (3.4%) controls with autosomal BCRs that directly disrupted a Tier 1 gene, which we reasoned would largely exclude the confounding influence of BCRs associated with pathogenic effects via direct gene disruption. We first tested whether direct disruption of noncoding genes could be responsible for pathogenic effects in DD cases and observed no difference in the fraction of cases versus controls that disrupted any subgroup of noncoding genes (Fig. 2). We next assessed whether pathogenic LRPEs could be predicted based on the absolute distance between disease genes and BCR breakpoints and observed no difference between DD cases and controls for proximity to any tier of genes (e.g., Tier 1 genes; Kolmogorov-Smirnov test; P = 0.159; Fig. 3A-B). These analyses confirm that linear distance to disease genes alone is insufficient to predict pathogenic LRPEs. However, when we tested the third model by comparing the fraction of cases to controls with a BCR breakpoint disrupting a TAD containing genes from each tier, we observed a significant effect for TADs containing Tier 1 genes (OR = 1.43; 95% CI = 1.68-3.18; Fisher’s exact test P = 0.033; Fig. 3C-D), supporting the role of 3D chromatin topology disruption in pathogenic LRPEs. Given these results, we next performed genome-wide analyses to define the TADs most strongly associated with DD phenotypes.

• Download figure
• Open in new tab

Fig. 3 Disruption of TADs, and not proximity to known dominant DD genes, is predictive of LRPEs

(A-B) Fraction of cases and controls with an autosomal BCR breakpoint in proximity to a Tier 1 dominant DD gene (Supplementary Table 3) when direct disruption of Tier 1 genes are included (A) and excluded (B). P value corresponds to a two-sample Kolmogorov-Smirnov test. (C-D) Fraction of cases compared to controls with a BCR breakpoint directly disrupting TADs⁴⁰ containing each of the four gene tiers when direct disruption of genes within that gene tier are included (C) and excluded (D).

We searched for specific TADs associated with risk for DDs by evaluating each autosomal TAD identified from a fetal lung fibroblast (IMR90) cell line⁴⁰ for an enrichment of case BCRs against a Poisson null model fit to the distribution of control breakpoints. Overall, we identified 26 recurrently disrupted TADs with suggestive evidence for association with DDs based on a Benjamini-Hochberg false discovery rate (FDR) q ≤ 0.1 (Fig. 4 and Supplementary Table 5). Five of these TADs surpassed a Bonferroni-adjusted genome-wide significance threshold of 2.2×10⁻⁵, including three known LRPE loci at MEF2C, FOXG1, and SOX9 (Fig. 5).^{11,14,52–54} These five TADs also remained significant when we compared BCRs from DD cases against simulated breakpoints, suggesting that our models were not simply capturing rearrangement hotspots. Consistent with the finding that most TAD boundaries are tissue-invariant,^28,29,40 tissue source had no impact on the genome-wide significant TADs (Supplementary Fig. 6). Given that we only removed cases and controls with a direct disruption of a Tier 1 gene, our TAD results represented a combination of true LRPEs, genic effects not previously associated with DDs, and other unknown mechanisms of disease etiology. For example, one of the genome-wide significant TADs was altered by three cases that all directly disrupted CDK6, a novel candidate DD gene identified from our exome-wide gene association analysis, suggesting that it likely represents a genic effect and not a LRPE.

• Download figure
• Open in new tab

Fig. 4 Genome-wide enrichment of TADs disrupted by BCR breakpoints from DD cases

(A) Genome-wide enrichment of BCR breakpoints in DD cases across 2,257 autosomal TADs.⁴⁰ P values correspond to enrichments against a Poisson null model fit to the distribution of control breakpoints. The genome-wide significance threshold of 2.21×10⁻⁵ (denoted by the red line) was determined by correcting for the total number of autosomal TADs tested and the blue line represents the Benjamini-Hochberg FDR<0.1 cutoff. Known dominant DD genes (Tier 1) contained within each genome-wide significant TAD are reported in parentheses. The GM12878 Hi-C maps³⁰ are shown for each of the six TADs significantly enriched for BCR breakpoints from DD cases: (B) 2p16.1 TAD containing BCL11A, (C) 5q14.3 TAD containing MEF2C, (D) 14q12 TAD containing FOXG1 and PRKD1, (E) 14q32.2 TAD containing BCL11B, (F) 17q24.3 TAD containing SOX9, and (G) 2q33.1 TAD containing SATB2. The annotation tracks under the HiC maps include BCR breakpoints from DD cases and controls, Tier 1 DD genes (blue), all other protein-coding genes from Gencode v19 (gray),⁴² VISTA enhancers (pink),⁶² and UCEs (green).⁶³

• Download figure
• Open in new tab

Fig. 5 Hi-C analysis of DD case with complex BCR disrupting the TAD containing MEF2C

(A) The 3D resequencing pipeline starts by using Hi-C data to call short variants (SNPs and indels) against a haploid reference. In this paper we used the DRAGEN software,⁷⁹ but similar results can be achieved with other publicly available variant callers such as GATK.⁸³ We then use the Hi-C alignment data as generated by Juicer⁸⁴ in conjunction with the 3D-DNA phaser to phase the variants and produce chromosome-length haploblocks. The phased variants enable the generation of molecule-specific contact maps, which in turn allow for molecule-specific annotation of SVs. Using the assembly tools from the 3D-DNA/Juicebox Assembly Tools ecosystem we then create assisted assemblies congruent with the annotated SVs and remap the contact data against the new reference to allow for phased diploid epigenetic analyses. (B) Molecule-specific Hi-C contact maps showing DNA-DNA interactions in the vicinity of the MEF2C promoter in LCLs derived from patient DGAP101: “normal” haplotype Hi-C data mapped to hg19 reference genome (left); haplotype with a chromothriptic chromosome 5 (middle, notice the numerous signal depletions along the diagonal corresponding to breaks and off-diagonal enrichments in the signal corresponding to fusion points); chromothriptic haplotype remapped against a reference that accounts for the chromotriptic rearrangements (right). The 1D tracks show the phased SNP density, highlight the syntenic regions between the three maps (rainbow colors are reserved for sequences in the vicinity of MEF2C in the “normal” reference, while hatching corresponds to sequences juxtaposed into the genomic segment of interest from elsewhere on chromosome 5 in the affected haplotype), as well as show the position of the promoter and 16 known enhancers⁶¹ in the ‘standard’ human reference as well as in the SV-corrected reference. (C) A dotplot of the whole chromosome 5 showing the correspondence between the affected and the normal molecules (100Kb synteny blocks are used, with direct synteny blocks colored red, and inverted blocks colored blue). The position of the MEF2C promoter is highlighted with dashed lines.

To identify candidate pathogenic LRPE loci among the 26 TADs with q ≤ 0.1, we required that each TAD: (i) be disrupted by a noncoding BCR breakpoint (e.g., does not disrupt a protein-coding gene) in at least 50% of cases contributing to the signal, (ii) contain a Tier 1 gene representing a plausible target gene, and (iii) have multiple cases disrupting the same TAD that present with phenotypes that are frequently observed in cases with direct disruption of the candidate target gene (Supplementary Table 6). This resulted in six candidate pathogenic TADs containing the known DD genes SATB2, MEF2C, FOXG1, SOX9, BCL11B, and BCL11A (Fig. 4). Supporting our statistical enrichments, five of the six significant TADs have been previously associated with pathogenic LRPE loci (SATB2, MEF2C, FOXG1, BCL11B and SOX9) based on individual case reports,^{11,14,52,54-56} suggesting that we are accessing bona fide LRPE signals. The novel candidate LRPE at 2p16.1 contains BCL11A, which encodes a zinc finger protein involved in the BAF SWI/SNF chromatin remodeling complex, and has been previously associated with an intellectual disability syndrome.⁵⁷ All four of the cases with noncoding BCRs disrupting the TAD containing BCL11A presented with DD or ASD.⁵⁷ Overall, these six TADs were disrupted by 30 cases (7.4%) and one control (0.32%; Fisher’s exact test; OR = 21.2; 95% CI 4.5-500.1; P = 7.12×10⁻⁷), suggesting that they represent highly-penetrant LRPE loci.

Hi-C analyses of individuals with noncoding BCRs disrupting the TAD containing MEF2C

We previously implicated the 5q14.3 locus as a putative pathogenic LRPE with MEF2C as the target gene based on a statistically significant enrichment of noncoding BCR breakpoints that all disrupted the same TAD containing MEF2C and observed down-regulation of this gene in multiple cases harboring the noncoding BCRs.¹¹ Based on these data, we hypothesized that the disruption of 3D topological organization could represent the underlying mechanism for this LRPE. To functionally validate this hypothesis, we generated high-throughput chromatin conformation capture (Hi-C) data from lymphoblastoid cell lines from five cases harboring BCRs disrupting the TAD containing MEF2C and developed a 3D resequencing workflow (see Supplementary Methods and Fig. 5A) to facilitate analysis of the resulting data. The goal of this workflow was to use Hi-C datasets to: (i) identify single nucleotide polymorphisms (SNPs), small insertions and deletions (indels), and SVs, (ii) phase these variants onto chromosome-length haploblocks, thereby reconstructing the end-to-end sequences of each molecule, and (ii) use the resulting diploid assembly to generate homolog-specific 3D contact maps.

We applied this workflow to two simple (one inversion and one translocation) and three additional complex BCRs from prior studies with noncoding breakpoints disrupting the TAD containing MEF2C.^11,14 Each rearrangement was genotypically distinct with different resultant derivative chromosomes. While complex BCRs were excluded from other aspects of this study, we included three in our Hi-C analysis because their impact on 3D genome organization has not been previously examined in a homolog-specific manner.^59,60 Comparing the results of rearrangement detection using Hi-C to those using liWGS, we found that the 3D resequencing pipeline detected 92.7% (n = 51/55) of the breakpoints found by the combination of both methods. Three of the four breakpoints missed by Hi-C were short segments (<10kb) of DNA that had been rearranged and inserted into a new position. The missed breakpoints were visible in the Hi-C map but had not been identified by the computational analysis. Conversely, liWGS detected 96.4% (n = 53/55) of the Hi-C identified breakpoints, failing to detect a breakpoint associated with a short interval, as well as a 78kb deletion. Crucially, Hi-C was able to reliably order and orient the rearranged sequences on each homolog, even when a breakpoint was missed (Supplementary Table 7, Supplementary Figs. 7A-E). By contrast, it is challenging to reliably order and orient the rearranged sequencing using liWGS data alone if breakpoints are missed. Taken together, these results demonstrate that Hi-C can be used to robustly generate both homolog-specific sequences and architectural maps.

In all five BCR cases we examined, 3D resequencing via Hi-C also revealed significantly altered 3D organization of the rearranged homolog (Supplementary Figs. 7A-E). Moreover, we observed reduced frequency of contact between the MEF2C promoter and 16 distal enhancers⁶⁰ (Supplementary Table 8), consistent with the dysregulation of MEF2C expression observed in the same cases. In the majority of the cases, the reduction in contact frequency appeared to result from the BCR greatly increasing the distance between the promoter and its enhancers (Supplementary Fig. 7). However, in one case with chromothripsis (DGAP101), the BCR breakpoints only had a modest effect on the linear distance between the MEF2C promoter and its enhancers. Instead, the promoter and enhancers were separated into distinct architectural domains through the creation of a new boundary, which likely prevented physical contact between the promoter and enhancers (Fig. 5B-C). The observed 3D remodeling suggests a reduced frequency of contact between sequences that influence MEF2C expression, providing a plausible explanation for how a noncoding BCR breakpoint can result in a DD phenotype through disruption of 3D genome architecture.

Genomic features predict TADs associated with pathogenic LRPEs

Motivated by our discovery of multiple genome-wide significant LRPE loci in DDs, coupled with our validation of functional changes in MEF2C and alterations to the 3D organization associated with the putative 5q14.3 LRPE, we exploited this unique BCR dataset to identify additional features contributing to pathogenic TAD disruptions that could be used for future LRPE predictions. As described above, we demonstrated that the genic content of TADs alone is insufficient to unequivocally predict the pathogenicity of an individual BCR, as 43.0% of BCRs in controls disrupted a TAD containing a Tier 1 gene. We also found that 9.4% of all TAD boundaries encompassing a Tier 1 gene were overlapped by at least one large polymorphic deletion in the genome aggregation database (gnomAD),⁶¹ which excludes adults with a history of early onset developmental conditions. Thus, additional genomic features beyond the presence of disease-associated genes are required to predict TADs preferentially disrupted by BCRs in DD cases.

To identify genomic features that characterize TADs intolerant to disruption, we annotated all autosomal TADs⁴⁰ in the genome with 54 features that can be broadly grouped into five categories: genes, cis-regulatory elements, primary sequence conservation, repetitive elements, and ‘other’ (Figs. 6A and Supplementary Methods). We defined 45 “positive” training TADs (disrupted by ≥2 BCR cases and zero BCR controls) and 261 “negative” TADs (disrupted by ≥1 BCR control and no BCR cases) and performed a univariate logistic regression for each of the 54 features, which identified 26 features at a FDR<0.05 (Supplementary Table 9). Next, given that many genomic features are highly correlated (Supplementary Fig. 8), we trained an elastic net regression on the positive and negative training TADs that included all 26 features from the univariate analysis and identified six features that were individually associated with case status after controlling for the effects of all other features: VISTA enhancers,⁶² ultraconserved elements (UCEs),⁶³ transposon-free elements,⁶⁴ TAD size,⁴⁰ the presence of at least one Tier 1 gene, and primary sequence conservation^65,66 (Fig. 6B).

• Download figure
• Open in new tab

Fig. 6 Identification of genomic features associated with pathogenic LRPEs

(A) An overview of the framework used to generate a model to predict pathogenic LRPEs across the genome based on disruption of 3D genome architecture. (B) Features identified to be potentially associated with TADs preferentially disrupted by BCRs from DD cases based on a BH-FDR<0.05 from a univariate logistic regression performed for each feature, ordered by effect size. (C) Evaluation of the LRPE model using an independent set of 372 TADs. Performance of TADs overlapping other common functional annotations are shown as a comparison. (D) To determine a cut-off score that could be used to identify TADs that are especially intolerant to disruption, we compared the fraction of cases to controls and simulated breakpoints that disrupt a TAD in each percentile and identified the inflection point at which a case enrichment (OR>1.5) begins to emerge.

We tested this model’s predictive accuracy on an independent set of 372 TADs (see Supplementary Methods for selection criteria) that were not used in our training data and determined that these six features alone were moderately predictive of LRPE pathogenicity (area under the receiver operating characteristics curve=0.633; Fig. 6C). We next defined a “LRPE pathogenicity” cutoff score of ≥0.43 (TADs ranked in the top 10th percentile from our model) based on the point in which case-control and case-simulation enrichments surpassed OR>1.5 (Fig. 6D). We demonstrated that 75.0% of previously established LRPE loci (Supplementary Table 10) in the human genome that were not represented in our training dataset surpassed this score, a highly significant enrichment when compared to all other TADs in the genome (Fisher’s exact test; P = 2.55×10⁻⁵). These orthogonal approaches collectively confirmed that this relatively simple six-feature regression model was able to prioritize TADs likely intolerant to disruption. While much larger cohorts will be required to power more sophisticated predictive models, these analyses demonstrate the potential of this approach to improve noncoding variant interpretation and shed greater light on the genomic features associated with noncoding mechanisms of disease.