
==== Front
Nucleic Acids Res
Nucleic Acids Res
nar
Nucleic Acids Research
0305-1048
1362-4962
Oxford University Press

39180401
10.1093/nar/gkae685
gkae685
AcademicSubjects/SCI00010
Narese/7
Narese/24
Computational Biology
Methods for constructing and evaluating consensus genomic interval sets
https://orcid.org/0009-0009-2762-9048
Rymuza Julia Department of Genome Sciences, School of Medicine, University of Virginia, Charlottesville, VA 22908, USA

Sun Yuchen Department of Genome Sciences, School of Medicine, University of Virginia, Charlottesville, VA 22908, USA
Department of Computer Science, School of Engineering, University of Virginia, Charlottesville, VA 22908, USA

https://orcid.org/0000-0002-1287-4931
Zheng Guangtao Department of Computer Science, School of Engineering, University of Virginia, Charlottesville, VA 22908, USA

https://orcid.org/0000-0002-7354-7213
LeRoy Nathan J Department of Genome Sciences, School of Medicine, University of Virginia, Charlottesville, VA 22908, USA
Department of Biomedical Engineering, School of Medicine, University of Virginia, Charlottesville, VA 22904, USA

Murach Maria Department of Genome Sciences, School of Medicine, University of Virginia, Charlottesville, VA 22908, USA
Department of Biochemistry and Molecular Genetics, School of Medicine, University of Virginia, Charlottesville, VA 22908, USA

Phan Neil Department of Genome Sciences, School of Medicine, University of Virginia, Charlottesville, VA 22908, USA
Department of Computer Science, School of Engineering, University of Virginia, Charlottesville, VA 22908, USA

https://orcid.org/0000-0001-9723-3246
Zhang Aidong Department of Computer Science, School of Engineering, University of Virginia, Charlottesville, VA 22908, USA
Department of Biomedical Engineering, School of Medicine, University of Virginia, Charlottesville, VA 22904, USA
School of Data Science, University of Virginia, Charlottesville, VA 22904, USA

https://orcid.org/0000-0001-5643-4068
Sheffield Nathan C Department of Genome Sciences, School of Medicine, University of Virginia, Charlottesville, VA 22908, USA
Department of Biomedical Engineering, School of Medicine, University of Virginia, Charlottesville, VA 22904, USA
Department of Biochemistry and Molecular Genetics, School of Medicine, University of Virginia, Charlottesville, VA 22908, USA
School of Data Science, University of Virginia, Charlottesville, VA 22904, USA
Department of Public Health Sciences, School of Medicine, University of Virginia, Charlottesville, VA 22908, USA
Child Health Research Center, School of Medicine, University of Virginia, Charlottesville, VA 22908, USA

To whom correspondence should be addressed. Tel: +1 434 924 7626; Email: nsheffield@virginia.edu
23 9 2024
24 8 2024
24 8 2024
52 17 1011910131
29 7 2024
05 7 2024
03 8 2023
© The Author(s) 2024. Published by Oxford University Press on behalf of Nucleic Acids Research.
2024
https://creativecommons.org/licenses/by/4.0/ This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

The amount of genomic region data continues to increase. Integrating across diverse genomic region sets requires consensus regions, which enable comparing regions across experiments, but also by necessity lose precision in region definitions. We require methods to assess this loss of precision and build optimal consensus region sets. Here, we introduce the concept of flexible intervals and propose three novel methods for building consensus region sets, or universes: a coverage cutoff method, a likelihood method, and a Hidden Markov Model. We then propose three novel measures for evaluating how well a proposed universe fits a collection of region sets: a base-level overlap score, a region boundary distance score, and a likelihood score. We apply our methods and evaluation approaches to several collections of region sets and show how these methods can be used to evaluate fit of universes and build optimal universes. We describe scenarios where the common approach of merging regions to create consensus leads to undesirable outcomes and provide principled alternatives that provide interoperability of interval data while minimizing loss of resolution.

Graphical Abstract

Graphical Abstract

National Human Genome Research Institute 10.13039/100000051 R01-HG012558 National Institute of General Medical Sciences 10.13039/100000057 R35-GM128636 NIH 10.13039/100000002
==== Body
pmcIntroduction

Advancements in high-throughput sequencing technologies have resulted in a vast amount of diverse epigenomic data that has given us tremendous insight into genome function. Epigenomic data are often summarized into genomic region sets stored in BED files. Through the work of hundreds of individual labs and projects such as ENCODE (1), the NCBI Gene Expression Omnibus (2) now contains almost 100 000 BED files (3). The volume of data has made integration challenging.

For an analysis that spans several genomic interval sets, one of the first steps is to define a consensus region set, or region universe, upon which the diverse sets can be interpreted (4–10). Such universe region sets have many common practical use cases. For example, they define genomic intervals for differential peak analysis (11); they form the regions of interest in a count matrix in single-cell epigenome analysis (11,12); they are used as a background for statistical region enrichment analysis (13–16); and they are a region vocabulary in vector representation approaches (17–20).

For many tools, the choice of universe is critical. It defines the features to which data will be projected. Currently, there are different ways of choosing a universe for analysis. Simple approaches include tiling the genome into fixed-size bins (12,15), or using intersection or union operations on a collection of region sets (11). Some methods have been developed to create better-fitting universes for specific downstream use cases (21,22). An alternative is to use a predefined universe from an external source; for example, the ENCODE consortium curated a registry of candidate cis-regulatory elements accessible through the SCREEN webserver (1), and the Ensemble Regulatory Build is a central, reusable source of regulatory region definitions (23). The choice of universe matters because universes can be a poor fit to data, and if a universe does not fit the data well, it can lead to incomplete or incorrect results (5,14,15). However, despite the importance of selecting a universe, it is often done ad hoc, and there are few approaches to assess the fit of a universe to a collection of region sets.

