Skip to main content
Download PDF

Comprehensive network modeling approaches unravel dynamic enhancer-promoter interactions across neural differentiation

Genome Biology volume 25, Article number: 221 (2024) Cite this article

Abstract

Background

Increasing evidence suggests that a substantial proportion of disease-associated mutations occur in enhancers, regions of non-coding DNA essential to gene regulation. Understanding the structures and mechanisms of the regulatory programs this variation affects can shed light on the apparatuses of human diseases.

Results

We collect epigenetic and gene expression datasets from seven early time points during neural differentiation. Focusing on this model system, we construct networks of enhancer-promoter interactions, each at an individual stage of neural induction. These networks serve as the base for a rich series of analyses, through which we demonstrate their temporal dynamics and enrichment for various disease-associated variants. We apply the Girvan-Newman clustering algorithm to these networks to reveal biologically relevant substructures of regulation. Additionally, we demonstrate methods to validate predicted enhancer-promoter interactions using transcription factor overexpression and massively parallel reporter assays.

Conclusions

Our findings suggest a generalizable framework for exploring gene regulatory programs and their dynamics across developmental processes; this includes a comprehensive approach to studying the effects of disease-associated variation on transcriptional networks. The techniques applied to our networks have been published alongside our findings as a computational tool, E-P-INAnalyzer. Our procedure can be utilized across different cellular contexts and disorders.

Background

Enhancers are cis-regulatory elements (CREs) essential to gene regulation. Ranging between ten to 10,000 base pairs (bps) in length [1], these CREs contain motifs, clustered recognition sites for transcription factors (TFs) [2]. Enhancers have proven arduous to locate despite their abundance in the non-coding genome. Enhancers interact with promoters to mediate their effects on transcription, often by regulating non-proximal genes from variable distances [3]; this contrasts with promoters, CREs that initiate transcription from known locations [4]. Enhancers are highly cell-type-specific regulatory elements, which adds another barrier to their identification and validation [5]. There is an emerging consensus that mutations associated with human disease in genome-wide association studies (GWAS) are concentrated in these enhancer regions [6], highlighting their importance to understanding gene regulation and their applications in medicine.

Previously, we collected high-throughput assay for transposase-accessible chromatin (ATAC-seq), chromatin immunoprecipitation (ChIP-seq) for H3K27ac, massively parallel reporter assay (MPRA), and RNA-seq datasets from seven time points during the differentiation of human embryonic stem cells (hESCs) into neural progenitor cells (NPCs) [7]. This was within a 72-h window, with samples collected at hours 0, 3, 6, 12, 24, 48, and 72. Using these datasets, we assembled seven time point-specific enhancer-promoter interaction networks (E-P-INs). ATAC-seq, indicating open chromatin regions, ChIP-seq for H3K27ac, acting as a marker for active enhancers, and RNA-seq, representing the transcriptome, combined with generalized Hi-C datasets, serve as the essential components for these networks. Input into the Activity-By-Contact (ABC) model, these datasets integrate to build networks of predicted enhancers and target promoters [8]. After constructing these networks, information about the enhancers within them was extracted using a series of novel methods: differential expression revealed time point specificity, Girvan–Newman clustering demonstrated the substructures through which enhancers regulate genes, variant analyses revealed enrichment for neurodevelopmental and neuropsychiatric disorders, and TF overexpression and MPRA datasets were used to validate our networks.

Overall, our analysis of these E-P-INs provides a framework for exploring gene regulatory programs and their dynamics across developmental processes, including a comprehensive approach for mapping disease-associated genetic variants to transcriptional networks.

Results

Reconstructing enhancer-promoter interactions with multi-omics datasets

Examining shifts in the intricate mechanisms of regulatory programs across seven stages of neural differentiation necessitated designing a systematic procedure to generate E-P-INs. Our E-P-INs are highly modular bipartite graphs populated with communities of enhancers and promoters; edges represent interactions between these elements. These networks were designed for cross-condition comparability, referencing the same atlas of enhancers and promoters (see the “Methods” section). Our E-P-INs capture the temporal programs of gene regulation within our seven time points exclusively. Ubiquitous CREs that are unchanging across time points were removed during filtration.

Eight E-P-INs were created using the ABC model: seven time point-specific networks and a single collapsed E-P-IN containing all edges found across these time points. Here, we report statistics from the collapsed network to summarize the structures of the seven time point-specific E-P-INs. Our E-P-INs contained 3,133 unique promoters and 1628 enhancers. This suggests that most interactions in these E-P-INs, consisting of differentially active and temporal CREs, involve multiple genes regulated through a single enhancer, with an average of 1.92 genes regulated through a single enhancer. The mean length of enhancers in our E-P-INs was 1047.83 bps; this aligns with previous literature [1, 9]. Examining sample composition, we found consistent interactions within our time point-specific networks. The smallest E-P-IN (48 h) contained 2134 interactions, totaling 13.02% of our collapsed network’s edges. The largest E-P-IN (24 h) included 15.68% of edges in the collapsed network, amounting to 2569 unique enhancer-promoter links. Additional file 2: Tab. S1 in “Supplementary information” contains our time point-specific and collapsed E-P-INs.

Enrichment in neural differentiation signals

Initial investigations of our E-P-INs showed characteristics relevant to early neural differentiation through enrichment and temporal analyses. Our E-P-INs demonstrated enrichment for biological processes pertinent to the cellular context and temporal changes from the start of neural induction to the 72-h time point.

Gene Set Enrichment Analysis (GSEA) was performed utilizing the 500 most frequently occurring promoters from each E-P-IN [10]. The Genomic Regions Enrichment Annotations Tool (GREAT) was used to examine enrichment in our enhancer regions [11, 12]. See the “Methods” section for an in-depth explanation of GSEA and GREAT’s methodologies. We report the results from our collapsed E-P-IN for brevity; time point-specific results are available in Additional file 3: Tab. S2 and Additional file 4: Tab. S3 in “Supplementary information”. The most significant results from GSEA and GREAT are reported in Fig. 1A and B, respectively. “GOMF_TRANSCRIPTIONAL_REGULATOR_ACTIVITY,” “primary neural tube formation,” and “ventral spinal cord interneuron specification” are examples relevant to neural differentiation. We note that GSEA results relate to transcription due to the number of promoters being tested and their recurrence across the networks. GREAT considers all enhancers in our E-P-INs. A concentrated GSEA testing set contributes to the difference in statistical magnitudes between GSEA and GREAT. These findings support our E-P-INs’ biological relevance.

Fig. 1
figure 1

