BRAIN-MAGNET: A novel functional genomics atlas coupled with convolutional neural networks facilitates clinical interpretation of disease relevant variants in non-coding regulatory elements

Adapted ChIP-STARR-seq approach

Compared to our previous application of ChIP-STARR-seq in human embryonic stem cells (ESCs) ⁹ that focussed on NCREs enriched for the TFs NANOG, OCT4, and the histone modifications H3K4me1 and H3K27ac (herein referred to as ESC ChIP-STARR-seq), we performed a number of updates on the experimental procedures, aiming to reduce experimental noise, increasing reproducibility and reducing experimental complexity and costs (Extended Data Fig. 1a).

First, to enrich for genomic regions with putative NCRE function in NSCs, we performed ChIP in duplicates for the TFs SOX2 and YY1, that are important in NSCs^13,14, and the histone modifications H3K27ac and H3K4me1, which are associated with putative NCREs¹⁵. ChIP- qPCR showed an enrichment at known NCREs surrounding the FGFR1, PAX6, SOX2 and NES loci (Extended Data Fig. 1b). Sequencing of the ChIP material (ChIP-seq) showed a good correlation between replicates (Pearson correlation coefficient > 0.88) (Extended Data Fig. 1c,d). Next, using the same ChIP-DNA, ChIP-STARR-seq plasmid libraries were generated from both ChIP replicates and pooled. To further decrease the number of transfections and simultaneously increase the number of individual plasmids covering each putative NCRE, theoretically reducing experimental noise, we pooled together the plasmid libraries derived from ChIP for the histone modifications H3K27ac and H3K4me1, and for the TFs SOX2 and YY1. Hereafter, we refer to these combined plasmid libraries as HIST and TF libraries, respectively, and the application of these libraries for ChIP-STARR-seq in NSCs or ESCs as NSC ChIP-STARR-seq or comparative ChIP-STARR-seq, respectively. Visual inspection of sequencing results to assess whether the pooled plasmid libraries are a good representation of the initial ChIP material, showed that peaks from ChIP-seq experiments were still recognisable in the sequencing tracks of the pooled plasmids (Extended Data Fig. 1e), and analysis of normalized sequencing reads confirmed a high correlation between initial ChIP-seq peaks and corresponding pooled plasmid libraries (Extended Data Fig. 1c,d). On average, 93.04% of HIST ChIP-seq peaks and 96.86% of TF ChIP-seq peaks were covered by more than 5 different plasmids and 46.91% of HIST peaks and 67.54% of TF peaks by more than 20 plasmids, with a mean of 27 and 43 distinct plasmids per HIST and TF ChIP- seq peak, respectively. Together this indicates an increased complexity of the generated plasmid libraries compared to our earlier study using ChIP-STARR-seq in ESCs⁹.

Second, to overcome the variability between samples that we previously observed in ESC ChIP-STARR-seq⁹, we implemented the following changes in the experimental strategy. We performed 5 independent transfections per plasmid pool, each starting with ∼3 million NSCs seeded 48 hours before the transfection (or ∼6 million ESCs transfected in suspension, see Methods). 24 hours post-transfection, we collected and pooled together approximately 30 million cells from these 5 independent transfections (Extended Data Fig. 1a). FACS analysis showed that 26.5% and 30% of NSCs transfected with the HIST or TF libraries, respectively, were GFP positive (Extended Data Fig. 1f).

Third, RNA was isolated from these pooled, unsorted cells and subsequently split and treated as two replicates, which were subjected to DNase treatment and cDNA synthesis. To limit PCR artefacts, the first PCR amplification step which suppresses residual plasmid contamination by the use of primers spanning a synthetic intron in the STARR-seq plasmid, was performed in 10 parallel reactions for each of the two replicates. These 10 PCRs per replicate were then pooled, purified and further amplified with the index primers required for sequencing, using 5 independent PCR reactions and two independent index primers for each replicate. Thus, we sequenced a total of four replicates of STARR-seq RNA for each pooled HIST or TF library, with an average of 20 million reads per sample (Extended Data Fig. 1a). As the two technical replicates amplified with independent index primers showed a high correlation (Pearson correlation coefficient >0.96) (Extended Data Fig. 1g-j), they were merged for further downstream data analysis.

Together this indicates that our revised experimental strategy generated highly complex plasmid libraries capturing genomic regions marked by H3K27ac, H3K4me1 or bound by SOX2 and YY1 in NSCs, which upon transfection and RNA sequencing results in highly correlated replicate data.

Differential NCRE activity analysis in NSCs and ESCs

To identify differences in NCRE activity between the two cell types, we next performed a differential NCRE activity analysis¹⁶, comparing NCREs captured in the same NSC-derived plasmid library with high activity in the one cell type and low activity in the other, and vice versa (Extended Data Fig. 4a). To this end, we first selected those NCREs that ranked in the top-10% of NSC ChIP-STARR-seq activity in NSCs and in the bottom-90% of comparative ChIP-STARR-seq in ESCs, and those NCREs that ranked in the top-10% of comparative ChIP-STARR-seq activity in ESCs and the bottom-90% in NSC ChIP-STARR-seq. With these settings, we identified 6,850 out of 14,818 NCREs (46.2%) overlapping between the top-10% of NSCs and ESCs (referred to as common-high) (Extended Data Fig. 4b), possibly explaining part of the overlap between TF motifs and TEs in NSC and ESC active NCREs (Fig. 1j and Extended Data Fig. 3c), and 7,968 out of 14,818 NCREs (53.8%) per cell type that were either NSC- or ESC-high. Visual inspection of the tracks surrounding the CHD2 locus shows both common-high NCREs and differentially active NCREs (Extended Data Fig. 4c). Investigating the distance distribution of these regions from the TSS, we observed that part of the NSC-high NCREs is located in closer proximity to TSSs, while common and ESC- high NCREs were spatially distributed in a more similar manner (Extended Data Fig. 4d). A deeper investigation revealed that all the three classes or NCREs are enriched mainly at intronic (NSC-high: 46.4%; ESC-high: 55.2%; common-high: 53.7%) and intergenic (NSC- high: 29.6%; ESC-high: 36.7%; common-high: 35.9%) regions with an increase at promoter- TSS regions for the NSC-high ones (NSC-high: 17.1%; ESC-high: 4.4%; common-high: 6.3%), in agreement with their tendency to be located in closer proximity to the TSS (Extended Data Fig. 4e). Next, we asked whether cell-type specific highly active NCREs were also more susceptible to LoF. As observed previously (Extended Data Fig. 2f,3b), the majority of the investigated regions are intolerant to LoF, with NSC-high NCREs showing a larger enrichment amongst intolerant regions, possibly indicating that alterations of these regions might be negatively selected in the human population or might be implicated in disease (Extended Data Fig. 4f). Investigating the pLI of the genes linked to these NCREs groups, we did not find a significant difference (Extended Data Fig. 4g).