Here, we address these limitations by introducing novel concepts for constructing and evaluating region universes. First, we introduce the idea of flexible genomic intervals, which represent region boundaries by intervals instead of points, allowing us to summarize many fixed regions into fewer flexible regions without loss of information. Next, we propose three methods for constructing flexible region universes: a coverage cutoff universe, a maximum likelihood universe, and a Hidden Markov Model (HMM). Finally, we propose three methods to evaluate the fit of a universe to a collection of region sets: (i) the base-level F10-score; (ii) a Region Boundary Distance score (RBD); and (iii) a likelihood model score that assesses the likelihood that the proposed universe was drawn from the given distribution of region sets.

To assess our universes and evaluation methods, we compared our methods against alternatives and predefined universes. We show that flexible universes can capture information from complex data collections into one well-defined universe. Moreover, we show how our assessment metrics provide complementary measures of assessing universe fit and we prove the relevance of these measures. We show that the union universe has many downsides and propose the HMM universe as a generally useful approach for defining well-fit universes. To demonstrate how these universes could affect downstream analysis, we conclude with an application of region set enrichment analysis, where we show how the results are affected by choice of universe. Overall, our results demonstrate the importance of considering region universe and provide promising new tools to construct better-fitting universes for a variety of use cases.

Materials and methods

Overview

To integrate genomic interval data, we first require a consensus set of intervals, or a universe. We may select a predefined universe from an external source or define one from a collection of input region sets using a consensus algorithm (Figure 1A). With a universe in hand, we can then ‘tokenize’ the original regions (Figure 1B). Tokenization redefines them into universe regions, normalizing differences in region boundaries to transform similar regions into a single representation. The most basic tokenization method is simple interval overlap. This approach works well if a universe approximates the original data well; otherwise, this may result in loss of precision. A universe may not be a good fit to a collection for a variety of reasons (Figure 1C); for example, (i) a region can be shifted; (ii) two neighboring regions may be merged, making them indistinguishable in downstream analysis (5); (iii) a universe may omit important intervals, leading to loss of information (5); or (iv) a universe may contain extraneous regions that do not reflect genome coverage, adding noise and compute time (14). If a universe is a poor fit, it can affect downstream analysis negatively; for example, a differential accessibility analysis wouldn’t even test a locus that had been dropped from the universe. It could also miss a significant locus if it were merged with an abutting locus that lacked differential signal. Or, a motif analysis based on universe regions that had been shifted could miss enriched sequence features that were present in the part of the region left out of the universe.

Figure 1. Overview of the concept of a universe. (A) A consensus algorithm takes a region set collection \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\mathbb {R} = [\mathcal {R}_1, \mathcal {R}_2, ...]$\end{document} and builds a consensus representation called universe \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\mathcal {U}$\end{document}. (B) We ‘tokenize’ raw regions into universe \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\mathcal {U}$\end{document} by redefining them as universe regions, creating a more uniform collection. (C) Universes may poorly represent region sets by shifting, merging, dropping or adding extraneous regions.

To address these issues, we developed three new approaches for constructing a universe that is a good fit to the original data: first the ‘coverage cutoff universe’; second, the ‘maximum likelihood universe’; and finally, an ‘HMM universe’. To assess them, we also developed three universe fit metrics. Finally, we applied these to a variety of real datasets.

Methods for building optimal universes

The coverage cutoff universe

The simplest example of a universe built from the data is a ‘union universe’, in which a collection of region sets is merged. This method is often done for differential analysis of ATAC-Seq data (11). The union universe by definition covers all bases from the original collection; however, it can also lead to very large regions, particularly if the number of input region sets is large. Another simple alternative is using an intersection operation, which would include only bases covered in every region set in the collection, but this has the opposite problem: it leads to very sparse universes.

We reasoned that a hybrid approach may achieve a better result. First, we conceptualize a collection of region sets as a coverage signal track across all input region sets. Then, similar to a peak calling approach, we choose a cutoff x such that universe includes only positions with coverage greater or equal to x (Figure 2A). Setting the cutoff to one corresponds to a union universe, and setting the cutoff equal to the number of input region sets corresponds to an intersection universe. Setting a cutoff in between the two balances these extremes and provides a tunable parameter that may be adjusted depending on the needs of downstream tasks. We call the resulting universe a coverage cutoff (CC) universe. A principled approach to selecting the cutoff is to use a simple likelihood model that calculates the probability of appearing in a collection. With this model, we can calculate an optimal cutoff according to Eq. (1) (see Supplementary methods for details).

(1) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{equation*} x \ge \frac{S_c}{g}. \end{equation*}\end{document}

Figure 2. Different approaches to building universes. (A) Coverage-based universes are derived from the genome coverage of a collection of region sets. Examples include intersection \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\mathcal {U}_{int}$\end{document}, coverage cutoff \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\mathcal {U}_{CC}$\end{document}, and union universe \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\mathcal {U}_{union}$\end{document}. (B) A flexible region in contrast to fixed region can represent boundaries of many variable regions. (C) The flexible coverage cutoff (CCF) universe is based on coverage of the genome by a collection. It uses two cutoff values: the lower defines flexible boundaries and the upper defines the region core. (D) A collection of genomic region sets is aggregated, and region starts, core (overlap), and ends are counted, creating signal tracks. (E) Maximum likelihood universe is derived from three signal tracks. Using a likelihood model, we build a scoring matrix that assesses the probability of each position being a given part of a flexible region. Next, we find the most likely path, which represents the maximum likelihood universe. (F) The HMM universe treats signal tracks representing genome coverage by different parts of a region as emissions of hidden states that correspond to different parts of flexible regions.

Here, Sc is a sum of genome coverage by collection and g is the size of the genome.

We realized that, in a sense, the CC universe is a point estimate of a more complex distribution of possible universes. We reasoned that we may gain some insight by modeling the boundaries of the consensus regions as intervals, rather than points. To do this, we developed a new concept of a genomic region we call a flexible region. In contrast to fixed interval that is defined by two fixed points (start and end), a flexible interval is defined by four (startstart, startend, endstart and endend) (Figure 2B). A flexible interval can model many region variations into one well-defined flexible interval.