Enrichment analyses demonstrate the E-P-INs’ abilities to capture neural differentiation. A GSEA outcomes from the collapsed network’s 500 most occurring promoters reveal biological processes, cellular components, and molecular functions essential to transcriptional regulation. B GREAT outcomes detailing enhancer enrichment show biological processes related to cell differentiation and neural development

Full size image

Next, Jaccard similarities between the CREs (i.e., enhancers and promoters) and enhancer-promoter interactions (E-P-Is) were computed to compare our seven time point-specific E-P-INs. Figures 2A, B, and C present Jaccard similarities in enhancers, promoters, and E-P-Is, respectively. Here, we show that these networks become more dissimilar with time. Additionally, this analysis revealed that time points before 48 h show high similarity compared to our 48-h and 72-h E-P-INs; these latter networks also show close resemblance. These results indicate our E-P-INs exhibit temporal trends of neural differentiation.

Fig. 2
figure 2

Temporal analysis revealed the E-P-INs exhibit changes related to neural differentiation. A Jaccard similarities detail the overlap between sets of promoters across our seven time point-specific networks, with clustered formation in the earliest five time points. B Jaccard similarities detailing differences in time point-specific sets of enhancers reveal temporal trends similar to those described for promoters. 2C) Jaccard similarities demonstrating changes in the E-P-Is across neural differentiation. Instances where an enhancer-promoter interaction changes across E-P-INs show more temporality than examining the presence/absence of either enhancers or promoters deemed dynamic

Full size image

Clustering reveals regulatory substructures

We identified four distinct substructures defining the modes through which enhancers regulate genes. Here, Girvan-Newman clustering was utilized to partition our bipartite graph into distinct communities [13]. The highly modular nature of our E-P-INs informed this approach, as enhancers and promoters tied together through gene regulation formed dense clusters in our networks. Figure 3 provides a graphical representation of our collapsed E-P-IN to illustrate the formations of these substructures. Substructures fall into two distinct classifications: redundant (R) and non-redundant (NR). The four substructures are a single enhancer regulating a single gene (1NR; Fig. 3A), multiple enhancers regulating a single gene (1R, Fig. 3B), a single enhancer regulating multiple genes (2NR; Fig. 3C), and multiple enhancers regulating multiple genes (2R; Fig. 3D). Additional file 1: Fig. S1 in “Supplementary information” contains a full view of our collapsed E-P-IN [14].

Fig. 3
figure 3

Graphical representations demonstrate the composition of gene regulatory substructures. A offers an example of a 1NR substructure. B shows a community of E-P-Is characterized as 1R. C demonstrates the 2NR substructure. D depicts a cluster labeled as 2R

Full size image

We implemented this clustering approach for each of the seven time point-specific E-P-INs and the collapsed E-P-IN. Figure 4A shows the four distinct substructures as a legend, and Fig. 4B displays the substructure distribution in the seven time point-specific E-P-INs, calculated by dividing the total number of communities by the number of communities allocated to each specific substructure. This distribution remains consistent, with slight fluctuations, across our seven time points. The 1NR and 2NR substructures constitute the majority of our time point-specific networks, averaging 81.9% composition across the seven E-P-INs. This skew towards NR substructures stems from the dynamic and temporal filtering used to reconstruct the networks.

Fig. 4
figure 4

Substructures across the seven time point-specific E-P-INs demonstrate minimal temporality. A is a legend that links substructures to their appropriate section on the chart. B shows the distribution of substructures across the seven time point-specific and collapsed E-P-INs. Redundant substructures comprise a smaller portion of these networks than those non-redundant

Full size image

We calculated average counts of enhancers, promoters, and E-P-Is in these regulatory substructures using our collapsed E-P-IN. In 1R communities, an average of 2.16 enhancers regulate a single gene; in 2NR communities, 3.29 genes, on average, are regulated through a single enhancer. Interestingly, the 2R substructure is denser, averaging 2.55 enhancers, 4.57 genes, and 7.72 E-P-Is per community. Additional file 1: Fig. S2 compiles histograms showing the distribution of these values. Our investigations could not reveal significant links between gene expression and a promoter’s number of edges with enhancers.

Enrichment for disease-associated variants

We hypothesized that our time point-specific E-P-INs would exhibit enrichments associated with neurodevelopmental and neuropsychiatric disorders [15]. Utilizing datasets containing 250,000 de novo autism spectrum disorder (ASD) single nucleotide variants (SNVs) from both patients (cases) and healthy siblings (controls) [16], 817 ASD-associated genes, 252 developmental delay-associated genes (DD) [17], and 34,914 GWAS-sourced SNVs from neurological/neuropsychiatric disorders (i.e., schizophrenia (SCZ), bipolar disorder (BD), major depressive disorder (MDD), Parkinson’s disease (PD), and Alzheimer’s disease (AD)), we implemented a permutation test to quantify an empirical p-value signifying enrichment; our null distribution in these permutation tests was the ABC model’s unfiltered outputs. Eight tests were employed on each E-P-IN: ASD-associated SNVs, ASD-associated genes, DD-associated genes, SCZ-associated SNVs [18], BD-associated SNVs [18], MDD-associated SNVs [19], PD-associated SNVs [20], and AD-associated SNVs [21].

With the significance threshold set at an empirical p-value = 0.05 and 1000 permutations ran, five of our eight categories demonstrated significant enrichment in our E-P-INs: ASD-associated genes and SNVs, DD-associated genes, and SCZ-associated and BD-associated SNVs. Additional file 5: Tab. S4 in “Supplementary information” contains a complete overview of the statistics for disease-associated variant enrichment. MDD-associated SNVs, PD-associated SNVs, and AD-associated SNVs were deemed non-significant. These findings suggest our E-P-INs capture enhancers relevant to neural differentiation; SCZ and BD are disorders with neurodevelopmental origins [22, 23]. Out of the 1592 enhancers considered in our collapsed E-P-IN, 140 contain either a case or control de novo variant for ASD, with case variants favored; 27 contain a GWAS-sourced SNV for a neurological/neuropsychiatric disorder. Of 3133 genes, 208 are associated with ASD, and 41 are associated with DD. ASD was the most significantly enriched among diseases/disorders, with a p-value < 0.001 in examinations of associated SNVs and genes. Figure 5A depicts the results of the ASD-associated de novo SNV permutation test, while Fig. 5B presents the ASD-associated gene permutation test results. Further confirming the enrichment of ASD, we devised an additional permutation test, utilizing the dynamic and temporal E-P-Is as the null distribution. We found that ASD-associated genes demonstrate significant interaction with enhancers containing an ASD-associated SNV, computed at an empirical p-value < 0.001. ASD enrichment arises from a common set of genes and SNVs across our seven time point-specific E-P-INs. Additional file 1: Fig. S3 in “Supplementary information” shows Jaccard similarities between sets of ASD-associated CREs in our networks.

