Regulatory genomic circuitry of human disease loci by integrative epigenomics

1Genome-wide association studies (GWAS) have been successful in discovering more than 100,000 genomic loci that contain common single-nucleotide polymorphisms (SNPs) associated with complex traits and disease-related phenotypes, providing a very important starting point for the systematic investigation of the molecular mechanism of human disease,. However, the vast majority of these genetic associations remain devoid of any mechanistic hypothesis underlying their molecular and cellular functions, as more than 90% lie outside protein-coding exons and probably have non-coding roles in gene-regulatory regions with circuitry that remains unresolved,.

2Large-scale experimental mapping– and integration of histone modification marks and DNA accessibility have helped to annotate diverse classes of gene-regulatory annotations, including distal-acting and tissue-specific enhancers and proximal-acting and mostly constitutive promoters,. These maps help to elucidate the molecular basis of complex traits by revealing preferential localization (enrichment) of trait-associated genetic variants in tissue-specific gene-regulatory elements,,– and by fine-mapping possible causal genetic variants in enriched annotations,–. However, these maps also have limitations: they miss many disease-relevant tissues, have variable quality, and are prone to experimental noise and methodological variation between protocols, laboratories, antibody lots, reagents, batches, computational processing pipelines, software versions and integration pipelines. Moreover, consortia that require common marks across samples often exclude samples that miss some marks or marks that are missing in some samples, thus reducing biological space coverage, and often only profile few marks in many samples, or many marks in few samples owing to cost limitations.

3Here we overcome many of these limitations and present a new human epigenome reference, EpiMap (for epigenome integration across multiple annotation projects) (Fig. 1a). We inferred chromatin-state annotations that combine multiple marks, and a high-resolution enhancer annotation that combines DNA accessibility and multiple chromatin enhancer states. We grouped enhancers into modules that show common activity patterns, and inferred candidate upstream regulators and enriched functions of downstream genes for each module on the basis of regulatory motif and gene ontology enrichments. We also inferred enhancer target genes using a machine learning approach. We integrated this high-resolution gene-regulatory circuitry with genetic association results, revealing traits with epigenomic enrichments, and predicting causal variants and tissue-specific target genes. We distinguished unifactorial, multifactorial and polyfactorial traits on the basis of the diversity of their enriched tissues, and partitioned the loci of polyfactorial traits according to their overlap in distinct enriched tissues, thus revealing their distinct biological processes and disease comorbidity patterns. We also distinguished monotropic versus pleiotropic loci, and found that top-scoring loci frequently have multiple predicted driver variants, converging through diverse pleiotropy patterns involving multiple enhancers with a common target gene, multiple genes in a common tissue, or multiple genes in multiple tissues. Our results demonstrate the utility of dense, rich, multidimensional, high-resolution epigenomic and regulatory circuitry annotations for gene regulatory studies, complex trait investigation and studies of disease locus mechanism, resulting in unprecedented scale, scope and coverage of biological space and disease complexity.

4We uniformly processed 3,030 observed– genomic tracks across 859 biosamples (406 ENCODE, 425 Roadmap Epigenomics and 28 Genomics of Gene Regulation (GGR) samples) that span 18 epigenomic assays, and computationally imputed 14,952 tracks (Fig. 1b, Supplementary Fig. 1, Supplementary Table 1), which are available for download and interactive visualization at http://compbio.mit.edu/epimap.

5Our imputed tracks matched held-out observed tracks, both visually across randomly selected regions (Extended Data Fig. 1a, b) and quantitatively with more than 85% peak recovery and more than 75% average genome-wide correlation for punctate marks (59% of tracks) genome-wide (Supplementary Fig. 2). Imputation was robust even with few supporting datasets, and performed best when target datasets showed more than 50% average correlation to their ten nearest datasets, which held for 98% of single-assay samples (Supplementary Fig. 3).

6Imputed data also matched independent post-data freeze experiments, outperforming ‘average signal’ and ‘nearest track’ benchmarks (in practice knowable only after generating the target track) for 96% of punctate marks and 77% of broad marks both genome-wide and specifically focusing on rare events (Extended Data Fig. 1c–h).

7Disagreement between imputed and observed tracks helped to flag 138 potentially problematic datasets, which independently also showed markedly lower quality control scores (Supplementary Fig. 2a–c) and revealed potential sample or antibody swaps (Supplementary Figs. 4, 5), some of which were independently flagged by the data producers. Subtraction of the imputed track signal from the observed track signal revealed 13 experiments with potential antibody cross-reactivity or secondary specificity (Supplementary Figs. 6–8). From subsequent analyses, we removed the 138 flagged datasets and 442 tracks based solely on assay for transposase-accessible chromatin with high-throughput sequencing (ATAC-seq) or low-quality DNase-seq data, resulting in 2,850 observed and 14,510 imputed marks across 833 biosamples used in the remainder of this work.

8The resulting compendium of 833 high-quality reference epigenomes, grouped into 33 tissue categories, represents a major increase in biological space coverage, with 75% (624 of 833) of biosamples corresponding to new biological specimens. Observed and imputed data co-clustered, but imputed datasets better captured the continuity between and within different sample categories, and more clearly revealed sample-type groups, probably driven by both their cleaner signal tracks and their sheer number. Moreover, the distances between imputed datasets were less affected by technical covariates (Supplementary Fig. 9) and were more consistent with sample groupings (Supplementary Fig. 10).