Next, we aimed at understanding whether these NCREs might regulate cell-type specific gene expression. We observed that NSC-high NCREs are linked to genes expressed at a significantly higher level in NSCs but also in ESCs (Extended Data Fig. 4h), whereas genes linked to ESC-high and common-high NCREs did not show major expression differences. Gene ontology analysis identified for all three categories of NSC-high, ESC-high and common-high NCREs terms related to transcriptional regulation (Extended Data Fig. 4i). Common- high NCREs were also enriched for terms related to nervous system development and generation of neurons and axons. Like previously seen for all NCREs active upon comparative ChIP-STARR-seq in ESCs, also ESC-high NCREs were enriched for nervous system development and generation of neurons, rather than other expected terms for pluripotent stem cells, despite that these genes are not highly expressed in ESCs. This observation further provides support for the hypothesis that these ESC-high NCREs, despite showing episomal activity in comparative ChIP-STARR-seq in ESCs when using plasmid libraries obtained from NSCs, might not be active at the endogenous chromatin in ESCs. Rather, these ESC-high NCREs might be primed for future activation at the endogenous locus, and this priming might not be recapitulated in an episomal plasmid-based system.

To further investigate this hypothesis, we assessed various chromatin profiling data at the common and cell-type specific highly active NCREs. As expected, since all the investigated putative NCREs are derived from H3K27ac, H3K4me1, YY1 and SOX2 ChIP in NSCs, all the common and differentially active NCREs (ESC-high, NSC-high) show central enrichment of these chromatin marks across NSC ChIP-seq data (Extended Data Fig. 4j,k). To further characterize the 7,968 NSC-high NCREs at their endogenous locus, we retrieved H3K4me3 and H3K27me3 from previously published data of H9-derived NSCs⁵⁴ and performed ATAC- seq in NSCs to assess chromatin accessibility. We observed that the NSC-high NCREs are enriched for accessible chromatin, as expected given the chromatin immunoprecipitation in this cell type, and a strong enrichment of the active chromatin mark H3K4me3 (Extended Data Fig. 4j,k) overlapping the enrichment of YY1 and H3K27ac. H3K4me3, as well as YY1, is highly enriched at promoters⁵⁵, suggesting this fraction of loci might partially include the 17% of NSC-high NCREs found in proximity to promoters (as previously discussed, the STARR-seq design does not allow to test promoter activity and thus despite being located in proximity to promoters these sequences are NCREs). Compatible with the immunoprecipitation of markers associated with NCRE activity used to identify ChIP-STARR- seq NCREs, we did not see enrichment of the inactive chromatin mark H3K27me3 at these NCREs.

To investigate the hypothesis that ESC-high NCREs are not enriched in active marks at the endogenous chromatin loci, we collected H9 ESC epigenome data from publicly available sources and from our previously generated H3K27ac, H3K4me1, OCT4 and NANOG ChIP- seq⁹. Furthermore, we generated ChIP-seq data for YY1 and SOX2 in ESCs. ESC-high NCREs, even though located in accessible chromatin in ESCs, show no enrichment of the active chromatin marks H3K27ac and H3K4me1 in ESCs (Extended Data Fig. 4j,k). The majority of the ESC-high NCREs from comparative ChIP-STARR-seq in ESCs show a very mild enrichment of SOX2 and NANOG binding, but interestingly, a subset of the ESC-high NCREs (n=1,100 out of the 7,968 ESC-high NCREs) had an accumulation of the histone marks H3K4me2 and H3K4me3, with especially the latter mark showing a bi-modal distribution with a valley over the NCRE centre. Primed NCREs are indeed characterized by methylation of H3K4 and no acetylation of H3K27³⁷. Therefore, this observation is compatible with the hypothesis that these NCREs, despite being highly active in the episomal setting, are not active at the endogenous locus but primed for activation at a later stage. Further pointing to this scenario, we noticed that ESC-high NCREs primed with H3K4 methylation in ESCs, gain H3K27 acetylation at the NSC stage, although to a lower level compared to NSC-high NCREs, compatible with the lower activity of these regions in NSCs (Extended Data Fig. 4j,k). Interestingly, these regions maintain H3K4me3 also at later stages of differentiation, in neurons and in astrocytes, where they also display H3K27ac. The same scenario described for ESC-high NCREs applies to common-high NCREs, for which we could observe a fraction of regions (n=1,200 out of 6,850 common-high) with H3K4me2/3 accumulation and no H3K27ac in ESCs, but with H3K27ac in NSCs. H3K4me3 regions have been previously reported to be located mainly at the TSS⁵⁵. In line with this, we observed that the majority of the ESC-high and common-high NCREs marked by H3K4me3 were located in proximity of a TSS, within 5kb down- or upstream (Extended Data Fig. 4l), but as outlined as the STARR- seq design does not allow to capture promoter activity these are thus NCREs in proximity of the TSS. Moreover, when investigating the enrichment of TF binding motifs, the TATA-box was among the most enriched motifs in the H3K4me2/3 positive subset of both the ESC-high and common-high NCREs (Extended Data Fig. 4m). Further supporting the hypothesis that this subset of regions is primed for activation at a later stage, we noticed that the majority of them (>75%) in NSCs belong to the NCRE activity categories 4 and 5 (Extended Data Fig. 4n).