A simple approach to constructing a flexible region universe is to define two cutoff values instead of one: the looser represents the cutoff for boundaries, and the stricter for the region core (Figure 2C). This way, positions with coverage between those two points will be assigned to flexible region boundaries and positions with coverage higher than the second cutoff will be assigned to core of the region. Using this idea, we built a confidence interval around the optimal cutoff value, which extends the CC universe into the coverage cutoff flexible (CCF) universe.

Maximum likelihood universe

While flexible intervals more naturally represent collections of overlapping region sets than fixed intervals, we reasoned that they still suffer from the possibility of merging neighboring regions when collections are large. To address this issue, we need information about not just the coverage of regions in a collection, but also the region start and end positions. To compute this information, we developed a fast plane-sweep algorithm (see Supplementary methods for details). With this tool we can quickly and efficiently calculate three tracks representing aggregate start, end, and coverage values of a region set collection at base-pair resolution (Figure 2D). We reasoned that a model that could incorporate all of these signals may improve universe resolution.

We can conceptualize a flexible universe as a path through the genome that assigns either start, core, or end state to each position. Using a universe scoring model and optimization algorithm, we can build the best path through the genome (Figure 2E). As a scoring model, we next developed a complex likelihood model, an extension of the simple likelihood model introduced earlier for CC universe, which considers not just the coverage (core), but also in region start and end signal tracks. This model describes for each position the probability of it being a region start, core, or end. We thus build the maximum likelihood universe (LH) in 3 steps: (i) compute the three signal tracks; (ii) use a likelihood model to build a scoring matrix; (iii) find the maximum likelihood path through the genome (see Supplementary methods for details).

Hidden Markov model universe

The maximum likelihood universe provides a simple and principled model for optimal flexible universes. However, a disadvantage is that it provides no tunability, since the likelihood scores are determined purely from the data. We reasoned that this may lead to results depending on input collection. To address this, we sought a more tunable model using a Hidden Markov Model (HMM).

An HMM models a hidden processes using (i) a matrix of transition probabilities between hidden states and (ii) emission probabilities of observations from hidden states. In our model, there are three observed sequences: the number of starts, overlaps, and ends at a given position. The hidden variable corresponds to the different parts of the flexible segment (Figure 2F and Supplementary Figure S1). We can tune transition probabilities, which can be chosen in a way that will prevent unnecessary segmentation, and emission matrix, which describes the relationship between observations and hidden states (see Supplementary methods for details).

Methods for evaluating universe fit

Having developed several new approaches to construct universes, we next sought to evaluate these universes and compare them to other common approaches. Because the choice of universe can dramatically affect downstream analyses, it is important to choose a universe deliberately. However, there are no well-established methods for assessing universe fit to data. Furthermore, different analyses may be better served by a different types of universe, indicating that there really is no generally optimal universe, but the idea of what makes a ‘good’ universe depends on the downstream analysis. For example, a differential accessibility analysis should prioritize sensitivity over specificity; in this case, it is not a major problem if the universe includes many regions present in only a few samples, since the cost of extra comparisons is lower. On the other hand, the cost of excluding a region that could be a significant differential locus would be high. In contrast, a word-based deep learning task that trains a model with input dimensions equal to the size of the universe may elect a more specific universe, at the cost of discarding some regions that are present in a few of the samples, because otherwise the training could be intractable. Thus, the question of universe optimality depends on the use case, and therefore, we require methods of evaluating universe fit that can be tuned to a research question. This problem is similar to the comparison of two generic interval sets, for which several methods have been developed (24), with two key differences: first, we want to compare a universe region set not only to one other region set, but to a collection of them; and second, the question is not symmetric: it is generally more important that a universe not miss information (regions), even at the cost of some extra regions – and the desired level of asymmetry can vary. Therefore, we developed three methods for assessing universe fit to data: (i) a base-level overlap score; (ii) a region boundary distance score and (iii) the universe likelihood.

Base-level overlap score

Our first metric is based on base coverage. We consider the universe as a prediction of whether a given genomic position is present in a given region set from the collection. Treating each region set from the collection as a query, we can then conceptualize matches and mismatches as true positives (covered in both universe and query), false positives (covered in universe, but not in query), or false negatives (covered in query, but not in universe) (Figure 3A). This allows us to calculate common classification evaluation measures such as precision and recall. Precision counts the number of true positives, so a low precision indicates presence of unimportant positions in the universe; recall measures how much of the universe is in a query. To combine precision and recall, we use the F10-score, a weighted version of the traditional F-score that pays 10 more times attention to recall than precision. This asymmetry captures our goal to prioritize sensitivity over specificity: by prioritizing recall, we indicate that it’s better to have a few extra, noisy regions that to exclude something important. The F10-score results in values between 0 and 1, with a perfect fit approaching 1. An alternative approach to base-level overlap score would be Jaccard similarity, however it does not account for asymmetry of the comparison.

Figure 3. Different approaches to assess how well the universe represents the data. (A) The base-level overlap measure considers the universe as a prediction of a region set and based on that it calculates number of false positives (FP), true positives (TP), and false negatives (FN), and from that derives recall (R) and precision (P), which are combined into the F10-score. (B) The region boundary distance (RBD) score assesses how well a universe represents start and end positions, by calculating distance from region set to universe, and from universe to region set; those two metrics are combined into a region boundary score by calculating their reciprocal, weighted harmonic mean. (C) Likelihood assessment uses a likelihood model based on signal tracks representing genome coverage by different parts of a region to calculate universe likelihood as a combination of likelihoods of all three signals tracks. D) A complete analysis example comparing a collection of region sets against 4 proposed universes: \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\mathcal {U}_1$\end{document} a precise universe, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\mathcal {U}_2$\end{document} a sensitive universe, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\mathcal {U}_3$\end{document} a fragmented universe, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\mathcal {U}_4$\end{document} well-fit universe. For each universe, all 3 metrics are calculated. The F10-score and RBD score assess individual region sets. The final score for a collection is their average. In contrast, the likelihood is calculated directly for the whole collection.