9Hierarchical and two-dimensional embedding clusterings of multiple marks using both genome-wide and relevant-region-specific correlation patterns (Extended Data Fig. 2a–c) grouped biosamples first by life stage (adult versus embryonic) and sample type (complex tissues versus primary cells versus cell lines), and second by distinct groups of brain, blood, immune, stem cell, epithelial, stromal and endothelial biosamples within them. Active marks (histone H3 lysine 27 acetylation (H3K27ac), histone H3 lysine 4 monomethylation (H3K4me1) and H3K9ac) primarily grouped biosamples by differentiation lineage (blood, immune, spleen, thymus, epithelial, stromal and endothelial) and tissue type (lung, kidney, heart, muscle and brain), while repressive marks (H3K27me3 and H3K9me3) captured life stage (pluripotent, induced pluripotent stem cell derived, embryonic and adult) (Extended Data Fig. 2b, Supplementary Fig. 11), consistent with previous studies,. Donor sex was not a primary factor in sample grouping.

10We annotated genome-wide locations of 18 chromatin states,, in all 833 biosamples using combinations of histone modifications, including multiple types of enhancer, promoter, transcribed, bivalent and repressed regions (Extended Data Fig. 3a), using a mixture of observed and scaled imputed data and excluding the 138 flagged observed datasets (Methods). Genomic coverage and mark frequencies remained stable across biosamples for most states (Extended Data Fig. 3b–d), but biosamples with fewer observed datasets showed more heterochromatin and Polycomb-repressed states, consistent with our previously noted lower imputation accuracy for broad marks.

11We annotated 2.1 million high-resolution active-enhancer regions by intersecting five active-enhancer states with 3.6 million accessible DNA regions from 733 DNase-seq experiments. These covered 13% of the genome cumulatively and 0.8% on average for biosamples individually (Fig. 2a, Supplementary Fig. 13a–g), and represent a more than twofold increase relative to the ENCODE 2020 release (Extended Data Fig. 4). Clustering biosamples by sharing of active enhancers captured biologically meaningful groupings (Extended Data Fig. 5, Supplementary Fig. 12).

12For each high-resolution active-enhancer region, we defined H3K27ac-based local activity levels across 833 biosamples and used them to group enhancers into 300 enhancer modules (Fig. 2a, b, Extended Data Fig. 6a–c, Supplementary Figs. 14, 15), including 290 tissue-specific modules (1.8 million enhancers, 88% of enhancers cumulatively, active in 2% of biosamples on average) and 10 broadly active modules (251,079 enhancers, 12% of enhancers, active across 77% of sample categories on average).

13Enhancer modules showed substantial high-resolution, tissue-specific gene ontology enrichments for neighbouring genes (Extended Data Figs. 6, 7, Supplementary Fig. 16), including ion channels (for brain modules); camera-type eye development (eye); neural precursor cell proliferation (neurosphere); endothelial proliferation, hemidesmosomes and digit morphogenesis (endothelial, stromal and epithelial); and organ development and morphogenesis (embryonic).

14We predicted 3.3 million tissue-specific enhancer–gene links by combining epigenomic–transcriptional correlation and genomic proximity, each gene linked to 13 enhancers and each enhancer to 1.5 genes on average, at a median distance of 42,359 bp. Links were approximately sixfold more specific than enhancers, and sample-specific links spanned larger distances than constitutive links (Extended Data Fig. 8a–f). Our links outperformed previous linking approaches, using both gene-set enrichment metrics and curated gold-standard datasets (Methods, Extended Data Fig. 8g, h), and greatly expanded the biosamples with predicted links (from 127 to 833).

15We predicted upstream regulators for 273 modules (91%), implicating 1,175 motifs grouped into 160 motif archetypes (Extended Data Figs. 6, 7, Supplementary Fig. 17), including 152 tissue-specific motif archetypes (enriched in 6 modules on average) and 8 broadly enriched (enriched in 53 modules on average). Specific motifs include: GATA and SPI1 in the blood and immune samples; NEUROD2 and RFX4 in the brain and peripheral nervous system,; KLF4 for digestive tissues; and TEAD3 for the placenta, myosatellite and epithelial cells.

16Broadly enriched motifs revealed highly connected, combinatorially acting master regulators, including HNF1A in the liver, kidney and pancreas (with NR5A2); AP-1 (also known as JUN) or JDP2 in immune, bone and cancer samples; and TEAD3, paired alternately with MYF6 (myosatellite), TFAP2A (placenta) and AP-1 (stromal) (Extended Data Fig. 7d, e).

17Motif enrichments often partitioned tissue categories into subgroups specific for developmental stage and tissue type (Fig. 2c, Extended Data Fig. 7a–c), including heart into embryonic heart (NFIX and E2F1), aorta and arteries (SRF and PAX5), and heart chambers (MEF2D and ESRRG); brain into embryonic (NFIX and NEUROD2), adult brain (RFX2 and SOX10), and astrocytes (NFE2L2 and JDP2); and haematopoietic cells into natural killer cells (ETV2), B cells (NFKB2 and SPIB), and multipotent progenitors (GATA1 and NFE2L2).

18We next used our 2.1 million enhancer annotations and their tissue specificity to interpret genetic variants associated with complex traits,,. We compiled a compendium of 803 well-powered GWAS with 10 or more significant loci and over 10,000 cases (15% of 5,454 GWAS publications) that capture over 70,000 GWAS loci (63% of the NHGRI-EBI catalogue).