Analysing further the TF motifs enriched in these regions, we noticed that the most enriched in the common-high subset positive for H3K4me2/3 in ESCs were YY1, p73 and p53, possibly explaining their activity in NSCs where these regions are marked by H3K27ac. We also identified the OCT4 and OCT4-SOX2-TCF-NANOG⁵⁶ motifs at this subset of ESC-high NCREs, possibly contributing to the episomal activity of these NCREs. Interestingly, the repressor protein ZNF281 was enriched at a small subset of ESC-high, H3K4me2/3 positive NCREs (Extended Data Fig. 4n). ZNF281 is a known repressor protein expressed at high levels in ESCs but downregulated upon neural differentiation⁵⁷, and it will therefore be interesting to investigate whether ZNF281 or related factors might contribute to explain silencing of some of the ESC-high NCREs at the endogenous chromatin.

BRAIN-MAGNET eliminates a clinically suspected candidate disease causing non- coding variant

We applied the NCRE atlas to a patient from our clinic, which was clinically suspected of Mowat-Wilson syndrome, a severe neurodevelopmental disorder linked to dysfunction of the ZEB2 gene. Previous extensive genetic investigations, including trio whole exome and genome sequencing did not reveal a protein coding mutation in ZEB2, or any other likely disease-causing variant, but identified a heterozygous chr2:144469837C>T variant in intron 2 of ZEB2. This variant overlapped with a highly active NCRE from NSCs (category 5), was absent in gnomAD v4, had a CADD score of 21.0, an ncER score of 99.6 and was predicted by the RegulomeDB v.2⁵⁸ to contain a TF binding motif (Extended Data Fig. 14a). According to LINSIGHT⁴⁴, an algorithm that aims to predict which non-coding variants are likely to have deleterious fitness effects, the chr2:144469837C>T variant was predicted to be likely pathogenic (score 0.98). In contrast, according to BRAIN-MAGNET, this nucleotide had a low cb score (65.7^th percentile within the NCRE (cb_each); 52.4^th percentile within all NCREs (cb_all)), while a different ONECUT1 TF motif located 36 nucleotides upstream is predicted to mostly impact on activity of this NCRE. To investigate this functionally, we cloned the NCRE in the STARR-seq plasmid and introduced by site directed mutagenesis either the variant identified in the patient or a 27 bp deletion centred on the ONECUT1 motif. Upon transfection in NSCs, the ONECUT1 deletion completely abolished the NCRE activity, while the patient variant did not majorly affect NCRE activity, in agreement with the BRAIN-MAGNET predictions (Extended Data Fig. 14b). In addition, clinical RNA-seq of patient derived fibroblasts did not reveal any ZEB2 expression changes and a diagnostic Episignature did not reveal evidence of a Mowat-Wilson syndrome specific DNA-methylation pattern (data not shown). Even though we cannot completely exclude that the intronic ZEB2 variant could have an effect in other cell types, based on the current functional testing the variant was clinically classified as a variant of unknown significance which was considered unlikely to contribute to the phenotype which remains unexplained. This illustrates how BRAIN-MAGNET can be used to clinically eliminate potential candidate causing non-coding variants identified during filtering of WGS data sets.

Cell culture

H9-derived human Neural Stem Cells (NSC, Gibco) were cultured as previously described⁵⁹. Briefly, cells were seeded onto Geltrex (Gibco) coated plates and cultured in KnockOut DMEM/F12 (Gibco) supplemented with 2 mM L-Glutamine (Gibco), 2% StemPro neural supplement (Gibco), 20 ng/ml EGF (Peprotech), 20 ng/ml b-FGF (Peprotech), 100 U/ml penicillin and 100 µg/ml streptomycin. Cells resulted negative at mycoplasma testing.

H9 ESCs were cultured as previously described⁵⁹ on Matrigel (Corning) coated plates in MTeSR-1 medium (STEMCELL Technologies).

Chromatin immunoprecipitation

For chromatin immunoprecipitation, 5x10⁷ NSC were harvested in 9 mL of PBS and dual cross-linked first with 2mM Di(N-succinimidyl) glutarate (Sigma) for 45 min followed by 10 minutes 1.1% buffered formaldehyde (50mM Hepes-KOH pH7.6, 100mM NaCl, 1mM EDTA pH8.0, 0.5mM EGTA pH8.0, 11% formaldehyde). Formaldehyde was then quenched with 125 mM glycine (Sigma) for 10 min. Cross-linked cell pellets were snap-frozen and stored at - 80°C.

ChIP was performed as described⁹. Briefly, cell nuclei were sonicated at 4°C with Diagenode Bioruptor for 45 cycles (30 sec on, 30 sec off), prior to overnight incubation with antibody- coupled beads. The antibodies used were rabbit-anti-H3K4me1 (ab8895, Abcam), rabbit-anti- H3K27ac (ab4729, Abcam), rabbit-anti-YY1 (ab109237, Abcam) and goat-anti-SOX2 (AF2018, R&D systems). Immunoprecipitated chromatin and total input control were de-cross- linked in a final volume of 100 µL, by addition of 3 µL of 10% SDS and 10 µL proteinase K (10 mg/ml, Roche) and 5µl RNase (11119915001, 500 ng/µl, Roche) and incubation overnight at 65°C at 1400 RPM on a shaking thermomixer block. Following Phenol:Chloroform:IAA (Sigma) extraction and ethanol precipitation, the concentration of ChIP-DNA was determined by Qubit (Invitrogen) and the sonication was assessed by gel-electrophoresis of total input DNA (target fragment size between 200 and 600 bp).

ChIP-qPCR

For qPCR, 10µl of ChIP material and total input material were diluted 6 times. 2 µL of DNA were used per qPCR reaction, using iTaq universal SYBR Green Supermix in a CFX96RTS thermal cycler (Bio-Rad), following manufacturer’s instructions. Fold enrichment was determined following the ΔΔct method. All data shown are averages of at least 2 biological replicates and 2 technical replicates.

ChIP-seq, ChIP-STARR-seq plasmid library preparation

ChIP-seq and ChIP-STARR-seq plasmid libraries were generated as described⁹. Briefly, 10 ng of ChIP DNA were end-repaired, dA tailed and adaptor ligated using NEBNext ChIP-seq library preparation kit (E6240, NEB), according to manufacturer’s instructions. Adaptor ligated DNA was eluted into 32 µL of 0.1xTE, 25 µL were used for ChIP-seq library preparation and 5 µl were processed for ChIP-STARR-seq plasmid libraries preparation.