Region boundary distance score

One disadvantage of the base overlap score is that it is unaware of region boundaries. A universe region that covers two abutting query regions would get a perfect score. This can be highly problematic in downstream applications; for example, in a differential analysis, lumping two distinct loci together could dilute differential signal. To address this, we sought a measure that would consider region starts and ends (Figure 3B). We calculate the distance between each boundary of each region in the query and the closest corresponding boundary in the universe. Universes with boundaries that are near the query boundaries would have shorter distances, indicating better universe fit. However, highly fragmented universes with many unnecessary boundaries would have very small distances from query to universe. To account for this, we also calculate the inverse distance: from boundaries in universe to the nearest boundaries in the query. Finally, we combine those two metrics into a region boundary distance score (RBD) by taking their reciprocal, weighted harmonic means. With this score we describe the universe’s ability to conserve information about starts and ends, with a score of one representing a perfect representation of boundary locations. However, we do not incorporate any information about collection coverage in this score.

For fixed universes, the start and end point are well-defined, however for flexible regions they are intervals. Therefore, for flexible universes, we modify the RBD score to set distance equal to zero if a boundary query region is inside the universe’s boundary interval.

Universe likelihood

Finally, we sought a metric that incorporates both information about region boundaries as well as genome coverage. We propose here a universe likelihood score (Figure 3C). We first calculate three signal tracks representing genome coverage by start, core, and end of the regions in the collection. Then for each signal track we make a separate model, which results in three separate models for different parts of a region. Each of these models describes the probability of a given position being a given part of a region, depending on the signal strength (see Supplementary methods for details). That results in a complex, probabilistic description of a region set collection. Next, we use this model to calculate the likelihood of the universe, which we can compare between universes. We make two versions of likelihood calculations, one suited for fixed universes and one for flexible universes. We use log likelihood, so our values range from minus infinity (low) to zero (high). Finally, to increase interpretability, we normalize the scores by subtracting the likelihood of an empty universe (one that contains no regions). Thus, a positive final score reflects a given universe that is more likely than having no regions at all, while a negative score means the universe is less likely than the empty universe.

To accommodate flexible universes, we adapted the likelihood score by calculating the boundary likelihood as if the whole flexible interval could contain a boundary position, rather than a single fixed point (see Supplementary methods for details).

Assessing region sets collections

Having developed three assessment methods, we can use them to compare competing universes to assess which universe is the best fit for a collection of region sets. We do this by computing the scores for each universe and comparing among universes (Figure 3D). The scores assess different aspects of the universe fit; the F10-score promotes sensitive universes over specific ones, RBD score penalizes sparse universes, and likelihood provides complex universe assessment. Although the likelihood score incorporates information about boundaries as well as how well the universe covers the collection, it can penalize sensitive universes because it is not intentionally biased toward asymmetry the way the previous scores are. Thus, by computing all these scores, we reason that we get a complete picture that can guide decisions for selecting a universe for a collection of region sets.

Evaluation on real data

Next, we developed an evaluation strategy to test our universe building and assessment methods on real data. We assembled five diverse collections of region sets representing different biological problems (Figure 4A): (i) CTCF ChIP small, a small random collection of CTCF region sets (n = 40) from the ENCODE database (1); (ii) CTCF ChIP large, CTCF ChIP-seq datasets (n = 877); (iii) TF ChIP, ChIP-seq experiments for diverse transcription factors (TFs) (n = 8503); (iv) B-LCL ATAC, a small set of ATAC-seq files from B-Lymphoblastoid cell lines (B-LCL; n = 400) from ChIP-Atlas (25); (v) a Random ATAC, random ATAC-seq results (n = 5000). These datasets vary in data type, collection size, and level of heterogeneity of input regions across region sets.

Figure 4. Overview of evaluation approach. (A) Five collections representing different biological problems used for assessment. (B) For each collection, we compared it to five data-driven universes and three predefined universes. The data-driven universes are tailored to the input collection, but the predefined universes do not vary by collection. (C) We assessed the fit of each universe to each collection using our three assessment methods.

We also assembled universes to assess. First, we obtained three universes that do not depend on analyzed data, which we call predefined universes: (i) the tiles universe, which bins the genome into non-overlapping 1000 bp tiles; (ii) the SCREEN universe, which consists of predefined cis-regulatory elements from ENCODE (1); and (iii) the Regulatory Build (RB) universe, consisting of pre-defined regulatory elements from Ensembl (23). In addition to these three external universes, we also built 5 data-driven universes that are specific to each region set collection. These include: (i) the union universe; (ii) the CC universe; (iii) the CCF universe; (iv) the LH universe; and (v) the Hidden Markov Model (HMM) universe. This led to 28 universes and 40 pairwise comparisons of universe-to-collection (Figure 4B). For each comparison, we computed our three assessment methods (Figure 4C). This gives us a comprehensive evaluation of both externally sourced and data-driven universes, tested on diverse query region set collections.

Results

Data overview

To explore the differences in our five region set collections, we first computed general coverage statistics (Supplementary Table S1). The smallest region set, B-LCL ATAC, contains ≈700 000 regions and covers 0.2% of the genome, whereas the largest, TF ChIP, contains ≈1.5 billion regions and covers 91% of genome. We also observed that the ATAC-seq collections have smaller regions on average than the ChIP-seq collections.

Universe overview

The universes also have very different characteristics, with some requiring additional filtering based on region likelihood and size (see Supplementary methods for details, Supplementary Figures S2 and S3). For example, for the CTCF ChIP large collection, the eight universes have different levels of precision and fragmentation (Figure 5A). The universes differed in average region size, number of regions, and percent of genome covered (Supplementary Figure S4, Supplementary Table S2). Having assembled the universes, we next computed our three assessment methods.