19We found 12,494 significant trait–tissue enrichments, enabling fine-mapping of candidate driver SNPs in over 12,956 loci (18%) from 226 traits in tissue-enriched enhancers (false discovery rate of <1%) (Extended Data Fig. 9a, Supplementary Figs. 18, 19). New biosamples captured the strongest GWAS enrichment in 65% of cases (146 of 226) and the only significant enrichment in 14% of cases (n = 32), and our annotations captured 15% more GWAS studies than DNase alone 226 versus 196 (Extended Data Figs. 9b, 10a–d).

20To capture common enrichments of similar biosamples, we also calculated trait enrichments for enhancer modules, resulting in approximately twofold fewer enriched traits (93 versus 226) and over 15-fold fewer SNPs in enriched annotations (Supplementary Fig. 20). Instead, using the hierarchical enhancer-sharing tree (Extended Data Fig. 2c) to reveal the appropriate tissue resolution of GWAS enrichment in comparisons of subtree-versus-parent enhancers captured similar numbers of traits (232 versus 226) and 57% more SNPs in enriched annotations (20,428 SNPs) (Extended Data Fig. 10e, Supplementary Fig. 18), representing an approximately fourfold increase from the 54 traits enriched in H3K27ac and the 58 traits enriched in H3K4me1 reported by the Roadmap Epigenomics project, although this increase is in part driven by the older GWAS catalog and the more stringent SNP pruning methodology used in Roadmap Epigenomics.

21Our epigenomic enrichments and enhancer–gene links yielded new biological insights on disease loci, with many compelling examples. For breast cancer GWAS, enriched in epithelial and cancer biosamples (Extended Data Fig. 9c), the highly localized rs17356907 genetic signal (P = 10−39, rank no. 12) localized precisely in a narrow epithelial and cancer enhancer nearest to USP44 but linked instead to NTN4, which is implicated in tumorigenesis and angiogenesis (Extended Data Fig. 9d). For depressed affect GWAS, maximally enriched in the cingulate gyrus (Extended Data Fig. 9e), the diffuse rs1261070 genetic signal (P = 6.0 × 10−12, rank no. 8) overlapped a broad set of enhancers between the CCDC68 and TCF4 promoters, all of which linked to TCF4, a transcription factor central to neurodevelopment, suggesting that multiple causal variants may contribute jointly to its dysregulation (Extended Data Fig. 9f). We have provided an interactive website for exploring more than 20,000 additional loci across more than 200 traits at http://compbio.mit.edu/epimap.

22We then studied trait–tissue, trait–trait and tissue–tissue epigenome GWAS co-enrichment patterns to gain insights into their complex interactions. First, we used the number of distinct tissue categories enriched in each trait (Extended Data Fig. 9a (left), Supplementary Data 1) to distinguish: 47 ‘unifactorial’ traits (21%) with most enriched nodes in only one tissue group (for example, QRS duration in heart, depressed affect in brain, proinsulin levels in endocrine) versus 179 ‘multifactorial’ traits (79%) enriched in seven tissue categories on average (for example, Alzheimer disease in immune cells and the brain; coronary artery disease in heart, stromal, endothelial, and adipose samples), of which 16 ‘polyfactorial’ traits (7%) enriched in 2 tissue categories on average (including waist-to-hip ratio in muscle, stromal, digestive, heart, and kidney samples).

23Second, we used trait co-enrichment patterns in the same tissues to cluster GWAS traits with similar properties. The resulting network (Fig. 3) showed a small number of densely connected communities of primarily unifactorial traits (for example, blood glucose in pancreas, heartbeat intervals in the heart, cholesterol in the liver, filtration in the kidney, immune traits in T cells and corpuscular volume in haematopoietic cells) with multifactorial connectors between them (for example, CAD between heart, endocrine and stromal; HDL and triglycerides between liver and adipose; lung function between lung, heart and digestive tissue; blood pressure between heart, endocrine, kidney and muscle and cell count between hematopoietic and digestive tissue) (Supplementary Figs. 21, 22). Many biologically meaningful similarities in this epigenomic co-enrichment-based network are missed by a network based on genetic overlap (517 edges, traits sharing 5% or more loci at a 10-kb resolution), which only captures 12% of epigenomic co-enrichment edges (332 of 2,787) (Supplementary Figs. 23–25).

24Third, we used co-enrichment properties of pairs of tissues in the same traits to distinguish ‘principal’ tissues (for example, immune cells, liver, heart, brain and adipose tissues) that showed consistently higher enrichments versus ‘partner’ tissues (for example, digestive, lung, muscle and epithelial tissues) for the same GWAS traits, suggesting that they have driver rather than auxiliary roles (Extended Data Fig. 11a). Specific principal–partner tissue pairs co-occurred more frequently than expected (Extended Data Fig. 11b), and revealed biologically meaningful traits where they probably co-act (Extended Data Fig. 11c), including: liver with adipose tissue (for cholesterol traits), with digestive tissue (for gamma-glutamyl transferase) and with blood cells (for C-reactive protein levels); and adipose tissue with endothelial cells (for waist-to-hip ratio), with heart tissue (for CAD) and with muscle tissue (for waist-to-hip ratio).

25We next used our epigenomic annotations to partition multifactorial trait SNPs into tissue-specific components, by studying functional and disease enrichments for distinct subsets of enhancer-overlapping SNPs in each enriched tissue (Fig. 4a–d, Extended Data Fig. 12, Supplementary Fig. 26).

26For example, the 375 CAD-associated SNPs lying in enriched tissue enhancers partitioned into: 235 heart-enhancer SNPs enriched in artery, cardiac and vessel morphogenesis; 179 adipose-enhancer SNPs in lipid homeostasis and axon guidance, consistent with adipose tissue innervation processes; 219 stromal-enhancer SNPs driven by cardiac fibroblasts and enriched in endothelin receptor signaling; and 112 embryonic stem cell-derived–muscle enhancer SNPs, enriched in vasculature, cardiac chamber and aorta development.