Fig. 5
figure 5

ASD enrichment across E-P-INs is substantial. A demonstrates significant enrichment for ASD-associated de novo SNVs in the collapsed network. B details the enrichment of ASD-associated genes, showing permutation test outcomes. C displays the number of communities within each substructure containing ASD-associated variation. The communities with variants are shaded red at the start of each substructure’s bar. D provides an example of a predicted enhancer region, binding ASD-associated TFs CHD2 and ZNF219, that regulated CHD8, SUPT16H, and RAB2B. Clinical reports provide cases where a mutation in chromosomal region 14q11.2, the location of this enhancer, involved these genes and caused ASD

Full size image

We mapped our ASD-associated de novo SNVs to the communities created with the Girvan-Newman clustering algorithm [13]. The distribution of SNVs was consistent with that of substructures in our E-P-INs; 8.42% of 1NR communities contained variation, 10.34% of 1R communities, 10.15% of the 2NR communities, and 15.75% of 2R communities. Figure 5C details the number of communities in each substructure that contain ASD-associated variants. Prior research has suggested enhancer redundancy as a counteracting mechanism to variation that would disrupt gene transcription [24]. Figure 5D gives an example of the susceptibility of non-redundant substructures, illustrating a case where an ASD-associated SNV fell within a predicted enhancer (chr14:21,993,966–21,995,169). Both TFs, ZNF219 [25] and CHD2 [26], have been previously found to be associated with ASD in transcriptional regulatory networks of neural differentiation through investigation of SNVs affecting the development of the disorder. ZNF219 and CHD2 were selected from a group of motif-predicted TFs. The enhancer regulates three genes (i.e., CHD8, SUPT16H, and RAB2B); these genes have been linked to ASD through differential expression, clinical variant analyses, and gene knockouts, respectively [27,28,29]. Mutations in chromosomal region 14q11.2, where the predicted enhancer is located, have been shown to involve these genes and cause ASD [30, 31].

Validating enhancers with transcription factor overexpression

Integration of TF data into our seven time point-specific networks enabled us to begin verification of E-P-Is and gain a more in-depth understanding of these enhancers’ role in gene regulation [32].

First, we used Find Individual Motif Occurrences (FIMO) [33] with motifs from ENCODE [34] and JASPAR [35] databases to find TFs predicted to bind to the sequences of our enhancers. Our p-value threshold for FIMO was 1e-5. Previous studies suggest that tens of TFs can bind to an enhancer [1]; this was in line with our findings. Enhancers in our E-P-INs averaged 12.30 unique TFs bound to them. Examining our original non-dynamic network, produced with the ABC model prior to the robust filtration process reported in “Reconstructing Enhancer-Promoter Interactions with Multi-omics Datasets” used to generate E-P-INs, we found a decrease in the number of unique TFs per enhancer (7.33). In total, 1183 unique TFs were predicted in our E-P-INs using FIMO. Next, we examined which TFs were more commonly predicted to bind to our enhancers; returning to our ABC network prior to filtration, CTCF was predicted to bind to the largest number of enhancers. Differing in our dynamic E-P-INs, SP9, a TF necessary for neuronal development [36], was the most common. SP9 binds to 1222 of the 1628 enhancers in our time point-specific networks.

Previously, we overexpressed TFs (i.e., BARHL1, IRX3, LHX5, OTX1, OTX2, and PAX6) [7] followed by RNA-seq in later-stage hESC-derived NPCs to detect their target genes. Using DESeq2 for differential expression analysis [37], we determined which genes were significantly affected following a TF’s overexpression. Using these results, we investigated whether enhancers predicted to bind these TFs regulated the differentially expressed genes (DEGs) downstream; downregulated DEGs were not considered in this analysis (see the “Methods” section). Using Fisher’s exact test, p-values were assigned to each TF-binding enhancer, quantifying their relation to the overexpressed TF’s regulon. We found enhancers predicted to bind each of the six overexpressed TFs where a substantial proportion of downstream promoters regulate DEGs. Here, the false discovery rate (FDR) values of the most significant enhancers are reported per TF: BARHL1 = 4.2e − 8, IRX3 = 1.8e − 4, LHX5 = 1.1e − 3, OTX1 = 4.9e − 2, OTX2 = 2.7e − 2, PAX6 = 1.0e − 6. A statically significant presence of downstream DEGs supports motif predictions from FIMO in these enhancers. In addition to this analysis, we designed a permutation test to determine if the proportion of downstream DEGs in enhancers predicted to bind overexpressed TFs was different than in enhancers that are not. Upregulated and downregulated DEGs were examined together. We found that all TFs scored an empirical p-value < 0.001.

PAX6 overexpression is sufficient to induce hESCs into neural lineage [38]. We report an additional enhancer (chr16:30,669,227–30670310) predicted to bind PAX6 here. This enhancer’s targets in our E-P-IN are enriched with DEGs downregulated during PAX6 overexpression. Only PHKG2, one of nine targets for this enhancer, was upregulated during PAX6 overexpression. In the PAX6-binding enhancer found to have a significant proportion of upregulated DEGs downstream (chr19:17438298–17439337), DDA1 was the sole downregulated target DEG. This suggests a pattern of enhancers regulating directional gene sets and a possible method of examining TF antagonism [39]. Figure 6A details the permutation test outcomes from PAX6. Figure 6B illustrates the most significant TF-binding enhancer predicted to bind PAX6 when examining upregulated DEGs.

Fig. 6
figure 6

Enhancers predicted to bind PAX6 more commonly regulated DEGs when PAX6 is overexpressed. A shows outcomes from permutation testing supporting the claim that enhancers binding PAX6 regulate DEGs more often than enhancers without PAX6 motifs. B details an enhancer region predicted to bind PAX6 that regulates upregulated DEGs. Genes shown in lighter colors are DEGs when PAX6 is overexpressed

Full size image

Validating enhancers with massively parallel reporter assays

In MPRAs, a synthetic construct, each comprising a candidate regulatory DNA sequence, a minimal promoter, and a unique transcribable DNA barcode, is introduced into cells; these candidates are assumed to be able to regulate the transcription of these barcodes. Following this insert, RNA and DNA sequencing are performed to evaluate each barcode’s RNA-to-DNA ratio [40]. MPRAs allow us to measure the activity of thousands of regulatory sequences at once. We used our previously collected perturbation MPRA datasets [7, 41, 42] to validate our predicted enhancers. Perturbation MPRA is a strategy for demonstrating the regulatory effects of TF-binding motifs. Here, nucleotides in a motif are randomly scrambled. We then analyze the change in regulatory activity of the DNA sequence being tested to see the effects of TF binding [42]. As perturbation MPRA facilitates the measurement of changes in transcriptional activity based on a sequence’s inclusion of a TFBS, we choose to use this dataset for this analysis. We examined E-P-IN enhancers for strict MPRA overlap, creating a subsetted network. This consolidated E-P-IN allowed for an in-depth investigation into TFs, as we computed a novel metric for TF × TF correlation using enhancer binding, genes regulated downstream, and time point-specific function.