Figure 5. Universes overview and results of base-level overlap score. (A) Example of universes assessed for the CTCF ChIP large collection, including the 3 constant external universes, and 5 data-driven universes built from the input collection. (B) Different universes represent genome coverage by the collection to a different extent, example from the Random ATAC collection. Collection \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\mathbb {R}$\end{document} consists of many different files, which are represented by the core signal track. Regions in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\mathcal {R}_1, \mathcal {R}_2, \mathcal {R}_3$\end{document} are best represented by CC, CCF and LH universes in terms of overlap. (C) Precision and recall distribution for each collection and universes assessment. (D) Average F10-score for each collection and universes assessment.

Assessment 1: base-level overlap F10-score

Region sets in a collection can differ widely; our first assessment method assesses fit by quantifying the degree of overlap between each region set and the universe. In example data from the Random ATAC collection, we observe that some universes cover many bases present in only few of the collection’s region sets, while other universes are more stringent (Figure 5B). To assess this globally, we first computed precision and recall for each comparison (Figure 5C). We observed that the tiling universe and union universe both have perfect recall for all tested collections, consistent with how these universes are constructed; the tiling universe covers the entire (mappable) genome, and the union universe by definition covers every base contained in the collection. In contrast, recall is lower for the more stringent data-driven universes; the CC, CCF and LH universes exclude positions with low coverage, especially for large collections; the HMM universe has higher recall in general, indicating that it contains most positions covered by the collection. Finally, the lowest recall scores are assigned to the external universes, SCREEN and RB, which is consistent with these universes being built from other data sources. This highlights an advantage of building bespoke universes tailored to a collection: recall is superior. On the precision side, the worst performer overall is the tiles universe, consistent with many tiles in the universe that do not reflect coverage in the collection. In contrast, SCREEN and RB had generally higher precision, especially for ChIP-Seq collections.

In general, data-driven universes tend to represent collections well with good precision and recall. This is most apparent for the B-LCL ATAC collection for which data-driven universes are much better than predefined. For likelihood universes (CC, CCF, and LH universes) built from large ChIP-Seq collections (CTCF ChIP large and TF ChIP) we observe the worst recall among data-driven universes. This is the consequence of Eq. (1), from which we observed that, for ChIP-Seq collections, the optimal cutoff value is higher. On the other hand, both the HMM universe and the union universe have high recall and low precision. In general, we see that universes with high recall have lower precision, reflecting the delicate balance between including complete information in universes without adding too much noise.

To propose a balance between precision and recall, we next calculated the F10-score, which assigns more weight to recall than precision (Figure 5D). For predefined universes, we see that the tiles universe scores well for ChIP-seq collections, which cover more of the genome, but poorly for ATAC-seq collections, which have lower coverage. In general, both RB and SCREEN are outperformed by data-driven universes, especially for the B-LCL ATAC collection, for which they contain too much noise. In general, for ChIP-Seq collections, the union universe is the best for these metrics, consistent with the weighting we chose that gives 10 times the weight to recall. Interestingly, for the TF ChIP collection, the HMM universe outperforms likelihood universes; on the other hand, for ATAC-Seq collections, the CCF universe outperforms both union and HMM universes. Overall, we conclude that computing precision, recall, and F10 provide useful insight into assessing universe fit. They provide a way to quantify the advantages of a data-driven universe and assess how much information is lost by an external universe.

Assessment 2: region boundary distance score

Next, we sought to address the major weakness we see in the base-level overlap score: that it does not consider region boundaries. Assessing boundaries is important because of how it affects downstream analysis. If two regulatory elements with distinct behavior are merged into a single region in a universe, then all downstream analyses will essentially evaluate the average of the two signals. But different universes have different sensitivities to boundary points (Figure 6A); for example, anecdotally, the union universe is not sensitive at all: it contains few boundaries, especially for larger collections. It will clearly merge together many neighboring regions, even if they have distinct patterns across input sets. The CC and CCF universes are more sensitive but still miss out on many boundaries for bigger collections; the LH universe is very sensitive to boundaries, but has other weaknesses (it tends to exclude positions with low coverage); all of those problem are solved with HMM universe, which is very sensitive and also is able to represent regions with low coverage. To assess boundaries globally, we turned to the region boundary distance score. First, for each region set, for each region, for each boundary, we calculated the distance to the nearest corresponding universe boundary (See Methods; Figure 6B). For all collections, the distance from collection to RB universe was very high. Similarly, the union universe performs very poorly in this metric, particularly for larger collections, consistent with intuition that the union regions lose boundary precision as the number of regions increases. We also computed the inverse: distances from query to universe (Figure 6C). Interestingly, for ATAC-Seq collections, we observe a small distance from query to universe but a high distance from universe to query. This indicates that all universes have many boundaries that are not present in the raw data, but at the same time boundaries present in the queries are well-reflected by universes.

Figure 6. Results of region boundary distance score. (A) Different universes represent region boundaries to a different extent. Three signal tracks provide summarized description of the whole collection. Both LH and HMM universe are most sensitive to region boundaries. (B) Distribution of median distances from query to universe. (C) Distribution of median of distance from universe to query. (D) Average RBD score for each collection and universe comparison.

To summarize both distance directions, we calculated their reciprocal, weighted harmonic mean, the Region Boundary Distance score (RBD) (Figure 6D). The average RBD score shows that the RB universe is a very poor fit to all collections, likely because it contains few, large regions. On the other hand, tiles universe has similar scores for all collections, which is good for big ChIP-Seq collections compared to other universes, but bad for ATAC-Seq collections. Interestingly, the SCREEN universe seems to be a good fit for all collections, and the best fit for ATAC-Seq collections, even outperforming the data-driven ones for this metric. Among data-driven universes, the HMM universe performs the best, with LH in second place for all collections except B-LCL ATAC. However, for this collection all data-driven universes have similar scores. Coverage-based universes (CC, CCF) perform well for ATAC-Seq collections, but not for large ChIP-Seq collections (CTCF ChIP large and TF ChIP) compared to other data-driven universes. As expected, the RBD score reflects the poor performance of the union universe for large ChIP-Seq collections; the merging leads to poor reflection of interval boundaries.