27These partitions also showed distinct co-associations (Fig. 4e). For example: heart, muscle and spleen enhancer CAD SNPs co-associated with high blood pressure and adipose and stromal enhancer CAD SNPs associated with waist-to-hip ratio and HDL cholesterol.

28Individual multifactorial trait loci included both single-tissue and multiple-tissue loci (Fig. 4f). Some CAD loci overlapped only heart enhancers (for example, EDNRA and TCF21), some almost fully adipose (for example, MIA3 and VAMP3), some lacked any enhancer overlaps (possibly acting at non-enhancer levels of regulation, or in uncaptured tissues or conditions), and many overlapped enhancers that were active in multiple tissues (for example, LDLR, APOE, SH2B3 and COL4A1), suggesting multiple mechanisms of action even at the single-locus level.

29For example, the stromal-only CAD-associated locus near the transcription facor TWIST1 (rs2107595, P = 1 × 10−24, rank no. 23) showed a strong stromal-specific signal and stromal-specific enhancer–gene links to TWIST1 (Fig. 4g). The heart-only CAD-associated 500-kb locus near EDNRA contains two separate associations, the transcription start site-centred rs6841581 (P = 5 × 10−24, rank 26) and the enhancer-centred rs4583018 (P = 8 × 10−15, rank 65), both in strong coronary artery enhancers and both linked to EDNRA through strong artery links, putatively reflecting multiple functional variants that converge on the same target gene in the same tissue (Fig. 4h).

30Even seemingly single-tissue loci sometimes showed pleiotropic signals: the 1-Mb rs3184504 locus (P = 5 × 10−30, rank no. 13; Fig. 4i) showed primarily coronary artery activity in a locus with a broad GWAS signal. However, in addition to coronary artery links, top enhancer-gene links also tied SH2B3, which regulates inflammation and cytokine signaling, to adipose tissue and linked PHETA1, which encodes an endocytic protein, to endothelial cells. These examples highlight that even individual loci may be pleiotropic, a property repeatedly found for many top-scoring loci.

31In this work, we presented a comprehensive map of the human epigenome, EpiMap, encompassing approximately 15,000 epigenomic tracks across 833 distinct biological samples that greatly expand the coverage of both embryonic and adult tissues and cells. We combined observed and imputed datasets across 18 epigenomic marks to jointly annotate and distinguish diverse classes of chromatin states, including enhancer, promoter, transcribed, repressed and quiescent regions. We extensively validated the high quality of our annotations and found that they outperformed stringent benchmarks, using both held-out and external experimental datasets for validation.

32We used this resource to assemble a comprehensive view of human genome circuitry across primary tissues, cells and cell lines, annotating 2.1 million high-resolution gene-regulatory regions; their activity patterns across 833 biosamples; their enriched regulatory motifs, motif combinations and putative upstream regulators that are responsible for their co-regulation; their enriched gene functions and biological pathways that they probably control; and their tissue-specific target genes. Our high-resolution enhancer annotations provide a highly concentrated view of the non-coding landscape, yielding many gene-regulatory insights but covering only 0.8% of the genome in each sample and only 13% total across all samples. Our linking revealed the high number of enhancers that control each gene and the high tissue specificity of long-range enhancer–gene links. Our upstream regulator analysis revealed a highly combinatorial and hierarchical view of gene regulation, with a small number of master regulators (for example, RFX2–RFX4, GRHL1, HNF1A and AP-1) interacting with diverse partners in different tissues to define tissue-specific gene-regulatory programs.

33Our work has also provided high-resolution molecular investigations of complex traits and human disease circuitry. We found statistically significant epigenomic enrichments for 232 GWAS traits implicating 20,000 SNPs in tissue-enriched enhancers, used trait and tissue co-enrichment patterns to annotate tissue partnerships and trait pleiotropy, and to partition disease SNPs into tissue-specific functional components. For individual GWAS loci, our work provides mechanistic insights at unprecedented scale. We have highlighted specific examples of GWAS investigations at varying levels of complexity, from the typically sought single-enhancer to single-gene, to multiple enhancers converging on a single target, to multiple genes and multiple tissues acting in pleiotropy in a single locus.