Our subsetted MPRA-based E-P-IN had strict inclusion criteria for edges. A functional MPRA-tested sequence had to intersect each included enhancer. Additionally, these sequences also needed to perturb a TF that matched a FIMO-predicted motif for that enhancer (see the “Methods” section). This E-P-IN contained 45 enhancers and 82 promoters, a ratio comparable to the collapsed E-P-IN (1.82 to 1.92). Applicable edges matching the E-P-INs’ criteria were found more often at some time points. A total of 195 edges in the 6-h network matched, while only 69 in the 72-h E-P-IN were found.

We devised a permutation test to investigate whether MPRA-validated enhancers were significantly enriched in our temporal, dynamic networks (see the “Methods” section). Using the ABC model networks as the null distribution, we computed an empirical p-value < 0.001 from 1000 permutations, showing that proven enhancers are more common in networks of E-P-Is after our filtration methods are utilized.

We devised a TF × TF correlation metric, incorporating enhancers that these TFs are predicted to bind, genes they regulate downstream, and the time points at which the TFs are functional. The metric is based on Jaccard similarities; scores between 0.0 and 1.0 are created for the three components that were previously mentioned. Additional file 6: Tab. S5 contains matrices for these component-specific scores and our TF × TF metric. Figure 7 shows a heatmap of similarity scores subjected to hierarchical clustering to reveal groups of correlated TFs. TF families, such as the NKX-homeodomain TFs [43] and IRF TFs [44], cluster together, supporting our approach’s ability to reveal TFs working in tandem. Analyses were performed on our MPRA-subset networks to provide interpretable results, detailing a group of enhancers binding TFs validated using experimental methods.

Fig. 7
figure 7

TF × TF matrix, stylized as a heatmap and organized using hierarchical clustering, displays our Jaccard-based similarity metric for transcription factors in our subsetted E-PN derived from MPRA datasets

Full size image

Discussion

Understanding gene expression and the architecture of E-P-Is has the potential to transform human health. Three distinct challenges confront researchers attempting to fill knowledge gaps regarding enhancers: locating enhancers, finding their target genes, and contextualizing their functions [45]. We present a comprehensive computational approach to navigate these questions.

The ABC model has proven to be a well-performing method for predicting E-P-Is, with standard use across the field. Additional layering with perturbation MPRAs allows for built-in motif and variant validation [41, 46]. Unraveling the TFs that bind enhancers is pivotal to their identification, as TF binding is linked to enhancer activity [1]. Developing models for E-P-I predictions using MPRA data to define candidate enhancers and score their activity could tremendously increase the strength of the E-P-I predictions, allowing for better capture of cell-type-specific elements. Our E-P-INs subset with MPRA data showed TF per enhancer counts, average bp lengths, and gene-to-enhancer counts in line with previous scientific findings. MPRA, paired with TF overexpression data and differential expression methods, offers a robust framework for the validation of E-P-Is, helping us understand concrete links between genes and these non-coding regions.

Contextualizing enhancer function is tied to temporal dynamics, shifting connections between enhancers and genes that remain dormant during specific stages in development. Using E-P-INs, as we proposed, allows comparison of CREs between different stages of development to be done with ease. Variant information can also be mapped directly into these E-P-INs and their substructures, allowing us to mold disease networks from model systems. Assessing the deleteriousness of the SNVs mapped to R substructures compared to those in NR could elucidate enhancer redundancy’s role in disease and development. This framework is generalizable, and with sufficient datasets, our analysis can be repeated across biological systems.

E-P-INAnalyzer has been developed alongside our analysis as a computational tool that allows other scientific researchers to model biological networks, perform clustering to reveal substructures, integrate TFs and variants, and search for network enrichment [47, 48]. This tool has been made available on our GitHub, found in the “Availability of data and materials” section. With the same ATAC-Seq and ChIP-Seq for H3K27ac needed to run the ABC Model and RNA-seq, E-P-INAnalyzer performs these analyses, returning condition-specific networks. TF and variant datasets are optional for building networks.

Conclusion

In this study, we have provided a robust, generalizable framework for comparing E-P-INs across multiple conditions and cellular contexts. We have developed techniques to verify computationally predicted enhancers using TF overexpression and MPRA. We have also demonstrated the use of these E-P-INs in understanding the mechanisms by which disease-associated variation perturbs regulatory networks. E-P-INAnalyzer, our computational tool, allows researchers to answer similar questions about their model system of interest; this tool creates E-P-INs and provides methods to understand them.

Methods

Reconstructing enhancer-promoter interactions with multi-omics datasets

Our E-P-INs’ architecture was composed using the Python [49] and R [50] programming languages. These languages enabled us to create comparable networks.

Assembling an E-P-IN is contingent on the ABC model [8], which uses epigenetic datasets to define candidate enhancers and predict their interactions with promoters. The ABC model relies on these straightforward biological notions: enhancers are distinguished through signals denoting potential for transcriptional regulation (ATAC-seq, ChIP-seq), and frequent promoter contact (Hi-C) indicates an edge. Here, we ran the ABC model using the authors’ recommended parameters. Peak-calling was performed using MACS2 on ATAC-seq [51]; these peaks were resized to 500 bps around their center. The 150,000 peaks containing the most read counts were maintained for downstream analysis. Hi-C data averaged from multiple cell types was utilized at a 5-Kb resolution (ENCFF134PUN). At this stage, the ABC model filtered ubiquitously expressed genes regardless of cell type. We generated ABC networks for each of the seven time points we sequenced.

We reconstructed time point-specific E-P-INs. First, a consolidated atlas was generated, correcting for sequencing disparities across samples. Enhancers located less than 100 bps away from each other were concatenated, leaving their time point-specific programs of regulation intact. These enhancers were then assigned searchable reference labels. Enhancers in this atlas were filtered according to two non-exclusive criteria. Regions that showed significant differential expression (p-values = 0.05) and temporal expression (p-value = 0.01) in our ATAC-seq and ChIP-seq H3K27ac were included in our E-P-INs. DESeq2 conducted differential expression [37] across possible pairings of time points, and ImpulseDE2 was used to determine temporal expression [52]. Regions in the 95th percentile of both ATAC-seq and ChIP-seq for H3K27ac counts were included in our E-P-INs regardless of dynamic expression. Promoters were processed similarly using our RNA-seq datasets. DEGs (p-value = 0.05) that demonstrated temporal expression (p-value = 0.01) were included in our E-P-INs. Promoters not included using the earlier mentioned threshold in the highest 5% of TPM values were included. This was done to capture highly expressed cellular processes during neural differentiation; 283 of the 19,825 candidate genes for our reconstructed E-P-INs are placed using the TPM threshold alone. Edges in our E-P-INs contain enhancers and promoters that both met the filtration criteria. Using this technique, we compiled our seven time point-specific E-P-INs and the collapsed network. The package pandas was utilized to store and manipulate these networks [53].