Assessment 3: likelihood calculation

Finally, we calculated the likelihood score for each comparison (Figure 7A). In general, predefined universes are a worse fit to the collections than an empty universe, with exception of SCREEN for CTCF ChIP large and TF ChIP collections. Among data-driven universes, the union universe performs very poorly, achieving negative scores for all collections. As expected, CC, CCF, and LH universes outperform the empty universe for all collections, as these were designed to optimize likelihood in some way. The HMM universe performs well overall; however, it is worse than empty universe for CTCF ChIP large and TF ChIP collections. A more detailed look reveals that the low HMM likelihood scores in these scenarios are driven by region coverage tracks, not boundaries, suggesting that for these collections, our current HMM parameterization may yield a universe with too much noise (see Materials and Methods; Supplementary Figure S5).

Figure 7. Results of universe likelihood and comparison between fixed and flexible scores. (A) Likelihood of each universe given collection. (B) Change of RBD score when we account for flexibility. (C) Change of likelihood, when we account for flexibility.

Comparing flexible to fixed universes

So far, our assessments have not taken into account that some universes can be flexible. We believe that flexible universes provide several advantages, and sought to assess them. Flexible intervals don’t quite fit into the standard 3-column BED format; however, they can be stored using optional fields of an extended BED format. In this approach, sequence name, start, and end of the flexible region are represented by the first three columns, and thickStart and thickEnd columns hold the information about end of the flexible start and start of the flexible end. To assess this, we applied our flexible-aware version of the RBD score. We observed that RBD score improves for all flexible universes (Figure 7B). The change is less significant for CCF universe; however for LH and HMM universes, the new score is close to one, with the HMM universe performing slightly better.

We computed a version of the likelihood score that considers universe flexibility. This score shows a significant improvement; for all flexible universes, all collection scores change by an order of magnitude (Figure 7C). Likelihood values that consider flexibility are similar for all universes; however, the HMM universe performs slightly worse for large ChIP-Seq collections (CTCF ChIP large, TF ChIP). This reflects that, unlike the CCF and LH universes, the HMM universe does not explicitly optimize likelihood. A more detailed look into likelihood showed that although the HMM universe has the best likelihood of cores of the regions, it performs less well for boundary positions (see Methods; Supplementary Figure S6).

Assessing across metrics

Since each metric assesses different aspects of universe fit, considering them independently limits the scope of assessment. For a holistic view, we summarized the scores of each metric into normalized heatmaps, allowing comparison within and across metrics (Figure 8). We observed several informative cross-score patterns: First, the F10-score is the most consistent metric across universes, indicating that all these universes cover the collections to similar extent (Figure 8A). In contrast, the RBD indicates more variation in how well universes represent region boundaries (Figure 8B). This is consistent with our intuition that matching boundaries is a more difficult task, since it requires ensuring large regions are split well. The disparity between F10-scores and RBD score demonstrates their combined utility. For example, the union universe has perhaps the overall best F10-scores across universes, but has low RBD score. Inversely, the SCREEN universe has the best RBD score for the B-LCL ATAC collection, but is significantly worse than any data-driven universe for F10-score. In likelihood scores, the CC, CCF, LH universes outperformed other universes to similar extent, while tiles, RB and union were universally poor fits for all universes (Figure 8C). The HMM has much better RBD score than CC, but a worse likelihood, reflecting that likelihood considers more than boundary positions. Additionally, the likelihood is stricter in boundary assessments: while for RBD score we take a median of actual values, for likelihood we use probability of a given position being a boundary. Comparison of flexible and fixed versions of RBD score and likelihood highlights the value of flexible regions (Figure 8D, E); RBD score for LH and HMM improved significantly after accounting for flexibility.

Figure 8. Results of universe comparison using different scores. (A) Row normalized F10-score of each universe given collection. (B) Row normalized RBD score of each universe given collection. (C) Row normalized likelihood of each universe given collection. (D) Row normalized flexible version of RBD score of each universe given collection. (E) Row normalized flexible version of likelihood of each universe given collection.

Application on downstream analysis

To demonstrate how universe affects downstream analyses and how our universe building and evaluation methods can be applied, we performed a region set enrichment analysis using LOLA, a tool for statistical region enrichment analysis (15). The goal is to take some demo region sets and then use them to search a database of region sets to find similar regions, and explore how the choice of universe affects the results. We constructed two different experiments. For our first experiment, we used the Random ATAC collection as a database. We used three predefined universes – tiles, RB and SCREEN – as well as data-driven universes built from the Random ATAC and B-LCL ATAC collections. To see how the universes based on rare cell types perform, we also added data-driven universes built from fifteen Glia ATAC-Seq files. We queried the database with fifteen files that were not present in the database and were not used for universe construction, representing three different cell types: five files from A549, five from B-LCL, and five from Glia samples. For our second experiment, we used 623 ChIP-Seq files representing different TFs to build a database and data-driven universes. We queried the database with thirty files representing different TFs: 10 files from EZH2, 10 files from POLR2A, 10 files from YY1. For both experiments, we assessed performance of the region enrichment analysis with R-precision (rPPV), which measures the precision based on the top R results, where R is equal to the number of correct results in the whole database.