34Beyond the specific examples highlighted in our figures, we have also provided a rich interactive supplementary website (Supplementary Fig. 27) for our study (at http://compbio.mit.edu/epimap), enabling detailed interactive exploration of functional and motif enrichments of 300 enhancer modules; motif–tissue networks and enrichments; GWAS enrichments for 232 traits against our biosample tree; GWAS-enriched tissue enhancer SNP overlaps and target gene predictions; and 20,000 disease locus visualizations with putative driver SNPs, enhancers, tissues and tissue-specific target genes. These can enable the generation of detailed hypotheses for future experimental follow-up in countless studies of gene regulation and disease.

35Our collection also has several limitations: tissue samples are not at single-cell resolution; we do not consider donor genotype or phenotype; imputation may result in increased homogeneity and miss rare sample-specific events; and we still miss many tissues, environmental and stimulation conditions, and developmental stages.

36Our work enables many future studies: hierarchical and multi-resolution tree-based analyses of gene regulation and GWAS; machine learning-based gene circuitry and combinatorial regulatory motif analyses,; more sophisticated network analyses of our tissue–trait, trait–trait and tissue–tissue relationships; and guiding the experimental prioritization, methodological development and validation experiments, which can continue to further our understanding of gene regulation and human disease circuitry.

37We analysed 3,030 datasets, including 2,329 epigenomic chromatin immunoprecipitation followed by sequencing (ChIP–seq) datasets, 635 DNase-seq datasets and 66 ATAC-seq datasets from ENCODE at https://www.encodeproject.org/, released as of 24 September 2018. These marks include tier 1 assays: DNase-seq, H3K4me1, H3K4me3, H3K27ac, H3K36me3, H3K9me3 and H3K27me3; tier 2 assays: ATAC-seq, H3K9ac, H3K4me2, H2AFZ, H3K79me2 and H4K20me1; tier 3 assays: POLR2A, p300, CTCF, SMC3 and RAD21; and tier 4 histone marks: 16 non-imputed histone acetylation marks, 4 methylation marks (H3K9me2, H3K79me1, H3K9me1 and H3K23me2), H3.3 and H3T11ph. We assigned unique sample IDs to each unique combination of: extended biosample summary, donor, sex, age and life stage, wherever each attribute was available. We removed samples with genetic perturbations and kept only samples with appropriately matched ChIP–seq controls. We provide a metadata matrix including the mapping between ENCODE accessions and our unique sample IDs (Supplementary Table 1; also at http://compbio.mit.edu/epimap). We mapped the 111 Roadmap biosamples and the 16 ENCODE 2012 biosamples to any of our biosamples with overlapping dataset accessions if the accessions were used in the flagship Roadmap epigenomics analysis. This mapping assigned 25 samples to ENCODE 2012 and 184 samples to Roadmap 2015, some of which were merged multi-donor samples in Roadmap, out of the final 833 samples that passed quality control. These were merged into 16 and 111 tissue types, respectively, in the Roadmap 2015 publication.

38We downloaded one alignment file per replicate, prioritizing filtered alignments aligned with BWA in hg19 whenever possible. We uniformly processed the ChIP–seq and DNase-seq datasets according to the processing pipelines established by the Roadmap Epigenomics Consortium. In brief, we filtered out improperly paired and non-uniquely mapped reads, truncated reads to 36 bp, filtered out a blacklist of low complexity and artefact regions (ENCODE accession ENCSR636HFF), and filtered reads against a mappability track of uniquely mappable regions for 36-bp reads. Truncating read lengths inevitably missed some repetitive regions that the longer reads could have helped resolve, but helped to avoid potential biases from alignment differences, as over two-thirds of the datasets had read lengths of 36 bp or lower (Supplementary Fig. 9). We converted.bam files to tagAlign, used liftOver to map GRCh38 alignments to hg19, and pooled all experiments within each ID and assay combination. We subsampled the pooled ChIP–seq datasets to a maximum of 30 million reads and the DNase-seq and ATAC-seq datasets to a maximum of 50 million reads. We used the SPP peak caller to estimate fragment length. In cases with extremely low fragment length in the ATAC-seq and DNase-seq datasets we used the average fragment length (73 bp) from the average of the rest of the tracks. We generated −log10 P value signal tracks against matched whole cell extracts for both the ChIP–seq and the accessibility datasets using the MACS2 and the SPP peak caller and cross-correlation analysis to identify the proper fragment length as in the Roadmap analysis.

39We carried out epigenomic imputation on 859 unique biosamples using ChromImpute for a total of 10,778 imputed datasets over 13 tier 1 and tier 2 assays using predictors trained on all 35 epigenomic assays across 859 samples. We also imputed 4,345 datasets for the five DNA-associated factors, using only the 35 epigenomic assays as features to train predictors with ChromImpute. We provide all imputed and processed observed tracks along with track sets for the 833 quality controlled samples at https://epigenome.wustl.edu/epimap.

40For imputation quality control and validation, we compared observed tracks to imputed tracks when both were available (that is, when at least two original observed datasets were available for that biosample). We calculated all imputation quality control metrics from the original ChromImpute publication, including genome-wide correlation, imputed and observed peak recovery (%), and the area under the receiver-operator characteristic curve (AUC) for all pairs of imputed and observed tracks. In addition to the quantitative metrics, we visually inspected the epigenomic predictions as part of our quality control. We showed (Extended Data Fig. 1b) three dense and varied regions of different resolutions (25 kb, 200 kb and 1.5 Mb) for each of two randomly chosen samples containing both observed and imputed tracks for each assay. We calculated the epigenomic profile quality metrics normalized strand cross-correlation coefficient (NSC), relative strand cross-correlation coefficient (RSC) and read depth for all datasets and compared these to the imputation quality control metrics (see the tables in Supplementary Table 1). We flagged low-quality tracks by detecting the elbow in the ranked correlation metrics, which we calculated as the point where the change in correlation exceeded 5% of the correlation. Validation on external datasets was carried out on 51 experimental tracks across eight marks and assays from ENCODE after our data freeze, similarly subsampled to 30 million (marks) and 50 million reads (accessibility), which we remapped from GRCh38 to hg19 and evaluated on fully remapped 200-bp bins (90.1%) in chromosome 1 (Extended Data Fig. 1c, d). For the data homogeneity analysis, we restricted the data to only biosamples in each mark with both observed and imputed data (Supplementary Fig. 10).

41To systematically identify both potential sample or antibody swaps and poor-quality experiments, we computed the correlation of each observed experiment against all 10,734 imputed tracks for histone marks and assays (all imputed tracks before removing samples by quality control). We then calculated the average correlation among the top 10 most similar tracks to each observed track. We flagged potential antibody swaps by comparing the average correlation against samples of the putative mark against those computed for other marks. We fitted a multivariate linear model to each mark comparison, flagged datasets with residuals greater than 3 standard deviations of the average correlation and visually confirmed seven antibody swaps (six low-quality tracks). Similarly, we flagged potential sample swaps by comparing the correlation between imputed and observed tracks against the average correlation in the top 10 tracks in the same mark. We fitted a multivariate linear model and flagged datasets with residuals greater than 3 standard deviations of the residuals distribution. We report 19 potentially swapped samples, of which 5 were also flagged as low-quality tracks (Supplementary Fig. 8).

42In addition to genome-wide quality control of imputed tracks, we also focused on the specific differences between observed and imputed tracks. For each observed mark, we generated a genome-wide ‘delta’ track, computed as the difference in signal intensity between the observed and the imputed data, rescaling imputed tracks to match the signal intensity properties of the observed tracks, as the observed tracks showed a general bias for higher intensity. Some of these ‘delta’ tracks showed surprisingly high correlations with ‘primary’ tracks of non-putative marks, indicating potential secondary antibody reactivities. To flag these reactivities, we compared the average correlation of each of the delta tracks to the top 10 closest imputed tracks for each mark. As with antibody swaps, we fitted a multivariate linear model in each mark combination to flag outliers. We flagged 19 tracks and reported 13 after visual inspection as potential secondary reactivities or single replicate swaps (for example, in the case of DNase-seq) (Supplementary Figs. 7, 8). We noted that some cases showed clear difference tracks that do not match available antibodies, suggesting that the secondary reactivity is not a common mark in our compendium.

43To evaluate the similarity of imputed and observed tracks across samples, we calculated the pairwise genomic correlations between all pairs of imputed and observed signal tracks. We hierarchically clustered the imputed or observed correlation matrix of each individual mark using Ward’s method. We averaged all imputed matrices for the six main marks (H3K27ac, H3K4me1, H3K4me3, H3K36me3, H3K27me3 and H3K9me3) to create a fused correlation matrix, which we similarly clustered. We plotted the hierarchically clustered tree for the fused matrix alongside the metadata information for each biosample using the circlize R package.

44In addition, we calculated mark-specific Spearman correlations that were restricted to relevant features within all observed and imputed tracks per mark. We mapped each of the 13 marks to its top state by emission probability in the ChromHMM 25-state model and any other states with emission probability over 80%. For ATAC-seq, we used the same region list as DNase-seq. For each mark, we averaged and reduced each 25-bp signal track to any 200-bp regions that were labelled as one of the states associated with the mark in any of the 127 imputed Roadmap biosamples under the 25-state model,. We calculated the Spearman correlation between sets of these region-restricted mark signal tracks and generated similarity matrices across all datasets for a mark. Using these Spearman correlation matrices on all observed and imputed signal tracks, we computed UMAP dimensionality reductions for each mark and assay using with the uwot R package with the default parameters, except for n_neighbours = 250, min_dist = 0.25 and repulsion_strength = 0.25.

45We computed epigenomic annotations on 3,533 imputed and 1,465 observed datasets for 6 marks on 833 samples using ChromHMM with the fixed 18-state model from Roadmap with the same mnemonics and colours. We used observed data wherever possible, except in cases with no observed data or where observed data were removed in quality control. The table of the signal tracks used to calculate the annotations is available as Supplementary Table 2. The observed data were binarized from signal tracks with a −log10 P value signal cut-off of 2. To binarize the imputed data and facilitate comparison with the observed data, we established mark-specific binarization cut-offs. We first separately calculated the overall probability distributions of all imputed and observed tracks for each mark. Then, for each mark, we set the imputed binarization cut-off value to the value of the quantile that matched the quantile in the observed data for the −log10 P value > 2 cut-off. We used liftOver to map all 833 (after quality control) ChromHMM annotations to GRCh38, using a stringent reciprocal mapping strategy, ensuring that all resulting GRCh38 regions were also 200 bp and non-overlapping, and we have provided these alongside hg19 annotations and as track sets at https://epigenome.wustl.edu/epimap/.

46We define active enhancers as the intersection of DHS regions with enhancer annotations and high H3K27ac signal (average signal of >2 in the region containing the DHS ± 100 bp). We defined DHS regions from an index list of 3,591,898 DHS element consensus locations in GRCh38, determined from 733 DNase-seq experiments, that we mapped using liftOver to 3,568,912 hg19 locations. We intersected the hg19 regions with the 833 imputed enhancer annotations (states 7, 8, 9, 10, 11 and 15 in the 18-state model). This resulted in 2,842,995 regions with at least one enhancer annotation in any biosample. Finally, we intersected this matrix with the H3K27ac signal in the ±100-bp region that encompassed each DHS from the same tissue-specific imputed and observed datasets used to calculate the ChromHMM annotations. This procedure resulted in 2,356,914 active-enhancer regions. We created an equivalent promoter element region using the promoter annotations (states 1, 2, 3, 4 and 14 in the 18-state model). We noticed that several regions shared both enhancer and promoter annotations. As a conservative cut-off, we assigned all regions to either enhancers or promoters if over 75% of its active occurrences were labelled as that type of element (Supplementary Fig. 13). This final thresholding procedure yielded 2,069,090 enhancers, 204,104 promoters and 122,358 dyadic elements (neither specifically promoter nor enhancer). The matrices and enhancer locations are available at http://compbio.mit.edu/epimap.

47For all images using tissue group order, including ChromHMM tracks and module heat maps, groups were ordered alphabetically within six major groups: tissue or organs (adipose, bone, digestive, endocrine, heart, kidney, liver, lung, mesenchymal, muscle, myosatellite, pancreas, placenta and EEM, reproductive, smooth muscle and urinary), other primary cells (endothelial, epithelial and stromal), blood and immune (blood and T cell, HSC and B cell, lymphoblastoid, spleen and thymus), nervous system (brain, eye, neurosphere and PNS), stem (embryonic stem cell-derived, embryonic stem cell and induced pluripotent stem cell), and other (cancer and other).

48To define enhancer modules, we clustered the binary enhancer matrix defined by intersecting enhancer annotations with DHS regions and with the average centred and flanking (±100 bp) H3K27ac signal above a −log10 P value of 2 using the k-centroids algorithm with the Jaccard distance and the number of clusters set to k = 300. The average module contained 6,897 enhancers, and the largest module (enumerating constitutive elements) contained 93,554 enhancer regions. In all heat map plots of module centres (and associated enrichment figures), we diagonalized the matrix by ordering each column in the heat map (module centres) by the biosample that contributed the maximal signal. All columns that had a signal over 25% in more than 50% of rows were shown first. We used this diagonalization procedure for all diagonalized heat maps. We coloured each module by the tissue group that contained its maximal signal. Modules highlight sample groupings and organize according to cell type and tissue. Major groups were ordered alphabetically within six major groups and samples were ordered within groups according to Ward method’s clustering of the Jaccard distance of the module centres matrix. We performed enrichment on the module centres against the metadata of included samples (signal over 25%) by the hypergeometric test, and show enrichments with −log10P > 2 (Fig. 2b).

49We performed gene ontology enrichments on each enhancer module using GREAT v3.0.0 for the biological process, cellular component and molecular function ontologies. We analysed and visualized the results in the same manner as in the Roadmap core paper. We only considered enrichments of 2 or greater with a multiple testing-corrected P < 0.01. For Fig. 4c, we reduced the gene ontology enrichment by modules matrix to terms with a maximal −log10P > 4 that were enriched in less than 10% of modules. The full enrichment matrix is shown in Supplementary Fig. 16. As in the case of the diagonalized module centres, we labelled each term according to the module containing its maximal signal. We used a bag of words approach (as described in Roadmap) to pick 36 representative terms out of 865 total terms for Extended Data Fig. 6b, such that each tissue group has at least one term and the rest are representatively allocated across groups.

50We performed motif enrichment analysis across enhancer modules as described in the Roadmap paper,. In brief, we measured the enrichment of 1,690 motifs consisting of the JASPAR (2018) core non-redundant vertebrate motifs, the HOCOMOCO v11 human motif set and the SELEX motifs by Jolma et al.. We computed the enrichments for each of the 1,690 motifs relative to a joint DHS and intergenic background, additionally controlled by 100 shuffled motifs for each motif. We reported the motif with the highest enrichment in any module for each of the 286 previously identified motif archetypes. We only reported motifs with a maximum log2-transformed fold change of at least 1, resulting in 160 motif archetypes (corresponding to 1,175 total motifs), which we show with their position weight matrix (PWM) logos against all 300 modules in Extended Data Fig. 6c.

51We predicted enhancer–gene links for each biosample using the Pearson correlation between gene expression and the histone mark activity of nearby enhancers (within 1 Mb) for six marks (H3K27ac, H3K4me1, H3K4me2, H3K4me3 and H3K9ac). We precomputed correlations between all genes and nearby enhancers across the 304 biosamples with paired expression data. A negative set of correlations for each enhancer was computed using random genes in a different chromosome. We predicted links for each biosample and ChromHMM enhancer state separately (states E7, E8, E9, E10, E11 and E15). Predictions were made by training an XGBoost classifier on the positive set of all valid links against their paired negative links, using precomputed correlations and distance to the transcription start site as features, and keeping all links with a probability above 5/7 (ref. ).

52We validated enhancer–gene links using curated gold-standard data in CD34, GM12878, HeLa and K562 cells (Extended Data Fig. 8). We compared four sets of correlation-based predictions (alone or with H3K27ac and H3K4me1 activity, and with and without distance-based rescaling) against distance alone, enhancer–gene links from Roadmap, and H3K27ac correlation and/or activity times distance (calculated using EpiMap tracks and enhancers in compared epigenomes). For methods without a threshold value, such as distance alone, only the nearest or highest score gene was used for each as a cut-off value for F1. In addition, we created a gene ontology-based gold-standard set of links from gene ontology terms that were enriched within enhancer clusters by GREAT. For each gene ontology term per cluster, we added enhancer–gene links for enhancers within 1 Mb of at least two genes in the gene ontology term. Negative link sets were constructed by taking physical and expression quantitative trait locus (eQTL) negative link sets that were also not enriched by gene ontology.

53We pruned the NHGRI-EBI GWAS catalogue (downloaded from https://www.ebi.ac.uk/gwas/docs/file-downloads on 3 May 2019) using a greedy approach: within each trait + PMID combination, we ranked associations by their significance (P value) and added SNPs iteratively if they were not within 5 kb of previously added SNPs (vs. 1 Mb SNP pruning used in Roadmap Epigenomics, which would instead result in 28% fewer SNPs if used here). We also removed all associations in the HLA locus (for hg19: chr6: 29,691,116–33,054,976). This reduced the catalogue from 121,000 to 113,000 associations. Finally, we reduced the catalogue to 803 unique GWAS (from 5,454 GWAS) with an initial sample size of at least 10,000 cases or individuals (wherever cases and controls were not annotated) and 10 GWAS SNPs after pruning. This resulted in 71,379 lead SNPs, which landed in 35,573 unique genome intervals when we split the genome into 10,000-bp intervals. In enrichment analyses, we considered enhancers as intersecting with GWAS SNPs if the SNP was within 2.5 kb of the enhancer midpoint.

54We performed the hypergeometric test to evaluate GWAS enrichments on flat epigenomes and on modules. For these flat enrichments, we compared the number of captured SNPs for each GWAS and enhancer set combination, using the 2 × 2 contingency table of: number of GWAS SNPs captured by the set; the total number of GWAS catalog SNPs captured by the set; the total number of GWAS SNPs; and the total number of GWAS catalog SNPs. We corrected these hypergeometric P-values using the p.adjust function in R with the BH method. To estimate empirical FDRs, we created 100 additional shuffled GWAS catalogs by resampling the trait names without replacement to break the genotype-phenotype associations. We calculated the average number of significant combinations in shuffled catalogs, divided by the number of significant combinations in the real GWAS catalog. Rarefaction curves were calculated on the flat epigenome enrichments by iteratively adding the sample that was either significantly enriched or the maximal enrichment for the most remaining GWAS until all GWAS were accounted for (Extended Data Fig. 10c, d).

55We constructed a tree by hierarchically clustering the Jaccard similarity of the binary enhancer-by-epigenomes matrix using complete-linkage clustering. Then, for each node in the tree, we calculated its consensus epigenomic set, defined as the set of all enhancers present in all leaves of the subtree, such that each node’s set was a superset of that of its parent. We then performed hyper-geometric enrichment tests for captured SNPs for each GWAS in the set of consensus enhancers in the node of interest. These tests were performed in the same manner as for epigenome-based and module-based enrichments, and P-values were corrected using the p.adjust function in R with the BH method.

56For the CAD example, gene ontology terms were calculated using the nearest gene of each enhancer hit by a lead SNP. We pruned genes to expressed genes by calculating the average RNA-seq profiles for each tissue group and excluded genes that had log2 FPKM < 2 in the average RNA-seq of each sample’s group. Of 833 samples, 341 samples have matched RNA-seq, which we list in addition to releasing the processed data at http://compbio.mit.edu/epimap. We kept only the gene ontology terms that were significant in 25% or less of nodes, and report the top two gene ontology terms per node in Fig. 4d and all gene ontology terms in Supplementary Fig. 26.

57For locus investigations (in NTN4, TCF4, TWIST1, EDNRA, and SH2B3), we found the nearest active enhancer to each lead SNP in each node (within 2.5 kb), plotted the H3K27ac signal in the 2.1 million enhancers only, and (1) directly mapped links that originated at one of the enhancers near a lead SNP in the top three enriched epigenomes or (2) any links in the locus present in at least half of the samples in one of the selected tissue groups.

58We assigned each internal node in the tree to a unique tissue if over 50% of the leaves of the subtree came from that tissue and as ‘multiple’ if the subtree was not the majority of one tissue. We assigned tissue labels to 641 of 832 (77%) internal nodes where the majority of leaves corresponded to a single group. Using these assignments, we created a tissue by GWAS matrix by adding the −log10 P values for each tissue node set from all of the GWAS enrichments on the tree. We binarized this matrix and computed the Jaccard similarity across tissues to calculate a tissue similarity matrix. To assess the significance of tissue overlap, we compared each overlap value against the overlaps from 10,000 permuted enrichments. We collapsed each permuted matrix into a tissue by the GWAS matrix to compute the overlaps under the null. We performed the permutations for each tissue against other tissues by shuffling the enrichment P values on the node by the GWAS matrix. Specifically, we (1) binarized the enrichment matrix, (2) fixed the column of the group of interest, (3) permuted the remainder of the matrix, keeping its row and column marginals the same, and then (4) calculated the cosine distance between the permuted and the original matrix of enrichments.

59To evaluate the cross-GWAS similarity, we normalized the tissue by the GWAS matrix for each GWAS to obtain the proportion of significance attributed to each tissue for each GWAS (Supplementary Fig. 21). We reduced the matrix to the 232 significant GWAS with at least 10,000 cases or individuals (wherever cases and controls were not annotated) and 10 GWAS SNPs after pruning, and with FDR <1%. We created a GWAS–GWAS network using the cosine distance matrix as an adjacency matrix, keeping 2,787 links with a cosine distance of 0.25 or less and 215 connected GWAS. We used the Fruchterman–Reingold algorithm to lay out the graph. We used the tissue by the GWAS matrix to colour links according to the maximum tissue in the product between each pair of nodes and to colour nodes according to the maximal tissue for each node (Supplementary Fig. 22).

60To compare the epigenetic network to trait genetic similarity, we binned SNPs in the GWAS catalogue into 10-kb windows starting from the beginning of each chromosome. We counted the number of intersecting bins between two traits and kept any trait pairs with Jaccard similarity of at least 1%. To compare this to the epigenetic network, we plotted only links in the epigenetic network that coincided with any SNP-sharing GWAS pairs. In addition, we plotted the heat maps of the tree enrichments distance matrix and the genetic similarity matrix side by side, first organized by hierarchically clustering the enrichments matrix and then by clustering the genetic similarity matrix (Supplementary Figs. 23–25).

61Further information on research design is available in the Nature Research Reporting Summary linked to this paper.

62Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at 10.1038/s41586-020-03145-z.