Enrichment in neural differentiation signals

We characterized our E-P-INs using two distinct forms of investigation. Using GSEA and GREAT, we performed enrichment analyses of enhancers and promoters, respectively [10,11,12]. We captured their biological processes, cellular components, and molecular functions. Using Jaccard similarities, we observed our E-P-INs morph over time.

GSEA is a functional-class-scoring method for pathway analysis. In GSEA, a predefined gene list is used to evaluate whether gene sets linked through shared biological functions are overrepresented in one of two compared conditions. It focuses on the collective patterns of these sets rather than changes in single genes. We used GSEA to examine the Gene Ontology (GO) collections; we provided the computational tool with the 500 most occurring promoters within our networks. Promoters with multiple predicted E-P-Is occurred more frequently in our E-P-INs. GREAT is an overrepresentation-analysis method. GREAT provided GO enrichment for the nearest genes to our enhancer regions. GREAT compares the proportions of gene sets in these enhancers to those in a background model, in this case, the entire human genome, to evaluate enrichment. In both cases, we report p-values from these tests. Visualizations are created using the Python packages Matplotlib [54] and Seaborn [55]. These visualizations use the − log(p-value) of GSEA and GREAT outcomes, positing the most significant GO in the chart’s longest and darkest bars.

Jaccard similarities are statistics used to measure overlaps between two sets. Jaccard similarities are a percentage of the intersection in two sets compared to the union. These scores serve as the foundation for our cross-network comparisons. All-inclusive sets of promoters, enhancers, and their edges from our E-P-INs are used to calculate Jaccard similarities, comparing the time point-specific E-P-INs across these three metrics.

Clustering reveals regulatory substructures

Revealing substructures tied to an enhancer’s regulation of genes necessitated using an algorithm for clustering built around betweenness centrality. This metric divides graphs using a path of their shortest connections, separating dense clustering in communities [13]. This approach, combined with our bipartite graphs’ high modularity, creates subnetworks between enhancers and promoters that we believe mimic regulatory profiles. Girvan-Newman clustering was employed using Python’s iGraph implementation [56].

After creating communities, we categorized them into four distinct categories: 1NR, 1R, 2NR, and 2R substructures. These categories represent redundant and non-redundant regulatory programs. 1NR entails a single gene regulated through a sole enhancer. 1R substructures see a single gene regulated by multiple enhancers. 2NR substructures have one enhancer responsible for regulating multiple genes. 2R communities consist of promoters and enhancers, all with redundancies, interacting in a large network. Distributions reported are percentages of substructures over the total communities. We reported average enhancer, promoter, and E-P-Is per community in each substructure.

Enrichment for disease-associated variants

Enrichment for disease-associated variation was assessed using in-house permutation testing and distinct sets of SNVs and genes. These permutation tests allowed for randomizable attempts and performed 1,000 iterations. SNV tests are employed using the BedTools Suite [57].

ASD-associated genes and DD-associated genes [17] both had an empirical p-value < 0.001 in their respective permutation tests. These tests determine whether our E-P-INs, built using dynamic and temporal CREs, demonstrate enrichment for these disease-associated genes. Using the ABC Model promoters before differential and temporal filtration as our null distribution, we extract unique gene lists equal in size to those of the E-P-IN we are testing in each iteration. The empirical p-value increases as higher counts of disease-associated genes are found in randomized null lists than in the list sourced from our E-P-IN. An analogous test is conducted for enhancers, as randomized lists of enhancers from our ABC model networks serve as a null distribution, intersecting with variants. Each iteration that counts more variant hits than our target E-P-IN increased the empirical p-value.

Additionally, we examined ASD-associated SNV enrichment through two additional methods. First, we calculated the percentage of communities in each substructure that contained an AS-associated de novo mutation. Next, we devised a permutation test to evaluate if ASD-associated SNV-containing enhancers interact with ASD-associated promoters with significance. We shuffled edges between enhancers and promoters in our network using filtered E-P-Is as the null distribution. We calculated an empirical p-value, which increased when more E-P-Is in our null distribution contained an enhancer with an ASD-associated SNV regulating an ASD-associated promoter. The permutation test employs 1000 permutations.

Validating enhancers with transcription factor overexpression

TF binding was predicted in our enhancers using FIMO [33], a computational tool in ENCODE’s MEME suite. We sourced position weight matrices from ENCODE [34] and JASPAR databases [35]. The consolidated atlas of enhancers was formatted as a BED file, converted to FASTA using the hg19 reference genome [58], and input into FIMO. Using a significance threshold of 1e − 5, we located TF binding sites (TFBSs) predicted to fall in our enhancers. TFs were then extrapolated from TFBSs.

We used Fisher’s exact test to measure the distribution of DEGs and non-DEGs in each enhancer that bound one of the six TFs for which we collected overexpression RNA-seq datasets. The Benjamini–Hochberg procedure was utilized for multiple testing corrections. FDR values below 0.05 were deemed statistically significant. Here, DEGs found downregulated (log2 fold change > 0) were not considered [7, 59,60,61,62,63]; these TFs function primarily as transcriptional activators. The p-value generated was used to find the enhancers’ relation to the overexpressed TF. DEGs were found at a p-value below 0.01 in our DESeq2 analysis.

We developed a permutation test to examine if enhancers predicted to bind overexpressed TFs regulated DEGs more frequently than enhancers that did not. We calculated the promotion of DEGs in both TF-binding and non-binding enhancers. Next, the DEG label is shuffled across genes, and the proportions are recalculated. The empirical p-value denotes the number of occurrences where the shuffled DEGs are found downstream of predicted non-TF-binding enhancers more often than enhancers thought to bind the overexpressed TFs. The permutation test iterates 1000 times before generating an empirical p-value. Here, the directionality of differential expression was not utilized; our analysis is focused on examining how these TFs drive neural differentiation. Examining genes impacted by co-binding or antagonistic TFs contributes to a better understanding of TFs’ role in our model system.

Validating enhancers with massively parallel reporter assays