For ChIP-seq library preparation, purified adaptor ligated DNA was PCR amplified with Illumina index primers (E7335 and E7500, NEB) according to manufacturer’s instructions. ChIP-seq libraries were assessed on an Agilent Tapestation. All sequencing occurred on an Illumina HiSeq 2500 platform, using 50bp single-end sequencing.

For details about the plasmid library preparation see Barakat et al⁹. Briefly, 5 µL of purified adaptor ligated DNA for each replicate were separately amplified with primers 147 STARRseq libr FW and 148 STARRseq libr RV¹¹, and used in a Gibson assembly reaction with an AgeI- SalI digested STARR-seq plasmid and a home-made Gibson reaction buffer (100mM Tris- HCl, 10mM MgCl2, 0.2 mM dNTP, 0.5U Phusion DNA polymerase (NEB), 0.16U 50 T5 exonuclease (NEB)). Gibson reactions for the ChIP replicates were pooled together, purified by ethanol precipitation and used for electroporation into electrocompetent MegaX DH10b E. coli bacteria (C640003, Invitrogen), according to manufacturer’s instructions. A total of 5 electroporation reactions per library were performed with 2 µl of DNA each. After recovery, bacteria from the 5 reactions were pooled together and 100 µL of a 1:100 and 1:10000 dilution was plated on Ampicillin containing Agar plates to enable estimation of the number of transformants. The remaining 5 mL of bacteria culture were incubated in 1 L of LB-media supplemented with Ampicillin over-night, and the plasmid DNA was isolated using a Maxiprep kit (Macherey-Nagel) according to manufacturer’s instructions.

Transfection of plasmid libraries

For each transfection in NSC ChIP-STARR-seq, approximately 3 million NSC were seeded in a 10 cm dish and 48 hours later transfected with 6.8 µg of plasmid library DNA and 27.5 µl Lipofectamine Stem Transfection Reagent (STEM00015, Invitrogen) following manufacturer’s instructions. In total, 5 dishes of transfected cells were used. For ESCs transfections in ESC ChIP-STARR-seq, approximately 6 million cells were seeded in a 10 cm dish and transfected in suspension with an identical transfection mix to NSCs in the presence of 5µM ROCK-inhibitor (Y27632, Millipore). After 24h of transfection, cells were collected, and RNA was extracted with Tri-reagent (Sigma) following the manufacturer’s instruction. A small fraction of cells was kept for estimating the fraction of GFP positive cells by flow cytometry.

ChIP-STARR-seq RNA and plasmid sequencing

ChIP-STARR-seq RNA libraries were prepared as described⁹, with minor modifications. Briefly, RNA was isolated with Trizol (Sigma), following manufacturer’s instruction and split into two replicates. The mRNA fraction was captured from each replicate of total RNA using Oligo (dT) 25 beads (61002, Life Technologies), DNaseI treated (18068-015, Life Technologies), and reverse transcribed with a GFP-mRNA specific primer (149 STARRseq rep RNA cDNA synth) following Superscript III protocol (18080-044, Life Technologies). To avoid plasmid DNA contamination, cDNA was PCR amplified in 10 parallel reactions for each replicate using primers spanning a synthetic intron (152 STARR reporter specific primer 2 fw and 153 STARR reporter specific primer 2 rv), as previously described¹¹. 35 ng of PCR amplified DNA were then used in a second round of PCR to add Illumina index primers (E7335, E7500, NEB). PCR was performed in 5 parallel reactions and using 2 independent index primers for each of the two replicates. In total, we sequenced 4 replicates of STARR-seq RNA for each pooled HIST or TF plasmid library. Corresponding plasmid libraries were similarly amplified using two independent Illumina index primers for each plasmid library. Quantity and quality of generated sequencing libraries was assessed on an Agilent Tapestation. All sequencing occurred on an Illumina HiSeq 2500 platform, using 50bp paired-end sequencing.

ATAC-seq

The Assay for Transposase Accessible Chromatin coupled to high-throughput sequencing (ATAC-seq) was performed as previously described⁶⁰ with slight modifications. Briefly, approximately 50.000 single cells were lysed in Resuspension buffer containing 0,1% IgePal, 0,1% Tween-20 and 0,01% Digitonin. Immediately after lysis, the buffer was washed out by Resuspension buffer with only Tween-20. Nuclei were centrifuged at 500xg for 10 min at 4°C. Nuclei were re-suspended in 25 μl 2x TD buffer (Illumina), 2.5 μl TDE1 (transposase, Illumina),

16.5 μl PBS, 0.5 μl Tween-20, 0.5 μl Digitonin, and incubated for 30 min at 37°C. Afterward, the sample was purified using the MinElute PCR Purification Kit (QIAGEN) according to the manufacturer’s protocol and eluted with 10 μl elution buffer. Samples were amplified using the NEBNext High Fidelity PCR master mix (New England Biolabs) and afterwards purified using a MinElute PCR purification kit (Qiagen) according to the manufacturer’s instructions. One microliter was loaded on an Agilent Technologies 2100 Bioanalyzer using a DNA 1000 assay to determine the library concentration and for quality check.

Cluster generation was performed according to the Illumina TruSeq SR Rapid Cluster kit v2 (cBot) Reagents Preparation Guide (www.illumina.com). Briefly, for sequencing libraries were pooled together to get a stock of 10 nM. One microliter of the 10 nM stock was denatured with NaOH, diluted to 10 pM and hybridized onto the flowcell. The hybridized products were sequentially amplified, linearized and end-blocked according to the Illumina Single Read Multiplex Sequencing user guide. After hybridization of the sequencing primer, sequencing- by-synthesis was performed using the HiSeq 2500 with paired end 50-cycle protocol followed by dual index sequencing. ATAC-seq libraries were sequenced to a depth of 20x10⁶ reads.

RNA-seq