Our results demonstrate clear impact of the universe on analysis performance (Figure 9A). In the first experiment, the data-driven universes were the best performers for the specific questions; Glia ATAC queries (gray dots) performed best under the tailored Glia data-driven universes, followed by the Random ATAC data-driven universes, and performed poorly with any pre-defined universes or with the B-LCL data-driven universes. The B-LCL queries (yellow dots) also performed best with the tailored B-LCL-driven universes, and performed reasonably well with predifined or Random ATAC data-driven universes, but poorly with the Glia data-driven universe. Finally, A549 samples performed equally well with predefined or Random ATAC data-driven universes, and poorly with the universes built on the other data types. This shows that using universes based on a specific cell type increase the performance of downstream analysis for this type. Our second experiment, based on TF data, shows a clear difference between the union universe and other more complex data driven universes. We also see an imbalance among the predefined universes, with the SCREEN universe outperforming tiles and RB in general. In conclusion, for ATAC-Seq data, rPPV is higher for data driven universes than predefined ones; for TF data rPPV is higher for more complex data-driven universes (CC, CCF, ML, HMM) than union universes. Overall, these results demonstrate that choosing the right universe is a complex task. In our experiments, the more complex data-driven universes (CC, CCF, LH and HMM) are always either comparable or superior to simpler data-driven or pre-defined universes, although the exact performance depends on the both the initial data and the downstream task. Most importantly, this analysis suggests that the assessment methods we present can be helpful in choosing the right universe.

Figure 9. Results of downstream enrichment analysis depends on universe. (A) R-precision (rPPV) of query files depending on universe. (B) Median of different universe, depending on the collection used for their construction.

Discussion

Many integrative epigenome analyses require the data to be defined on a set of consensus regions, or universe. This universe is critical for analysis because it determines the precision of regions assessed. Our experiments highlight how different universes can have different levels of fit to collections and may therefore be useful for different tasks. Despite the importance of this choice, few approaches have been developed to aid analysts in building appropriate universes or assessing the fit of an existing universe. In this study, we have addressed these issues by presenting several novel methods to build universes from collections of region sets, as well as new ways to assess the fit of a universe to a collection of region sets. We also introduced the concept of flexible segments, and proposed several methods for constructing universes that can use either traditional fixed boundaries or flexible interval boundaries.

In general, data-driven universes outperformed predefined ones. However, the data-driven universes also have a weakness: by definition, they change with the underlying collection, and therefore cannot support an integrative analysis that spans collections. If results need to be compared across collection, then a shared universe is required. There are a few options for analysis in this case, all of which are facilitated by our work: First, a custom data-driven universe that spans all included collections can be built. Since we have described several ways to create well-performing data-driven universes, it would be easy to just design a bespoke universe for a given comparative analysis. However, this may not always be possible or convenient, and at some point, an external universe may be preferred. In this case, a trade-off is required: fit of the universe must be sacrificed to increase interoperability with other collections. Our assessment methods now provide a principled way to assess this trade-off and inform research decisions.

Among the non-data-driven universes we tested, SCREEN performed well for the collections we tested. However, there are almost certainly other collections or use cases for which the tiles, RB, or other external universes would be a better choice. Along these lines, we propose that our methods for building universes can be used in the future to create predefined universes from large collections, thereby creating even better global universes that can be re-used for integrative analysis. In the future, a centralized repository of universes, built using different methods and for different target use cases, could be a useful resource; a given collection of region sets could be represented into different universes based on the balance of fit and need for integration.

Based on our results, we propose that the union universe, though widely used, does not represent ChIP-Seq data well, particularly for large collections. Instead, we propose the HMM universe as a good all-around option that solves many of the issues with simpler universes. It has the highest sensitivity to boundary positions and good recall. It also provides adjustable parameters; by setting emission and transition probabilities, users may adjust model sensitivity and keep it consistent across collections. Still, the final choice of universe should consider the needs of downstream analysis. Our results show that assessing universe fit is a complex question, with many features to optimize. Even with our assessment metrics, it is difficult to claim an optimal universe for a given collection; instead, the answer depends on the downstream analysis priorities. For example, in general, the CCF and LH universe represent properties of a whole collection well, but they exclude infrequent regions; thus, they may be useful for NLP analysis of the genome but could lead to losing information about rare cell types in single cell analysis. In contrast, the union universe by definition covers all bases found in the region set collection, and therefore has a great F10-score; however, it also merges regions, which is reflected by poor RBD score and likelihood scores, indicating that it would not be a good fit for an application that requires high region resolution. Therefore, multiple perspectives must be considered for a holistic assessment of universe fit.

One advantage of our new universe construction methods is that they naturally create flexible universes. Flexible universes are a new way of summarizing information from large collections with less loss of information. We showed that proposed approaches of making flexible universes improve results over inflexible universes. We see flexible regions as a powerful new concept that can modify our current way of thinking about universes. Furthermore, we expect that using them for differential peak analysis, statistical region enrichment analysis and NLP approaches has potential to improve results. Flexible regions will become more useful as we and others develop the necessary tooling to work with them; for example, we will require tokenization methods that can project a traditional region set into a flexible universe quickly and accurately, while considering the universe flexibility.

In conclusion, this research provides new concepts, methods, and insight that will help researchers to determine the best analysis path for many types of genomic region analysis.

Supplementary Material

gkae685_Supplemental_File

Data availability

Software is available at https://github.com/databio/geniml.

Supplementary data

Supplementary Data are available at NAR Online.

Funding

National Human Genome Research Institute [R01-HG012558]; National Institute of General Medical Sciences [R35-GM128636]. Funding for open access charge: NIH.

Conflict of interest statement. N.C.S. is a consultant for InVitro Cell Research, LLC.
==== Refs
References