We recently published perturbation MPRA datasets [7, 41, 42], systematically perturbing binding motifs across hundreds of enhancers and validating their functionality at the same time points in early neural differentiation as our E-P-INs’ model. We used this data to validate TFs predicted to bind by FIMO and enhancers in our networks [7]. This validation comes with stringent criteria, through which we validated a small subset of our E-P-INs.

Our MPRA-filtered E-P-INs overlap with perturbation-based datasets. MPRA-tested regions need to intersect with candidate enhancers included in assembled E-P-INs; this intersection needs to be in the proper context; the enhancer must have regulatory activity simultaneously when the MPRA is perturbing a motif. Finally, the perturbed motif also needs to be supported using our FIMO analysis.

We created a permutation test to investigate whether MPRA-validated enhancers were enriched in our filtered networks. We calculated the number of enhancers overlapping MPRA-tested sequences in 1000 random instances of our ABC model networks; we calculated an empirical p-value to show the count of these permutations that surpassed the number of overlaps in our collapsed E-P-IN.

Using this smaller network, we derived a similarity metric that is Jaccard-based. We generated sets of enhancers, downstream genes, and time points that TFs bind to, regulate, and are functional. These sets create individual Jaccard scores between our TFs, which are averaged at equal weights < into a composite score. These scores are then hierarchically clustered with SciPy [64].

Availability of data and materials

Computation methods and datasets necessary to reproduce this study, including the tool E-P-INAnalyzer, are available through GitHub [48] and Zenodo [47]. This codebase is distributed under the GNU General Public License v3.0.