For RNA sequencing of ESCs, we used our previously published data⁵⁹. For RNA sequencing of NSCs, RNA was isolated from two independent cultures following TRI reagent (Sigma) manufacturer’s instruction. RNA was further purified using column purification (Qiagen, #74204). mRNA capture, library prep and barcoding were performed according to standard procedures of the Erasmus MC Biomics facility. Libraries were then single end sequenced, 50bp, on an Illumina HiSeq2500 machine at a depth of approximately 20 million reads per sample.

Data processing

The FASTQ files of RNA-, ChIP-, ATAC- and STARR-seq data were trimmed using Trimmomatic (version 0.39)⁶¹ to remove possible adaptor contamination and remove low quality reads. The trimmed RNA- and ATAC-seq data were mapped to the GRCh38/hg38 human genome build using the HISAT2 aligner (version 2.2.1)^62,63. For RNA-seq, aligned reads were counted for each gene using htseq-count (version 0.12.4)⁶⁴ and further processed using edgeR⁶⁵ to normalize the data and perform downstream analysis. Bowtie2 (version 2.4.2) (with --very-sensitive parameter)⁶⁶ was used to align both ChIP- and STARR-seq data against the GRCh38/hg38 human genome build. Only properly and uniquely mapped reads, with mapping quality more than 30 (MAPQ >=30), were kept followed by removing possible duplicated reads (only for ChIP-seq data) using Picard’s MarkDuplicates (version 4.0.1.1) (http://broadinstitute.github.io/picard/). The aligned reads were converted to bigwig format using ”bamCoverage” function with --binSize 10 and --normalizeUsing CPM parameters, the read coverage was then computed using “multiBigwigSummary” and pearson correlation between replicates was calculated using “plotCorrelation” via deepTools (version 3.4.3)⁶⁷. Peak calling was performed using MACS2 (version 2.1.2) (with -q 0.05)⁶⁸ for each replicate of ESCs and NSCs separately, and took sequencing data of total input as background. Peaks overlapping blacklist regions defined by the ENCODE project (Accession: ENCFF356LFX, Date: 05.05.2020) were removed. The overlapping peaks between replicates were merged as common peaks using intersectBed (version2.30.0)^69,70.

NCRE activity analysis

To generate a non-redundant set of scaffold regions to assess activity of NCREs, ChIP-seq and RNA-plasmid peaks (H3K4me1, H3K27ac, YY1 and SOX2 of ChIP-seq and histone and TF plasmid DNA-seq libraries of STARR-seq) were merged after extending the peak summit to a width of 1000 bp (500 bp from each side). All regions larger than 1000 bp were split in half until all regions were at most 1000 bp in length, avoiding large regions and preserving high genomic resolution for later analysis. All further analysis for both NSC ChIP-STARR-seq and comparative ChIP-STARR-seq was performed on these scaffold regions.

To define active regions, reads under the defined scaffolds in both NSC and comparative ChIP-STARR-seq in ESCs were counted by FeatureCount⁷¹. Read counts for each scaffold region were normalized by the total number of aligned reads in each library, after removing low coverage regions with less than 20 reads in at least two samples. The average number of normalized reads between replicates was taken for computing NCRE activity followed by the ratio of RNA-seq and DNA-seq (plasmid) reads as measurement of NCRE activity. The scaffold regions were ranked based on the average activity of RNA-seq samples in 5 different, equally sized, categories and the top-10% regions were defined as the most active NCREs in each cell line.

As NSC ChIP-STARR-seq we define those experiments that used NSC-derived plasmid libraries transfected in NSCs. As comparative ChIP-STARR-seq we refer to those data obtained from the transfection of NSC-derived plasmid libraries in ESCs. When discussing ESC ChIP-STARR-seq, we refer to our previous study that generated data obtained from transfections of ESC-derived plasmid libraries in ESCs⁹. Furthermore, we defined cell-specific regions as the overlapping regions between active regions (top-10%) of a given cell line and regions with less activity (bottom90%) of another cell line (i.e., from NSC ChIP-STARR-seq or comparative ChIP-STARR-seq in ESCs). We also considered the active regions that were in the top-10% of both NSC ChIP-STARR-seq and comparative ChIP-STARR-seq as common highly active regions.

NCRE-gene interaction

To integrate gene expression data, NCREs were linked to their target genes using HiC data and distance to the closest TSS. Previously published HiC data from the germinal zone (GZ) at gestation weeks 17–18 of 3 human foetal brains were used¹⁷ to link NSC NCREs to target

genes. This data provides 10 kb resolution bins for gene loop interactions and 40 kb resolution for TADs. Moreover, to link NCREs identified in comparative ChIP-STARR-seq in ESCs to their target genes, we used published HiC data from H9 human ESCs with 25 kb resolution bins for both gene loop interactions and TADs³⁶.

NCREs were then intersected with pre-calculated significant interactions using intersectBed to define NCRE-gene interaction. NCRE-gene interactions located within the same TAD were considered for downstream analyses. The coordinates of HiC data were liftedOver to the GRCh38/hg38 human genome build before intersecting with regions. In addition, NCREs were annotated for target genes based on the distance to nearest TSS using the “annotatePeaks.pl” function of HOMER (version 4.11)⁷². Finally, all NCRE-gene interactions either using HiC or distance to closest TSS were considered for further analysis. In total, we found 28,514 target genes (coding and non-coding) for NSCs and 20,347 target genes for ESCs, of which 19,739 target genes are identified based on the closest TSS and the other target genes are based on HiC.

We only included protein-coding genes in our analysis. Gene expression levels were plotted and statistical comparison was performed using Wilcoxon signed rank test in R. The Online Mendelian inheritance in Man (OMIM) gene list (updated 21-04-2022) was downloaded from the OMIM database to further assist in annotation of target genes.

Gene ontology analysis

The web interface of Enrichr^18,19 (https://maayanlab.cloud/Enrichr/) was used for functional enrichment analysis using the default settings, and the whole genome was set as background. The geom_point function of ggplot2 (version 3.5.0) was used to visualize the GO enrichment.

Functional enrichment

NCREs were annotated for genome features using the “annotatePeaks.pl” function of HOMER (version 4.11)⁷². Motif enrichment analysis was performed using the “findMotifsGenome.pl” function of HOMER ⁷² with default setting and whole genome as background. Significant enriched motifs were selected based on p-value ≤ 0.01.

Sequence characteristics analysis

To determine whether different DNA sequence features distinguish defined regions between NSCs and ESCs, we considered the following features: (i) the non-coding essential regulation (ncER) score (https://github.com/TelentiLab/ncER_datasets/; updated 06-03-2019)²⁰; (ii) GC content, as determined by the GCcontent R packages based on BSgenome. Hsapiens.UCSC.hg38 (version 1.4.3); (iii) conservation score for each NCRE, as derived from the gscores R packages based on phastCons100way.UCSC.hg38 (version 3.7.2)²¹; (iv) Orion scores²²; (v) CADD scores²³; (vi) Haploinsufficiency scores²⁷ and (vii) probability of loss-of- function intolerance (pLI) score⁷³. The overlaps between DNA sequence features and regions were defined using intersectBed. The coordinates of data were converted to the GRCh38/hg38 human genome build by liftOver prior to intersections. The scores of the genes linked to the regions were plotted to show pLI distribution. Statistically significant differences between groups were determined using the Wilcoxon signed rank test in R.

Transposable element enrichment

The RepeatMask (GRCh38/hg38, updated 03-09-2021) was downloaded from the UCSC table browser and intersected to the regions. To determine enrichment of transposable elements in NCREs, we followed a strategy previously used for investigating active NCREs in human embryonic stem cells⁹. The number of overlaps of each type of repeat (n_overlaps) with all regions (n) was used to calculate the relative frequency (f_all = n_overlaps/n). Multiplication of the relative frequency with the number of regions (n_test, e.g., ESC, NSC etc.) in any tested group yields the expected frequency (E). This number was compared with the actual observed frequency in the subgroups (f_test = (n_overlap, test)/n_test = O) to calculate the observed versus expected ratio (O/E). We considered repeats with O/E > 2 as enriched. For the subsequent data interpretation, we only focused on transposable elements that were present multiple times (n_overlap > 15).

Epigenome profiling of NCREs

To find the enrichment of endogenous marks around NCREs in ESC and NSC, we collected various ChIP-seq data from previously published studies^9,54,74 and from the ENCODE project⁷⁵.

The FASTQ files of ChIP-seq data was re-analyzed along with ATAC-seq data (from this study) based on what we described in the “Data processing” section.

Heatmaps for differentially active regions across the endogenous marks were created using deeptools (version 3.4.3)⁶⁷. We first converted bam files of endogenous marks into bigwig format using the “bamCoverage” function with --binSize 10 and --normalizeUsing RPKM parameters. The bigwig files were then used to compute the number of reads across 6 kb centered on NCREs using “computeMatrix” with --referencePoint center, --upstream 3000 and --downstream 3000 parameters. Finally, the line plots and heatmaps were generated using the “plotProfile” and “plotHeatmap” functions, respectively.

Predicting NCRE activity using DeepSTARR

In our initial attempts to predict NCRE activity based on DNA sequence composition (Extended Data. 6a), we used the previously published DeepSTARR framework⁵⁰ As this model requires 249 bp input sequences, we first generated 249 bp regions centered on the NCRE regions. Subsequently, we applied the prediction function using the default settings of DeepSTARR, as instructed on the author’s Github repository (https://github.com/bernardo-de-almeida/DeepSTARR). Finally, we calculated the Pearson correlations between DeepSTARR output and ChIP-STARR-seq experimental results.

Construction of BRAIN-MAGNET

To generate BRAIN-MAGNET, we first removed from the analysis scaffold of 148,198 genomic regions assessed in ChIP-STARR-seq a small number of patch sequences resulting in a final number of 148,114 genomic sequences. These ranged from 500 to 1000 base pairs in length, and were augmented with their reverse complement to yield a total of 296,228 examples. A random 10% of the augmented sequences were held out for both validation and testing purposes. The DNA sequences were extracted from the human genome using bedtools getfasta (version 2.30.0) and their reverse complements were obtained using revseq of EMBOSS (version 6.6.0.0)⁷⁶. The input sequences were then converted into an one-hot encoded matrix of length 1000, with each nucleotide being represented by a binary vector (A = [1,0,0,0], C = [0,1,0,0], G = [0,0,1,0], T = [0,0,0,1]). Sequences shorter than 1000 bp were padded with zeros. The corresponding sequence activities served as outputs of the model. We subsequently used input and output in BRAIN-MAGNET, a single task convolutional neural network, to predict either activities of NCREs identified by NSC ChIP-STARR-seq or by comparative ChIP-STARR-seq in ESCs. The model is composed of three convolutional layers (filters = 128, 256, 512; size = 11, 9, 7), each followed by batch normalization, a ReLU non-linearity, and MaxPool (size = 2). After the convolutional layers, there are two fully connected layers, each with 1024 neurons, followed by batch normalization, a ReLU non- linearity, and dropout where the fraction is 0.4. The final layer is mapped to either NSC or comparative ChIP-STARR-seq outputs.

Hyperparameters were manually adjusted to optimize the performance on the validation set. The model was implemented and trained in Pytorch (version 1.13.1 with CUDA version 11.7)⁷⁷ using the Adam optimizer with a learning rate of 0.01, mean squared error (MSE) as the loss function, a batch size of 128, and early stopping criterion with a patience of 20 on 100 epochs. Model training, hyperparameter tuning, and performance evaluation were conducted on different sets of genomic regions.

The performance of the model was evaluated separately for NSC and ESC predictions on the held-out test sequences. We used the Pearson and Spearman correlation coefficient (R and Rho) across all bins for a quantitative genome-wide evaluation (Extended Data Fig. Sb). Furthermore the area under the precision-recall and receiver operating characteristic curve (PR and ROC; calculated using pr.curve and pr.roc functions from R package PRROC (version 1.3)⁷⁸ were used for NCRE classification in the various activity groups (top-10% vs. the other defined NCRE activity classes) (Extended Data Fig. 6c). To benchmark BRAIN- MAGNET to existing variant prioritization scoring methods, we first retrieved LINSIGHT⁴⁴, CADD²³, GERP⁴⁵ and ncER²⁰ scores for all nucleotides for which BRAIN-MAGNET provides cb scores. We subsequently calculated the Pearson correlation between BRAIN-MAGNET cb scores and the outputs of the other variant prioritization scoring methods (Extended Data Fig. 6d). To reduce the computational burden, for this analysis we used the previously used test set of NCREs, containing randomly selected 10% of our NCRE atlas. Furthermore, we similarly calculated and plotted the correlation between BRAIN-MAGNET cb scores and LINSIGHT, CADD, GERP and ncER score for all nucleotides of the functionally investigated motifs and mutations studied in this work (Extended Data Fig. 6e), and visualized all these scores as heatmaps for the individual studied loci.

Motif discovery using TF-MoDISco-lite

DeepExplainer⁴², a DeepSHAP implementation of DeepLIFT, was used to compute contribution (cb) scores for all nucleotides of the assessed NCREs, which reflects the individual contribution of each nucleotide to the NCRE activity in either NSCs or ESCs. 100 dinucleotide-shuffled versions of each input sequence were used as reference sequences (the fork from https://github.com/kundajelab/shap/blob/master/shap/explainers/deep/deep_pytorch.py). The obtained hypothetical importance scores were then multiplied by the one-hot encoded matrix of the sequences to derive the actual nucleotide cb scores, which were visualized using the ggseqlogo (version 0.2)⁷⁹. The cb scores of each NCRE can be visualized using the R-shinny app for which the code is available in our github repository: https://github.com/ruizhideng/NSC-ChIP-STARR-seq_paper.

The hypothetical importance scores, along with the corresponding one-hot encoded matrix of the sequences, were subjected to TF-MoDISco-lite (https://github.com/jmschrei/tfmodisco-lite), a more efficient and user-friendly version of TF-MoDISco⁴¹, to identify the most relevant motifs associated with changes of each nucleotide. The TF-MoDISco-lite algorithm in default settings was used to extract and cluster similar seqlets across all selected sequences, the identified seqlets were matched against JASPAR 2022 CORE vertebrate non-redundant database (https://jaspar.genereg.net/download/data/2022/CORE/JASPAR2022_CORE_non-redundant_pfms_jaspar.zip) by Tomtom (version 5.5.1)⁸⁰ within TF-MoDISco-lite⁴¹. To further facilitate the interpretation of the cb scores across the genome for downstream analysis, we first converted the negative cb scores to positive values allowing to include the impact of enriched motifs at sequences with negative cb scores. Subsequently, we calculated normalized cb scores taking either a normalisation of all NCREs (cb_all) or normalisation within each NCRE (cb_each) into account, and calculated percentile scores for these. All raw and normalized cb scores used throughout the study are (given the file size) available via a defined USCS browser session: https://genome.ucsc.edu/s/BarakatLab/BrainMagnet_NSC_ESC_cb_scoreshg38

GWAS analysis

To test the utility of BRAIN-MAGNET cb scores for the functional interpretation of GWAS loci, we made use of a recent MPRA data set provided by Guo et al, which tested 2,221 non-coding variants associated with 10 neuropsychiatric disorders by MPRA, with an average of 5 SNPs per disease risk loci being tested⁴⁶. Comparing the activity of matched risk versus reference allele pairs, the authors identified 892 SNPs that showed differential activity in their MPRA depending on the genotype of the given SNP, which the authors refer to as daSNVs. We retrieved the daSNVs and non-daSNVs from Supplementary Data 3 and 5 of Guo et al ⁴⁶ and intersected these with NCREs and cb scores from our data set. We subsequently plotted the distribution of cb scores (>60 cb_all) over the daSNV and non-daSNVs and additionally plotted the GERP, LINSIGHT, CADD and ncER scores for the same SNPs. Furthermore, we used the BRAIN-MAGNET cb scores to fine map two previously identified GWAS loci for which multiple SNPs in LD were assessed in the Guo et al study and present in multiple NCREs from our data set.

Genomic variant analysis

For the gnomAD variant analysis, the curated VCFs of gnomAD v4 were downloaded from https://gnomad.broadinstitute.org/downloads#v, which contains 909,084,110 short variants mapped to the GRCh38 genome build obtained from whole-genome sequencing of 76,215 individuals. bcftools (v1.8) was used to keep high-quality variants (“PASS” flag) and rare variants (AF < 0.1%). tabix (v1.9) was used to extract the variants overlapping our NCREs.

The total number of high-quality rare variants was calculated without chromosome Y and mitochondrial DNA (chromosome M), as these sequences were not covered in our NCRE analysis scaffold. The VISTA enhancers, ENCODE_cCREs, fiveUTR, threeUTR and codingExon regions (GRCh38/hg38, updated 14-02-2024) were downloaded from the UCSC table browser. DAE and nDAE regions were retrieved from Yousefi et al²⁶. The total number of filtered rare variants (MAE <0.1) genome wide (variant_genome) was normalized by the size of the genome (size_genome) and used to calculate the expected frequency (E = n_variant/size_ genome). The number of filtered rare variants located in the target regions (listed in Extended Data Fig. 5a) (variant_region) was normalized by the size of the target region (size_region) to calculate the observed frequency (O = variant_region/size_region). The O/E ratio was calculated by comparing the observed frequency over the expected frequency.

For the de novo variant analysis from the Genomics England 100,000 Genomes Project (GEL), de novo variants of GEL data version 18 (main-programme_v18_2023-12-21) were retrieved by Labkey from the de novo_flagged_variants table in the GRCh38 genome build. Next, high confidence de novo variants (stringent_filter = 1) were kept. Then, solved cases, labelled with “yes”, “unknown” and “partically” in the gmc_exit_questionnaire table, were removed. Finally, 4,415 individuals with unexplained neurology and neurodevelopmental disorders were kept, for which 13,276 de novo variants overlapped with our NSC NCREs. Also we analysed 4,558 genomes from individuals with a variety of other rare disease categories (including cardiovascular disorders, hearing an ear disorders, ultra-rare disorders, dermatological disorders, gastroenterological disorders, renal and urinary tract disorders, respiratory disorders, endocrine disorders, dysmorphic and congenital abnormality syndromes, growth disorders, skeletal disorders, ophthalmological disorders, metabolic disorders, hematological and immunological disorders and rheumatological disorders), harbouring 13,564 de novo variants overlapping NSC NCREs. The de novo variants were further annotated by the BRAIN-MAGNET cb scores, and the cb_each distribution of NCREs were plotted.

For the GEL analysis using variants overlapping with high-confidence motifs, we used 20 bp sequences with coordinates centred on the high-confidence motifs called by TF-MoDISco (q- values < 0.05) from the NSC NCRE category 4 and 5. An aggregation of gVCFs (AggV2, only available in version 10) comprising 78,195 genomes from GEL harbouring over 722 million variants was intersected with these high-confidence motif coordinates. Only variants marked with the “PASS” flag, indicating high quality, were retained, resulting in 11,788 variants overlapping 2,980 motifs. We next retained only variants in unsolved individuals (55,379 genomes, 9,787 variants) and subsequently individuals with neurodevelopmental disorders (10,536 genomes, 4,259 variants). Next we retained only variants with allele count (AC) <= 10 (2,700 genomes, 2,684 variants). Finally, only variants located in NCREs with target genes recorded as OMIM disease genes were kept, resulting in a final list of 705 variants from 824 genomes that overlapped 440 motifs.

For screening of an independent WGS data set, we first used 330 trio WGS from individuals affected by rare disorders from the Institute of Medical Genetics and Applied Genomics (IMGAG) database from Tübingen, Germany. Removing those individuals in which already a definitive molecular diagnosis was found (n=81) left 249 individuals, of which 60 harbored a candidate disease explaining variant. In those 249 individuals, a total of 16,437 de novo variants was found. Intersecting these variants with the NCRE atlas and filtering for variants with a BRAIN-MAGNET cb score >95th percentile resulted in 70 de novo variants that were filtered for potential clinical phenotype matches by HPO terms and manual curation. Finally, 3,971 singleton WGS of undiagnosed individuals were screened for rare deletions overlapping with NCREs that contained nucleotides with a cb score >95th percentile. Among the 12,172 deletions called in this data set, 1,584 were rare (filtered by <10 exact in-house calls and <50 larger in house deletion calls). For 117 of those deletions, the associated HPO terms of the target gene of the deleted NCRE matched the patient HPO terms. We then manually curated for matching clinical phenotypes, focusing on deletions <20 kb in size not containing additional coding sequences, resulting in the 8 NCRE deletions, of which 3 are potentially disease- explaining and 5 represent carrier deletions with currently no second coding hit identified.

Experimental validation of wild type and mutant NCREs in vitro

For experimental validation, we randomly selected 15 NCREs of which the target gene was expressed in NSCs. NCREs were amplified from genomic DNA and cloned into the STARR- seq reporter plasmid (kind gift of A.Stark)¹¹ as previously described⁹. Subsequently, site- directed mutagenesis was used to generate NCRE deletions or point mutations in the generated plasmids. NSCs were then transfected with wild type or mutant NCRE plasmids using polyethylenimine (PEI, Sigma) or Lipofectamine™ Stem Transfection Reagent (Thermo Scientific) respectively. Twenty-four hours post transfection cells were collected, stained with Hoechst dye and the NCRE activity was measured by FACS analysis (20,000 cells per sample). GFP-positive cells within the mCherry-positive population were quantified to assess NCRE activity compared to an empty STARR-seq vector control. Two independent transfection experiments were performed, each in duplicates. Statistical analysis was performed using a one-way ANOVA test followed by multiple comparison test (Fisher’s LSD test). Calculations were conducted in GraphPad Prism (version 8).

For CRISPRi experiments, we first generated a HEK293 cell line stably expressing ZIM3- KRAB-dCas9⁸¹ upon doxycycline induction, by transduction with pLX303-ZIM3-KRAB-dCas9 (a kind gift from Mikko Taipale, Addgene plasmid # 154472). gRNAs targeting the RAB7A NCRE were cloned into a pgRGFP plasmid (Addgene #82695⁸², a kind gift of Allan Mullen) and transfected into the ZIM3-KRAB-dCas9 expressing cells using polyethylenimine (PEI, Sigma). 48h post transfection cell were collected for RNA isolation (Tri reagent, Sigma). cDNA was generated with the iSCRIPT cDNA synthesis kit (Bio-Rad) and the change in expression of RAB7A was quantified by qRT-PCR (iTaq universal SYBR Green Supermix) (Sigma), performed in a CFX96RTS thermal cycler (Bio-Rad), as previously described⁵⁹. The ΔΔct method was used to calculate the fold change in gene expression, using TBP expression as housekeeping normalisation control. Statistical analysis was performed using one-way ANOVA test and calculations were conducted in GraphPad Prism (version 8).

Transgenic NCRE assays in zebrafish

Zebrafish (Danio rerio) were raised and maintained under standard conditions, as previously described²⁶. Adult and larval fish were kept on a 14h/10h light–dark cycle at 28°C. Larvae were kept in HEPES-buffered E3 medium. Media was refreshed daily and at 24 hours post fertilization (hpf) 0.003% 1-phenyl 2-thiourea (PTU) was added to prevent pigmentation. All zebrafish experiments were performed in compliance with Dutch animal welfare legislation and only larvae <5 days post fertilization (dpf) were used for experimentation. The wild type and mutant RAB7A NCRE used in the in vitro experiments were transferred by Gibson assembly between the AscI and PacI sites of an E1b-GFP-Tol2 enhancer assay plasmid (a kind gift from Ramon Birnbaum⁸³) containing an E1b minimal promoter followed by GFP, using the following transfer primers: Transfer_fw: 5’- AGATGGGCCCTCGGGTAGAGCATGCACCGG-3’ and Transfer_rv: 5’-TCGAGAGATCTTAATGGCCGAATTCGTCGA-3’. Constructs were injected into fertilized eggs at the 1 cell stage (25 pg per egg) using standard procedures, with or without 25 pg transposase mRNA. At least 50 eggs were injected per construct in at least two different injection experiments. GFP expression was observed and annotated at 24 hpf, using a fluorescent Leica M165FC stereomicroscope. Embryos representative of each category were subsequently imaged on a Leica SP5 intravital imaging setup with a 20 × /1.0 NA water-dipping lens. Images were analyzed using imageJ (FIJI).