1. Moore J.E. , PurcaroM.J., PrattH.E., EpsteinC.B., ShoreshN., AdrianJ., KawliT., DavisC.A., DobinA., KaulR.et al . Expanded encyclopaedias of DNA elements in the human and mouse genomes. Nature. 2020; 583 :699–710.32728249
2. Barrett T. , WilhiteS.E., LedouxP., EvangelistaC., KimI.F., TomashevskyM., MarshallK.A., PhillippyK.H., ShermanP.M., HolkoM.et al . NCBI GEO: Archive for functional genomics data sets—update. Nucleic Acids Res. 2013; 41 :D991–D995.23193258
3. Xue B. , KhoroshevskyiO., GomezR.A., SheffieldN.C. Opportunities and challenges in sharing and reusing genomic interval data. Front. Genet. 2023; 14 :1155809.37020996
4. Kruczyk M. , UmerH.M., EnrothS., KomorowskiJ. Peak finder metaserver - a novel application for finding peaks in ChIP-seq data. BMC Bioinformatics. 2013; 14 :280.24059758
5. Lun A.T.L. , SmythG.K. De novo detection of differentially bound regions for ChIP-seq data using peaks and windows: Controlling error rates correctly. Nucleic Acids Res. 2014; 42 :e95.24852250
6. Lun A.T.L. , SmythG.K. Csaw: a bioconductor package for differential binding analysis of ChIP-seq data using sliding windows. Nucleic Acids Res. 2015; 44 :e45.26578583
7. Smith J.P. , SheffieldN.C. Analytical approaches for ATAC-seq data analysis. Curr. Protoc. Hum. Genet. 2020; 106 :e101.32543102
8. Fan H. , AtiyaH.I., WangY., PisanicT.R., WangT.-H., ShihI.-M., FoyK.K., FrisbieL., BuckanovichR.J., ChomiakA.A.et al . Epigenomic reprogramming toward mesenchymal-epithelial transition in ovarian-cancer-associated mesenchymal stem cells drives metastasis. Cell Rep. 2020; 33 :108473.33296650
9. Smith J.P. , CorcesM.R., XuJ., ReuterV.P., ChangH.Y., SheffieldN.C. PEPATAC: an optimized pipeline for ATAC-seq data analysis with serial alignments. NAR Genom. Bioinform. 2021; 3 :lqab101.34859208
10. Miller H.E. , MontemayorD., AbdulJ., VinesA., LevyS.A., HartonoS.R., SharmaK., FrostB., ChédinF., BishopA.J.R. Quality-controlled r-loop meta-analysis reveals the characteristics of r-loop consensus regions. Nucleic Acids Res. 2022; 50 :7260–7286.35758606
11. Yan F. , PowellD.R., CurtisD.J., WongN.C. From reads to insight: a hitchhiker’s guide to ATAC-seq data analysis. Genome Biol. 2020; 21 :22.32014034
12. Chen H. , LareauC., AndreaniT., VinyardM.E., GarciaS.P., ClementK., Andrade-NavarroM.A., BuenrostroJ.D., PinelloL. Assessment of computational methods for the analysis of single-cell ATAC-seq data. Genome Biol. 2019; 20 :241.31739806
13. Simovski B. , KanduriC., GundersenS., TitovD., DomanskaD., BockC., Bossini-CastilloL., ChikinaM., FavorovA., LayerR.M.et al . Coloc-stats: a unified web interface to perform colocalization analysis of genomic features. Nucleic Acids Res. 2018; 46 :W186–W193.29873782
14. Kanduri C. , BockC., GundersenS., HovigE., SandveG.K. Colocalization analyses of genomic elements: approaches, recommendations and challenges. Bioinformatics. 2019; 35 :1615–1624.30307532
15. Sheffield N.C. , BockC. LOLA: Enrichment analysis for genomic region sets and regulatory elements in R and Bioconductor. Bioinformatics. 2016; 32 :587–589.26508757
16. Nagraj V.P. , MageeN.E., SheffieldN.C. LOLAweb: a containerized web server for interactive genomic locus overlap enrichment analysis. Nucleic Acids Res. 2018; 46 :W194–W199.29878235
17. Gharavi E. , GuA., ZhengG., SmithJ.P., ZhangA., BrownD.E., SheffieldN.C. Embeddings of genomic region sets capture rich biological associations in low dimensions. Bioinformatics. 2021; 37 :4299–4306.34156475
18. Gharavi E. , LeRoyN.J., ZhengG., ZhangA., BrownD.E., SheffieldN.C. Joint representation learning for retrieval and annotation of genomic interval sets. Bioengineering. 2024; 11 :263.38534537
19. Zheng G. , RymuzaJ., GharaviE., LeRoyN.J., ZhangA., BrownD.E., SheffieldN.C. Methods for evaluating unsupervised vector representations of genomic regions. NAR Genom. Bioinform. 2024; 6 :lqae086.39131817
20. LeRoy N.J. , SmithJ.P., ZhengG., RymuzaJ., GharaviE., ZhangA., BrownD.E., SheffieldN.C. Fast clustering and cell-type annotation of scATACdata with pre-trained embeddings. NAR Genom. Bioinform. 2024; 6 :lqae073.38974799
21. Samb R. , KhadraouiK., BelleauP., DeschênesA., Lakhal-ChaiebL., DroitA. Using informative multinomial-dirichlet prior in a t-mixture with reversible jump estimation of nucleosome positions for genome-wide profiling. Stat.Appl. Genet. Mol. Biol. 2015; 14 :517–532.26656614
22. Granja J.M. , CorcesM.R., PierceS.E., BagdatliS.T., ChoudhryH., ChangH.Y., GreenleafW.J. ArchR is a scalable software package for integrative single-cell chromatin accessibility analysis. Nat. Genet. 2021; 53 :403–411.33633365
23. Zerbino D.R. , WilderS.P., JohnsonN., JuettemannT., FlicekP.R. The Ensembl Regulatory Build. Genome Biol. 2015; 16 :56.25887522
24. Chikina M.D. , TroyanskayaO.G. An effective statistical evaluation of ChIPseq dataset similarity. Bioinformatics. 2012; 28 :607–613.22262674
25. Zou Z. , OhtaT., MiuraF., OkiS. ChIP-atlas 2021 update: a data-mining suite for exploring epigenomic landscapes by fully integrating ChIP-seq, ATAC-seq and bisulfite-seq data. Nucleic Acids Res. 2022; 50 :W175–W182.35325188