References

  1. Li L, Wunderlich Z. An enhancer’s length and composition are shaped by its regulatory task. Front Genet. 2017;8:63.

    Article  PubMed  PubMed Central  Google Scholar 

  2. Schoenfelder S, Fraser P. Long-range enhancer-promoter contacts in gene expression control. Nat Rev Genet. 2019;20(8):437–55.

    Article  CAS  PubMed  Google Scholar 

  3. Pennacchio LA, Bickmore W, Dean A, Nobrega MA, Bejerano G. Enhancers: five essential questions. Nat Rev Genet. 2013;14(4):288–95.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  4. Le NQK, Yapp EKY, Nagasundaram N, Yeh HY. Classifying promoters by interpreting the hidden information of DNA sequences via deep learning and combination of continuous FastText N-Grams. Front Bioeng Biotechnol. 2019;7:305.

    Article  PubMed  PubMed Central  Google Scholar 

  5. Heinz S, Romanoski CE, Benner C, Glass CK. The selection and function of cell type-specific enhancers. Nat Rev Mol Cell Biol. 2015;16(3):144–54.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  6. Nasser J, Bergman DT, Fulco CP, Guckelberger P, Doughty BR, Patwardhan TA, et al. Genome-wide enhancer maps link risk variants to disease genes. Nature. 2021;593(7858):238–43.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  7. Inoue F, Kreimer A, Ashuach T, Ahituv N, Yosef N. Identification and massively parallel characterization of regulatory elements driving neural induction. Cell Stem Cell. 2019;25(5):713–727.e10.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  8. Fulco CP, Nasser J, Jones TR, Munson G, Bergman DT, Subramanian V, et al. Activity-by-contact model of enhancer-promoter regulation from thousands of CRISPR perturbations. Nat Genet. 2019;51(12):1664–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  9. Kvon EZ, Waymack R, Gad M, Wunderlich Z. Enhancer redundancy in development and disease. Nat Rev Genet. 2021;22(5):324–36.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  10. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A. 2005;102(43):15545–50.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  11. McLean CY, Bristor D, Hiller M, Clarke SL, Schaar BT, Lowe CB, et al. GREAT improves functional interpretation of cis-regulatory regions. Nat Biotechnol. 2010;28(5):495–501.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  12. Tanigawa Y, Dyer ES, Bejerano G. WhichTF is functionally important in your open chromatin data? PLoS Comput Biol. 2022;18(8): e1010378.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  13. Girvan M, Newman MEJ. Community structure in social and biological networks. Proc Natl Acad Sci USA. 2002;99(12):7821–6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  14. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13(11):2498–504.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  15. Ma Y, Bendl J, Hartley BJ, Fullard JF, Abdelaal R, Ho SM, et al. Activity-dependent transcriptional program in NGN2+ neurons enriched for genetic risk for brain-related disorders. Biol Psychiatry. 2023;S0006–3223(23):01426–9.

    Google Scholar 

  16. An JY, Lin K, Zhu L, Werling DM, Dong S, Brand H, et al. Genome-wide de novo risk score implicates promoter variation in autism spectrum disorder. Science. 2018;362(6420):eaat6576.

    Article  PubMed  PubMed Central  Google Scholar 

  17. Fu JM, Satterstrom FK, Peng M, Brand H, Collins RL, Dong S, et al. Rare coding variation provides insight into the genetic architecture and phenotypic context of autism. Nat Genet. 2022;54(9):1320–31.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Bipolar Disorder and Schizophrenia Working Group of the Psychiatric Genomics Consortium. Electronic address: douglas.ruderfer@vanderbilt.edu, Bipolar Disorder and Schizophrenia Working Group of the Psychiatric Genomics Consortium. Genomic dissection of bipolar disorder and schizophrenia, including 28 subphenotypes. Cell. 2018;173(7):1705–1715.e16.

    Article  Google Scholar 

  19. Coleman JRI, Peyrot WJ, Purves KL, Davis KAS, Rayner C, Choi SW, et al. Genome-wide gene-environment analyses of major depressive disorder and reported lifetime traumatic experiences in UK Biobank. Mol Psychiatry. 2020;25(7):1430–46.

    Article  PubMed  PubMed Central  Google Scholar 

  20. Simón-Sánchez J, Schulte C, Bras JM, Sharma M, Gibbs JR, Berg D, et al. Genome-wide association study reveals genetic risk underlying Parkinson’s disease. Nat Genet. 2009;41(12):1308–12.

    Article  PubMed  PubMed Central  Google Scholar 

  21. Jansen IE, Savage JE, Watanabe K, Bryois J, Williams DM, Steinberg S, et al. Genome-wide meta-analysis identifies new loci and functional pathways influencing Alzheimer’s disease risk. Nat Genet. 2019;51(3):404–13.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  22. Sanders B, D’Andrea D, Collins MO, Rees E, Steward TGJ, Zhu Y, et al. Transcriptional programs regulating neuronal differentiation are disrupted in DLG2 knockout human embryonic stem cells and enriched for schizophrenia and related disorders risk variants. Nat Commun. 2022;13(1):27.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  23. Nishioka M, Takayama J, Sakai N, Kazuno A-A, Ishiwata M, Ueda J, et al. Deep exome sequencing identifies enrichment of deleterious mosaic variants in neurodevelopmental disorder genes and mitochondrial tRNA regions in bipolar disorder. Mol Psychiatry. 2023;28(10):4294–306.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Wang X, Goldstein DB. Enhancer domains predict gene pathogenicity and inform gene discovery in complex disease. Am J Hum Genet. 2020;106(2):215–33.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  25. Moen MJ, Adams HHH, Brandsma JH, Dekkers DHW, Akinci U, Karkampouna S, et al. An interaction network of mental disorder proteins in neural stem cells. Transl Psychiatry. 2017;7(4):e1082.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  26. Carvill GL, Heavin SB, Yendle SC, McMahon JM, O’Roak BJ, Cook J, et al. Targeted resequencing in epileptic encephalopathies identifies de novo mutations in CHD2 and SYNGAP1. Nat Genet. 2013;45(7):825–30.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  27. Weissberg O, Elliott E. The mechanisms of CHD8 in neurodevelopment and autism spectrum disorders. Genes (Basel). 2021;12(8):1133.

    Article  CAS  PubMed  Google Scholar 

  28. Bina R, Matalon D, Fregeau B, Tarsitano JJ, Aukrust I, Houge G, et al. De novo variants in SUPT16H cause neurodevelopmental disorders associated with corpus callosum abnormalities. J Med Genet. 2020;57(7):461–5.

    Article  PubMed  Google Scholar 

  29. Kato Y, Shirai R, Ohbuchi K, Oizumi H, Yamamoto M, Miyata W, et al. Hesperetin ameliorates inhibition of neuronal and oligodendroglial cell differentiation phenotypes induced by knockdown of Rab2b, an autism spectrum disorder-associated gene product. Neurol Int. 2023;15(1):371–91.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Prontera P, Ottaviani V, Toccaceli D, Rogaia D, Ardisia C, Romani R, et al. Recurrent ∼100 Kb microdeletion in the chromosomal region 14q11.2, involving CHD8 gene, is associated with autism and macrocephaly. Am J Med Genet A. 2014;164A(12):3137–41.

    Article  PubMed  Google Scholar 

  31. Terrone G, Cappuccio G, Genesio R, Esposito A, Fiorentino V, Riccitelli M, et al. A case of 14q11.2 microdeletion with autistic features, severe obesity and facial dysmorphisms suggestive of Wolf-Hirschhorn syndrome. Am J Med Genet Part A. 2014;164(1):190–3.

    Article  CAS  Google Scholar 

  32. Panigrahi A, O’Malley BW. Mechanisms of enhancer action: the known and the unknown. Genome Biol. 2021;22(1):108.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  33. Grant CE, Bailey TL, Noble WS. FIMO: scanning for occurrences of a given motif. Bioinformatics. 2011;27(7):1017–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  34. Kheradpour P, Kellis M. Systematic discovery and characterization of regulatory motifs in ENCODE TF binding experiments. Nucleic Acids Res. 2014;42(5):2976–87.

    Article  CAS  PubMed  Google Scholar 

  35. Rauluseviciute I, Riudavets-Puig R, Blanc-Mathieu R, Castro-Mondragon JA, Ferenc K, Kumar V, et al. JASPAR 2024: 20th anniversary of the open-access database of transcription factor binding profiles. Nucleic Acids Res. 2024;52:D174–82. https://doi.org/10.1093/nar/gkad1059.

    Article  PubMed  Google Scholar 

  36. Zhang Q, Zhang Y, Wang C, Xu Z, Liang Q, An L, et al. The zinc finger transcription factor Sp9 Is required for the development of striatopallidal projection neurons. Cell Rep. 2016;16(5):1431–44.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  37. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550.

    Article  PubMed  PubMed Central  Google Scholar 

  38. Zhang X, Huang CT, Chen J, Pankratz MT, Xi J, Li J, et al. Pax6 is a human neuroectoderm cell fate determinant. Cell Stem Cell. 2010;7(1):90–100.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  39. Tepper RG, Ashraf J, Kaletsky R, Kleemann G, Murphy CT, Bussemaker HJ. PQM-1 complements DAF-16 as a key transcriptional regulator of DAF-2-mediated development and longevity. Cell. 2013;154(3):676–90.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  40. Ashuach T, Fischer DS, Kreimer A, Ahituv N, Theis FJ, Yosef N. MPRAnalyze: statistical framework for massively parallel reporter assays. Genome Biol. 2019;20(1):183.

    Article  PubMed  PubMed Central  Google Scholar 

  41. Kreimer A, Ashuach T, Inoue F, Khodaverdian A, Deng C, Yosef N, et al. Massively parallel reporter perturbation assays uncover temporal regulatory architecture during neural differentiation. Nat Commun. 2022;13(1):1504.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  42. Liu J, Ashuach T, Inoue F, Ahituv N, Yosef N, Kreimer A. Optimizing sequence design strategies for perturbation MPRAs: a computational evaluation framework. Nucleic Acids Res. 2024;52(4):1613–27.

    Article  PubMed  PubMed Central  Google Scholar 

  43. Stanfel MN, Moses KA, Schwartz RJ, Zimmer WE. Regulation of organ development by the NKX-homeodomain factors: an NKX code (Noisy-le-grand). Cell Mol Biol. 2005;Suppl 51:785–99.

    Google Scholar 

  44. Yanai H, Negishi H, Taniguchi T. The IRF family of transcription factors: Inception, impact and implications in oncogenesis. Oncoimmunology. 2012;1(8):1376–86.

    Article  PubMed  PubMed Central  Google Scholar 

  45. Claringbould A, Zaugg JB. Enhancers in disease: molecular basis and emerging treatment strategies. Trends Mol Med. 2021;27(11):1060–73.

    Article  CAS  PubMed  Google Scholar 

  46. Koesterich J, An JY, Inoue F, Sohota A, Ahituv N, Sanders SJ, et al. Characterization of De Novo promoter variants in autism spectrum disorder with massively parallel reporter assays. Int J Mol Sci. 2023;24(4):3509.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  47. DeGroat W, Kreimer A. E-P-INAnalyzer. Zenodo. 2024. https://doi.org/10.5281/zenodo.13119427

  48. DeGroat W, Kreimer A. E-P-INAnalyzer. GitHub 2024. https://github.com/KreimerLab/E-P-INAnalyzer

  49. Van Rossum G, Drake FL. Python 3 Reference Manual. Scotts Valley: CreateSpace; 2009.

    Google Scholar 

  50. R Core Team. R: A Language and Environment for Statistical Computing. Vienna: R Foundation for Statistical Computing; 2021. Available at: https://www.R-project.org/.

  51. Zhang Y, Liu T, Meyer CA, Eeckhoute J, Johnson DS, Bernstein BE, et al. Model-based Analysis of ChIP-Seq (MACS). Genome Biol. 2008;9(9):R137.

    Article  PubMed  PubMed Central  Google Scholar 

  52. Fischer DS, Theis FJ, Yosef N. Impulse model-based differential expression analysis of time course sequencing data. Nucleic Acids Res. 2018;46(20):e119.

  53. McKinney W. Data Structures for Statistical Computing in Python. Presented at the Python in Science Conference. Austin; 2010. p. 56–61. https://doi.org/10.25080/Majora-92bf1922-00a.

  54. Hunter JD. Matplotlib: A 2D graphics environment. Comput Sci Engine. 2007;9(3):90–5.

    Article  Google Scholar 

  55. Waskom ML. seaborn: statistical data visualization. J Open Source Software. 2021;6(60):3021.

    Article  Google Scholar 

  56. Csardi G, Nepusz T. The Igraph Software Package for Complex Network Research. InterJournal, Complex Systems. 2005;1695(11).

  57. Quinlan AR, Hall IM. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics. 2010;26(6):841–2.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  58. Church DM, Schneider VA, Graves T, Auger K, Cunningham F, Bouk N, et al. Modernizing reference genome assemblies. Plos Biol. 2011;9(7):e1001091.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  59. Bulfone A, Menguzzato E, Broccoli V, Marchitiello A, Gattuso C, Mariani M, et al. Barhl1, a gene belonging to a new subfamily of mammalian homeobox genes, is expressed in migrating neurons of the CNS. Hum Mol Genet. 2000;9(9):1443–52.

    Article  CAS  PubMed  Google Scholar 

  60. Dou Z, Son JE, Hui CC. Irx3 and Irx5 - novel regulatory factors of postnatal hypothalamic neurogenesis. Front Neurosci. 2021;15:763856.

    Article  PubMed  PubMed Central  Google Scholar 

  61. Pillai A, Mansouri A, Behringer R, Westphal H, Goulding M. Lhx1 and Lhx5 maintain the inhibitory-neurotransmitter status of interneurons in the dorsal spinal cord. Development. 2007;134(2):357–66.

    Article  CAS  PubMed  Google Scholar 

  62. Huang B, Li X, Tu X, Zhao W, Zhu D, Feng Y, et al. OTX1 regulates cell cycle progression of neural progenitors in the developing cerebral cortex. J Biol Chem. 2018;293(6):2137–48.

    Article  CAS  PubMed  Google Scholar 

  63. Kim EA, Noh YT, Ryu MJ, Kim HT, Lee SE, Kim CH, et al. Phosphorylation and transactivation of Pax6 by homeodomain-interacting protein kinase 2 *. J Biol Chem. 2006;281(11):7489–97.

    Article  CAS  PubMed  Google Scholar 

  64. Virtanen P, Gommers R, Oliphant TE, Haberland M, Reddy T, Cournapeau D, et al. SciPy 1.0: fundamental algorithms for scientific computing in Python. Nat Methods. 2020;17(3):261–72.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

Download references

Acknowledgements

Not applicable.

Review history

The review history is available as Additional file 7.

Peer review information

Tim Sands was the primary editor of this article and managed its editorial process and peer review in collaboration with the rest of the editorial team.

Funding

This work was supported in part by the National Institute of Mental Health grant R00MH117393-05 (A.K.).

Author information

Authors and Affiliations

  1. Center for Advanced Biotechnology and Medicine, Rutgers, The State University of New Jersey, 679 Hoes Lane West, Piscataway, NJ, 08854, USA

    William DeGroat & Anat Kreimer

  2. Institute for the Advanced Study of Human Biology (WPI-ASHBi), Kyoto University, Kyoto, Japan

    Fumitaka Inoue

  3. Department of Electrical Engineering and Computer Sciences and Center for Computational Biology, University of California, Berkeley, 387 Soda Hall, Berkeley, CA, 94720, USA

    Tal Ashuach

  4. Department of Systems Immunology, Weizmann Institute of Science, 234 Herzl Street, Rehovot, 7610001, Israel

    Nir Yosef

  5. Chan-Zuckerberg Biohub, 499 Illinois St, San Francisco, CA, 94158, USA

    Nir Yosef

  6. Department of Systems Immunology, Ragon Institute of MGH, MIT, and Harvard Institute of Science, 400 Technology Square, Cambridge, MA, 02139, USA

    Nir Yosef

  7. Department of Bioengineering and Therapeutic Sciences, University of California, 513 Parnassus Ave, San Francisco, CA, 94143, USA

    Nadav Ahituv

  8. Institute for Human Genetics, University of California, 513 Parnassus Ave, San Francisco, CA, 94143, USA

    Nadav Ahituv

  9. Department of Biochemistry and Molecular Biology, Rutgers, The State University of New Jersey, 604 Allison Road, Piscataway, NJ, 08854, USA

    Anat Kreimer

Authors
  1. William DeGroat
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Fumitaka Inoue
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Tal Ashuach
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Nir Yosef
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Nadav Ahituv
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Anat Kreimer
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

A.K. conceived and designed the analysis, collected data, wrote the manuscript, and developed analysis ideas and tools. W.D. performed analysis, wrote the manuscript, and developed analysis ideas and tools. T.I., T.A., N.Y., and N.A. contributed analysis ideas. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Anat Kreimer.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Supplementary figures S1 – S3 with descriptions.

13059_2024_3365_MOESM2_ESM.xlsx

Additional file 2: Supplementary table S1 - E-P-INs containing enhancers, promoters, and edges. Tables containing seven time-point-specific E-P-INs, a collapsed E-P-IN, and an MPRA-validated E-P-IN.

Additional file 3: Supplementary table S2 - GSEA findings in E-P-INs.

Additional file 4: Supplementary table S3 - GREAT findings in E-P-INs.

13059_2024_3365_MOESM5_ESM.xlsx

Additional file 5: Supplementary table S4 - Disease-associated variant enrichment in E-P-INs. Empirical p-values from our permutation tests for time-point-specific and collapsed E-P-INs.

13059_2024_3365_MOESM6_ESM.xlsx

Additional file 6: Supplementary table S5 - TF similarities. Three submatrices detailing enhancer, gene, and time point similarities between TF pairs in an MPRA-validated E-P-IN and the composite TF by TF similarity score matrix.

Additional file 7. Peer review history.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

DeGroat, W., Inoue, F., Ashuach, T. et al. Comprehensive network modeling approaches unravel dynamic enhancer-promoter interactions across neural differentiation. Genome Biol 25, 221 (2024). https://doi.org/10.1186/s13059-024-03365-w

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s13059-024-03365-w

Keywords

Genome Biology

ISSN: 1474-760X

Contact us