 Methodology article
 Open Access
 Published:
JOA: Joint Overlap Analysis of multiple genomic interval sets
BMC Bioinformatics volume 20, Article number: 121 (2019)
Abstract
Background
Nextgeneration sequencing (NGS) technologies have produced large volumes of genomic data. One common operation on heterogeneous genomic data is genomic interval intersection. Most of the existing tools impose restrictions such as not allowing nested intervals or requiring intervals to be sorted when finding overlaps in two or more interval sets.
Results
We proposed segment tree (ST) and indexed segment tree forest (ISTF) based solutions for intersection of multiple genomic interval sets in parallel. We developed these methods as a tool, Joint Overlap Analysis (JOA), which takes n interval sets and finds overlapping intervals with no constraints on the given intervals. The proposed indexed segment tree forest is a novel composite data structure, which leverages on indexing and natural binning of a segment tree. We also presented construction and search algorithms for this novel data structure. We compared JOA ST and JOA ISTF with each other, and with other interval intersection tools for verification of its correctness and for showing that it attains comparable execution times.
Conclusions
We implemented JOA in Java using the fork/join framework which speeds up parallel processing by taking advantage of all available processor cores. We compared JOA ST with JOA ISTF and showed that segment tree and indexed segment tree forest methods are comparable with each other in terms of execution time and memory usage. We also carried out execution time comparison analysis for JOA and other tools and demonstrated that JOA has comparable execution time and is able to further reduce its running time by using more processors per node. JOA can be run using its GUI or as a command line tool. JOA is available with source code at https://github.com/burcakotlu/JOA/. A user manual is provided at https://joa.readthedocs.org
Background
Genomic interval intersection is a major component of analysis pipelines for Next Generation Sequencing (NGS) technologies such as RNASeq, ChIPSeq, and exome sequencing. There are various existing tools that perform interval intersection [1–5] and many other genomic analyses. UCSC Genome Browser has been continuously improved since its first launch. Lately, the Data Integrator feature is released in UCSC Genome Browser, which allows users to combine and extract data from multiple tracks (up to 5 tracks), simultaneously [1]. BEDTools is a toolset developed for genomics analysis tasks such as comparison, manipulation, and annotation of genomic features in BAM, BED, GFF and VCF formats [2]. Pybedtools extends upon BEDTools and provides a Python interface for BEDTools [3].
BEDOPS is a highly scalable and easilyparallelizable genome analysis toolkit, which enables tasks to be easily split by chromosome for distributing wholegenome analyses across a computational cluster [4]. BEDOPS requires sorted interval set files before intersection. GROK utilizes region algebra for genomic region operations and provides R, Python, LUA, command line interfaces, and C++ API extension [5]. NCList uses a novel data structure named Nested Containment List (NCList) for interval databases [6]. NCList keeps intervals in lists such that intervals in each list have containment relation. However, no list fully contains one another, which implies that if lists are ordered on their start positions, then they are also sorted on their end positions. GLANET stands for genomic loci annotation and enrichment tool and provides joint overlap analysis of 2 to 3 tracks at most [7].
Layer et al. propose a novel parallel “slicethensweep” algorithm for nway interval set intersection. Their algorithm requires the input intervals not to be contained fully in one another [8]. Tabix indexes tabdelimited files and converts a sequential access file into a random access file [9].
In this paper, we generalize the problem of genomic interval intersection (finding common overlapping intervals) from 2 or 3 interval sets to n interval sets. We implement segment tree based construction and search algorithms in parallel for each interval set and chromosome. We accept n interval sets in Browser Extensible Data (BED) format and we apply a divide and conquer algorithm for finding common overlapping intervals. The intervals in a single set may overlap, may contain one another, and may be unsorted. We divide the interval sets into two sets, recursively, until one or two interval sets remain as base cases. In case of two interval sets, we construct a segment tree from one of the interval sets and we then search the intervals of the other interval set as query intervals on the constructed segment tree. In case of one interval set, we return all of its intervals as overlapping intervals. Results coming from the base cases are combined as follows: one of the results is used as query intervals, a segment tree is built on the other result, query intervals are searched on the segment tree, and overlapping intervals are returned. This process is repeated recursively until all of the n interval sets are processed.
We also propose a composite data structure, indexed segment tree forest which leverages on natural binning of segment trees and indexing the segment tree nodes at a certain depth. The same divide and conquer algorithm is also applied on the indexed segment tree forest (ISTF). However, this time when a segment tree is constructed, we proceed by converting this segment tree to an indexed segment tree forest by cutting this tree at a certain depth and indexing the segment tree nodes at that cutoff level. Our aim is to reduce the search space and search on the shorter indexed segment tree forest rather than one tall segment tree and eventually to reduce the query time. The overview of the proposed methodology is depicted in Fig. 1.
Results
Execution time comparison
We compared JOA with other tools using their latest available versions such as GROK v1.1.1, BEDTools v2.27.1, and BEDOPS 2.4.35. To compare the tools, we have used real and semisynthetic datasets. As an additional case study, we compared JOA segment tree and JOA indexed segment tree forest using 141 ENCODE Dnase hypersensitive sites. Details of these datasets can be found in Availability of data and material.
During comparisons, JOA is run both with segment tree (ST) and indexed segment tree forest (ISTF) options, GROK is run through its python API and intersectionL method is utilized, BEDTools is called by its intersect a b utility, and for BEDOPS, –intersect set operation is used. BEDTools is run with its “multiinter” parameter when there are more than 2 input files. BEDTools gives a number to each input file starting at 1, when run with the “multiinter” parameter and produces a comma separated list of the input file numbers that are jointly overlapping in the fifth column of its output file. Therefore, from BEDTools’ output file, we counted the number of rows that has all the file numbers separated by comma in the fifth column and calculated the number of jointly overlapping intervals.
All of the runs are carried out on Centos 6.6, Intel Xeon Gold 6132 (skylake), 2.6GHz ×28 processors (14 cores/socket). All the provided execution times in Tables 1, 2, 3 and 4 are averaged over 50 runs.
There exists other tools and formats such as BGT, GQT, SeqArray and BGEN [10–13] that requires joint overlap analysis. BGT is a new, compact file format for efficiently storing and querying wholegenome genotypes of tens to hundreds of thousands of samples. BGT search alleles/variants in the input bgt format file within a given site or list of sites [10]. Genotype Query Tools (GQT) is a new indexing strategy that expedites analyses of genome variation datasets in VCF format based on sample genotypes, phenotypes and relationships [11]. SeqArray is a storageefficient highperformance new data format for WGS variant calls which comes with SeqArray software package including key R functions [12]. BGEN is another novel binary file format for imputed genotype and haplotype data which are supported by many tools[13]. BGT, SeqArray and BGEN introduce their own new data format and carry out analysis on their new data formats except GQT, which carries out analysis on the VCF format. However, JOA does not support these data formats but designed for text based formats including “.bed” and “.pk’ formats. Therefore, we did not compare JOA against them.
Execution time comparison of tools using benchmarks datasets
JOA ST and ISTF are run for both (1 node, 1 processor per node) and (1 node, 8 processors per node) settings. Average execution times are listed in Table 1. Running times are comparable and usage of more processors per node decreased the running time of JOA because of parallel implementation of ST and ISTF algorithms. For “5 small files”, BEDTools “multiinter” parameter usage increased its running time as it has a different implementation which is also reflected in its output. Number of overlapping intervals found for each benchmark dataset is listed in Table 5. JOA and BEDTools found exactly the same number of overlapping intervals whereas GROK and BEDOPS found different number of overlapping intervals. This might be because of a different design, implementation or a case which is handled differently such as nested intervals.
Execution time comparison of tools using semisynthetic datasets
To compare the tools on a controlled dataset, we uniformly sampled 500bp long intervals from the human genome and generated bed format input. In the first scenario, we kept the number of intervals constant at 100,000 in each file and we increased the number of files by twofold from 2 to 512. In the second scenario, we kept the number of files constant at 2 and we increased the number of intervals in each file by twofold from 1,000,000 to 16,000,000.
In the firs scenario, we could not run GROK for more than 256 input files since calling a python function with more than 256 arguments was not allowed. Also, in the first scenario, BEDTools is called with the “multiinter” parameter when there are more than two input files.
Average execution times for both scenarios are provided in Tables 2 and 3. In the first scenario, when the number of input files is greater than 2, BEDTools average execution times are extremely high because of the “multiinter” option usage. In both scenarios, JOA has reduced execution times when 2 processors per node (ppn) is used instead of 1 ppn.
The number of overlapping intervals found for the first scenario and the second scenario are provided in Tables 6 and 7, respectively. JOA and BEDTools found the same number of overlapping intervals for the semisynthetic datasets as in the case of benchmark datasets. However, GROK and BEDOPS found different number of overlapping intervals.
Execution time comparison of JOA ST versus ISTF using ENCODE data
To show the parallel processing advantage of JOA, we designed and ran an additional case study for JOA using 141 interval sets of Dnase hypersensitive sites. We found and supplied all jointly overlapping intervals for these 141 interval sets in Additional file 1: Table S1. This additional case study verified that JOA can handle a significantly larger number of interval files and JOA is able to attain better running times with the help of its parallel processing ability as shown in Table 4.
JOA segment tree versus indexed segment tree forest detailed execution time comparison
We used the last four benchmarks (2.3M, 2.3M), (6.4M, 6.4M), (9.2M, 9.2M), and (6.4M, 9.2M) from benchmarks datasets in Table 1. We compared JOA ST versus JOA ISTF with respect to read, construction and query running times which are averaged over 50 runs. We show that ST and ISTF perform comparably with respect to read, construction and query running times as it is shown in Fig. 2.
JOA memory usage
Segment tree requires \(\mathcal {O}(nlog{n})\) storage for n given intervals. Since indexed segment tree forest is based on indexing of segment tree nodes at a certain depth, it has similar memory requirements. JOA is implemented in Java with fork/join framework which speeds up parallel processing by using all available processor cores. However, parallel processing requires loading all the data into memory. Also, segment tree and indexed segment tree forest construction requires memory allocation of many temporary data structures which can be only deallocated automatically by the Java garbage collector.
We calculated the memory usage of JOA ST and ISTF for the semisynthetic datasets. The memory usage in megabytes (MBs) are presented in Additional file 2: Tables S2 and S3. All the memory usage in the first and second scenarios are as expected; however, for the 2 interval sets of 8M intervals each, and for the 2 interval sets of 16M intervals in the second scenario, JOA’s memory consumption is high. This can be the due to the loading all interval sets into memory, the additional temporarily used data structures, or the Java garbage collector. As long as Java Virtual Machine has enough memory, garbage collector may not take action to deallocate memory and follow manually written memory deallocation statements in the code.
Discussion
Parallel implementation of JOA ST and JOA ISTF algorithms require all the input intervals to be loaded into the memory. Moreover, JAVA garbage collector deallocates memory automatically without taking manually written memory deallocations into account. Therefore, although JOA ST and JOA ISTF have comparable and sometimes lower execution times, their memory footprint may be high.
Regarding JOA ISTF, we analyzed how the percentage parameter affects the cutoff depth and how the cutoff depth affects the execution time (Figs. 3 and 4). Depending on these analyses, we set the percentage parameter to 0.5%. In addition to that, we provided a detailed analysis on the average height of the binary search trees (BSTs) constructed from hash indexes. We showed that average height of BSTs is always less than or equal the height of the original segment tree.
Furthermore, for JOA ISTF, we analyzed the number of hash indexes for varying preset values using human genome, chromosome 1 intervals of the 5 small files’ which is provided as the second benchmark dataset in Table 1. The smaller the preset value, the higher the number of different hash indexes, and vice versa (Fig. 5).
Moreover, we computed the mean and standard deviation of the number of segment tree nodes assigned to the same hash index as we vary the preset value. The smaller the preset value, the less the number of segment tree nodes assigned to a hash index (Fig. 6). This shows that there is a trade off between number of hash indexes and the mean number of segment tree nodes assigned to a hash index as we change the preset value. This is an inherent advantage of the proposed indexed segment tree forest construction algorithm and allows for generalizing to varying preset values.
Conclusion
In this paper, we presented efficient methods for parallel joint overlap analysis of n interval sets. The proposed segment tree and novel indexed segment tree forest solutions are optimized with a divide and conquer algorithm design and implemented as a tool named JOA. We showed that JOA ST and ISTF have comparable execution times. We compared JOA with other state of the art tools such as GROK, BEDTools, BEDOPS and demonstrated the efficacy of JOA with respect to the execution time. Especially, when the number of processors per node is increased, JOA had less execution time than the other tools. Furthermore, we also verified that JOA is able to identify all the overlapping intervals.
To verify the parallel processing advantage of JOA, we designed and ran an additional case study for JOA in which we found the jointly overlapping intervals of 141 interval sets of ENCODE Dnase hypersensitive sites.
Methods
Segment tree
Similar to interval trees, segment tree data structure is another wellknown space partitioning tree and a data structure for storing intervals. Its space complexity is \(\mathcal {O}(nlog{n})\) and it can be constructed in \(\mathcal {O}(nlog{n})\) time for n given intervals. Finding all intervals in the segment tree containing a query point q_{x} requires \(\mathcal {O}(log{n}+k)\) time for n given intervals and k hits [14].
Let \(I := {[x_{1} : x^{\prime }_{1}], [x_{2} : x^{\prime }_{2} ],..., [x_{n} : x^{\prime }_{n}]}\) be a set of n given intervals on the real line. Let p_{1},p_{2},...,p_{m} be the list of distinct interval endpoints (low and high endpoints), sorted from left to right. We simply partition the real line induced by these endpoints p_{i}. We name the regions in this partitioning as elementary intervals. Thus, the elementary intervals from these endpoints p_{1}, p_{2}, …, p_{m−1}, p_{m} are, from left to right, [p_{1}:p_{1}],(p_{1}:p_{2}),[p_{2}:p_{2}],...,(p_{m−1}:p_{m}),[p_{m}:p_{m}].
To this end, we build a binary search tree T whose leaves correspond to these elementary intervals induced by the endpoints of the intervals in I in an ordered way: the leftmost leaf corresponds to the leftmost elementary interval, and so on. We denote the elementary interval corresponding to a leaf μ by Int(μ).
An internal node of T represents the intervals that are the union of intervals of its two children: Int(ν) corresponding to node ν is the union of the intervals Int(ν_{leftChild}) and Int(ν_{rightChild}) in the subtree rooted at ν. The interval represented by the parent of leaf nodes, Int(ν), is the union of the elementary intervals Int(ν_{leftChild}) and Int(ν_{rightChild}). Note that elementary intervals are different from the input intervals, I.
Each internal node ν in T represents an interval Int(ν) whereas each leaf node μ represents an elementary interval, Int(μ). Each node stores a subset of input intervals, i.e., the canonical subset, I(ν), where I(ν)⊆I. This canonical subset of node ν stores the intervals [x:x^{′}]∈I such that Int(ν)⊆[x:x^{′}] and Int(parent(ν))⫅̸[x:x^{′}]. This implies that if interval [x:x^{′}] is in the canonical subset of node ν, I(ν), then Int(ν) is completely contained in [x:x^{′}] whereas Int(parent(ν)) is not contained in [x:x^{′}]. The constructed balanced binary tree T is called a segment tree. This way of construction ensures nonoverlapping, totally consecutive intervals for the nodes at any depth, from left to right. In other words, it provides a natural binning at any depth of the tree. In Fig. 7, we exemplify how we store 4 intervals in the segment tree leaves and internal nodes which are constructed from the endpoints of these 4 given intervals. Moreover, in Fig. 7, the arrows from the nodes point to their canonical subsets.
Segment tree construction complexity analysis
To construct a segment tree for n given intervals, we proceed as follows: We sort the endpoints of n intervals in O(n logn) time and define elementary intervals at each end point and between each consecutive endpoints. We then construct a binary tree on these elementary intervals, where each interval is the union of its left and right child’s intervals (or elementary intervals). This process continues until the extent of all of the intervals are represented by the root node. The construction can be done bottomup in linear time. In the last phase, n given intervals are stored in the canonical subset of the nodes, when a node’s interval, Int(ν), is fully contained in an input interval. As a result, an interval can be attached to more than one node and the number of intervals attached to nodes decreases as we go up in the tree; since, a node’s interval, Int(ν), becomes larger as we get closer to the root node.
Segment tree query
A query is processed starting at the root node. If the query point q_{x} overlaps with the node’s interval, Int(ν), the associated intervals stored at that node, I(ν), are output and the query continues on the left or right child of that node, visiting one node per level of the tree. The time complexity of a segment tree query is O(logn+k) where n is the number of intervals and k is the number of overlapping intervals in the segment tree for the query point q_{x} [14].
Indexed segment tree forest
After analyzing constructed segment trees for real data sets, we observed that nodes at the top of the segment tree (approximately top two thirds of the segment tree) do not store any intervals or store only a few intervals in their canonical subsets. In other words, input intervals are mostly stored in the bottom nodes of the segment tree.
Keeping the whole segment tree with a significant number of nodes without any or with a few intervals might be unnecessary. Furthermore, passing through all these nodes for each query in order to find overlapping intervals may increase the query time. Instead of having one tall segment tree, we can cut the segment tree at a certain depth close to the bottom of the tree and have as many short segment trees rooted at segment tree nodes present at this cutoff depth plus the segment tree nodes with no children above this cutoff depth.
Indexed segment tree forest parameters
Cutoff depth
We decide on the cutoff depth by considering the total number of intervals stored in canonical subsets of nodes at the cutoff depth and above the cutoff depth. We find a depth, d, such that the number of intervals stored from root node to the nodes at this depth is greater than 0.5% of the total number of intervals stored in the segment tree. We then choose d+1 as the cutoff depth. The closer the cutoff depth to the bottom of the tree, the more segment trees there will be in the forest. In other words, cutoff depth determines the number of short segment trees in the forest.
Moving intervals that were stored in the nodes above the cutoff depth: All the intervals attached to the nodes that are above the cutoff depth must be distributed to the nodes at the cutoff depth. Specifically, if an interval is attached to a node above the cutoff depth, then this interval must be attached to its offspring nodes at the cutoff depth. If there is no offspring node at the cutoff depth then there can be two cases. Case 1: If the node storing the interval has no offspring, then we directly add this node to our segment tree forest. Case 2: If the node has offspring/s above the cutoff depth, then we attach the interval to its lowest offspring nodes and add these lowest offspring nodes to our segment tree forest, with a higher priority given to the node closer to the cutoff depth in order to keep the order consistent between the intervals of the nodes. Note that we do this extra work for a small number of nodes.
Linking segment tree nodes at the cutoff depth to each other: To ensure fast access between consecutive segment tree nodes at the cutoff depth, we connect segment tree nodes to each other through forward and backward pointers. We call these nodes linked nodes (See Fig. 8).
As mentioned earlier, we cut the segment tree at the cutoff depth and keep the segment tree nodes at this cutoff depth and the nodes with no children above the cutoff depth in an indexed segment tree forest. Each segment tree node at this cutoff depth is in fact a root of a segment tree in the forest and we compute its hash index using a hash function for each segment tree node and we store [index,segment tree node] tuples in a map.
Cutoff depth decision, percentage parameter
We make use of percentage parameter in cutoff depth decision. This percentage parameter determines the number of intervals that needs to be moved w.r.t. the total number of intervals stored in the segment tree. As this percentage parameter increases, the number of intervals that we accept to move increases, and we cut the original segment tree closer to the leaves. We analyzed how this percentage parameter affects the cutoff depth and how the cutoff depth affects the construction and query time. For this purpose, we used an input of 2.3 million intervals and varied the percentage parameter from 0.5 to 10.0. We observed that as we increase the percentage parameter, the number of intervals that need extra movement increases and the cutoff depth decreases, as expected (See Fig. 3). In this context, cutoff depth is defined as depth from the leaf level. However, construction and query time increases slightly, most probably because of more individual nodes are added to the indexed segment tree forest data structure (See Fig. 4). Based on these results, we decide on cutoff depth by considering movement of 0.5% of the total intervals. As a result, we internally set the percentage parameter to 0.5%.
Hash function, preset value
To index the nodes at the cutoff depth we use one universal hash function as shown in Eq. 1. By using this hash function, we index these short segment trees and we access each short segment tree in \(\mathcal {O}(1)\) time instead of \(\mathcal {O}(cutoff)\) time. The preset value in the hash function determines the number of different hash indexes that one can have.
Lower preset values result in many hash indexes with fewer segment trees with the same index, therefore, fewer collisions. But, lower preset values may result in sequential search of more than one segment tree, which is definitely not preferred. Conversely, higher preset values result in fewer hash indexes with more segment tree nodes with the same index, which implies more collisions. To efficiently handle collisions, we construct a binary search tree (BST) from the segment tree nodes with the same index and store the root of this BST in the hash map. Parent nodes of these linked nodes in the BST constitute the artificial nodes as shown in Fig. 9.
This collision handling strategy implies that each segment tree may be reached in \(\mathcal {O}(height(BST))\) time instead of \(\mathcal {O}(height(OriginalSegmentTree))\) time. As long as the height of BST formed from these segment trees with the same index is less than the height of the original segment tree, search in indexed segment tree forest will be still less than search in one tall segment tree.
Theoretically, when all the nodes at the cutoff depth have the same hash index; the same segment tree will be constructed for them. Therefore, the height of the hash BST will be always less than or equal to the height of the original segment tree. To exemplify this situation, we carried out an analysis by using inputs of 2.3 million intervals. We observed that average height of the hash BST is less than the average height of the original segment tree for each chromosome. We also averaged over all chromosomes and showed that the average height of hash BSTs is 12.9 and the average height of the original segment tree is 19.4 for varying percentages from 0.5 to 10 and for the preset value of 1000000 as shown in Table 8.
Query in indexed segment tree forest
For each query interval, we compute its lowIndex and highIndex using its low and high endpoints, respectively. We start searching on a linked node pointed by the lowIndex if it exists, otherwise we find the lowerIndex (highest index lower than lowIndex) and start searching at the node shown by the lowerIndex and continue searching forward. If it is not possible, we start searching on the linked node pointed by the highIndex if it exists, if not, we compute higherIndex (lowest index higher than highIndex) and search the node pointed by higherIndex and continue searching backward. If there is no node pointed by higherIndex, it means that there is no overlapping interval with the query interval. The pseudocode of the indexed segment tree forest search algorithms are provided in the Additional file 2.
How to guarantee that at most two additional index searches are enough?
As it is shown in Fig. 10, we first compute lowIndex_{i} and highIndex_{j} using query low and high endpoints, respectively. Then we search for the segment trees pointed by one of these indexes in the order of lowIndex_{i}, lowIndex_{i−1}, highIndex_{j} or highIndex_{j+1}.
Here we present why we may need to consider only two more segment trees pointed by the indexes lowIndex_{i−1} and the highIndex_{j+1} (Fig. 10).
From the preserved order between intervals of consecutive nodes we know that
As a result of inequality Eq. 8, lowNode_{i−2}.interval and query interval can not overlap. Therefore we may need to look at only one more index preceding the lowIndex_{i} and search for the segment tree pointed by that index and forward. In the same manner, we may need to consider only one more index subsequent to the highIndex_{j} and search for the segment tree pointed by that index and backward.
Real datasets: benchmarks
We defined six benchmarks, which are listed in Table 1 for execution time comparisons of the tools.
The first two benchmarks include five files containing DNaseHypersensitivity hotspot peaks from mouse fetal tissues: fAdrenalDS12528, fAdrenalDS15123, fAdrenalDS17319, fAdrenalDS17677 and fAdrenalDS20343. They contain 193835, 188966, 137386, 132500, and 195098 intervals, respectively. The first benchmark that is called “2 small files” is a subset containing the first two of the five files, whereas the second benchmark named as “5 small files” is the set of all five files.
The last four benchmarks come from three files consisting of histone modifications from human cell lines: HepG2, Dnd41 and A459 including 2,362,386 (2.3M), 6,473,749 (6.4M), and 9,218,913 (9.2M) intervals, respectively. The last benchmark is a subset containing the last two of these files. Sources of these benchmark datasets are provided in Availability of data and material.
Abbreviations
 NGS:

Next Generation Sequencing
 JOA:

Joint Overlap Analysis
 ST:

Segment Tree
 ISTF:

Indexed Segment Tree Forest
References
 1
Speir ML, Zweig AS, Rosenbloom KR, Raney BJ, Paten B, Nejad P, Lee BT, Learned K, Karolchik D, Hinrichs AS, Heitner S, Harte RA, Haeussler M, Guruvadoo L, Fujita PA, Eisenhart C, Diekhans M, Clawson H, Casper J, Barber GP, Haussler D, Kuhn RM, Kent WJ. The UCSC Genome Browser database: 2016 update,. Nucleic Acids Res. 2016; 44(D1):717–25. https://doi.org/10.1093/nar/gkv1275.
 2
Quinlan AR, Hall IM. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics. 2010; 26(6):841–2. https://doi.org/10.1093/bioinformatics/btq033.
 3
Dale RK, Pedersen BS, Quinlan AR. Pybedtools: a flexible python library for manipulating genomic datasets and annotations. Bioinformatics. 2011; 27(24):3423–24. https://doi.org/10.1093/bioinformatics/btr539.
 4
Neph S, Kuehn MS, Reynolds AP, Haugen E, Thurman RE, Johnson AK, Rynes E, Maurano MT, Vierstra J, Thomas S, Sandstrom R, Humbert R, Stamatoyannopoulos JA. BEDOPS: highperformance genomic feature operations. Bioinformatics. 2012; 28(14):1919–20. https://doi.org/10.1093/bioinformatics/bts277.
 5
Ovaska K, Lyly L, Sahu B, Janne OA, Hautaniemi S. Genomic region operation kit for flexible processing of deep sequencing data. IEEE/ACM Trans Comput Biol Bioinforma. 2013; 10(1):200–6. https://doi.org/10.1109/TCBB.2012.170.
 6
Alekseyenko AV, Lee CJ. Nested Containment List (NCList): a new algorithm for accelerating interval query of genome alignment and interval databases. Bioinformatics (Oxford, England). 2007; 23(11):1386–93. https://doi.org/10.1093/bioinformatics/btl647.
 7
Otlu B, Firtina C, Keleş S, Tastan O. Glanet: genomic loci annotation and enrichment tool. Bioinformatics. 2017; 33(18):2818–28. https://doi.org/10.1093/bioinformatics/btx326.
 8
Layer RM, Quinlan AR. A parallel algorithm for n way interval set intersection. Proc IEEE. 2015; PP(99):1–10. https://doi.org/10.1109/JPROC.2015.2461494.
 9
Li H. Tabix: fast retrieval of sequence features from generic TABdelimited files. Bioinformatics. 2011; 27(5):718–9. https://doi.org/10.1093/bioinformatics/btq671.
 10
Li H. Bgt: efficient and flexible genotype query across many samples. Bioinformatics. 2016; 32(4):590–2. https://doi.org/10.1093/bioinformatics/btv613.
 11
Layer RM, Kindlon N, Karczewski KJ, ExAC EAC, Quinlan AR. Efficient compression and analysis of large genetic variation datasets. Nat Methods. 2015; 13:63–5. https://doi.org/10.1038/nmeth.3654.
 12
Zheng X, Gogarten SM, Lawrence M, Stilp A, Conomos MP, Weir BS, Laurie C, Levine D. Seqarray—a storageefficient highperformance data format for wgs variant calls. Bioinformatics. 2017; 33(15):2251–7. https://doi.org/10.1093/bioinformatics/btx145.
 13
Band G, Marchini J. Bgen: a binary file format for imputed genotype and haplotype data. bioRxiv. 2018. https://doi.org/10.1101/308296. https://www.biorxiv.org/content/early/2018/05/02/308296.full.pdf.
 14
de Berg M, Cheong O, van Kreveld M, Overmars M. Computational Geometry: Algorithms and Applications. 3rd ed. 2008. Berlin: Springer; 2010. https://www.amazon.com/ComputationalGeometryApplicationsMarkBerg/dp/3540779736.
Acknowledgements
Not applicable.
Funding
No funding was received for the study.
Author information
Affiliations
Contributions
BO conceived the idea of indexed segment tree forest, designed and implemented the parallel algorithms for segment tree and indexed segment tree forest, conducted execution time comparison analysis. TC supervised the project and execution time comparison analysis. BO and TC wrote the manuscript. All authors read and approved the final manuscript.
Corresponding author
Ethics declarations
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Competing interests
The authors declare that they have no competing interests.
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Additional information
Availability of data and materials
JOA’s source code is available at https://github.com/burcakotlu/JOA/. You can download JOA executable jar, joa.jar from http://burcak.ceng.metu.edu.tr/joa/. A user manual is provided at https://joa.readthedocs.org.
Benchmark Datasets: The first two benchmarks “2 small files” and “5 small files” are obtained from five files: fAdrenalDS12528, fAdrenalDS15123, fAdrenalDS17319, fAdrenalDS17677 and fAdrenalDS20343. They are all in bed format and available for download from http://burcak.ceng.metu.edu.tr/joa/. The last four benchmarks (2.3M,2.3M),(6.4M,6.4M), (9.2M,9.2M) and (6.4M,9.2M) come from three files: 2.3M (wgEncodeBroadHistoneHepg2CtcfStdAlnRep1), 6.4M (wgEncodeBroadHistoneDnd41H3k04me2AlnRep2), 9.2M (wgEncodeBroadHistoneA549H3k04me2Dex100nmAlnRep2). They are downloaded as bam format files from http://hgdownload.cse.ucsc.edu/goldenpath/hg19/encodeDCC/wgEncodeBroadHistone/ and converted into bed format files.
SemiSynthetic Dataset: All files are in bed format. Compressed simulated_data.tar.gz file can be downloaded from http://burcak.ceng.metu.edu.tr/joa/.
ENCODE Data: For additional case study of JOA, we downloaded 141 interval sets of ENCODE Dnase hypersensitive sites (29 of them are in pk format and 112 of them are in bed format) from http://hgdownload.cse.ucsc.edu/goldenpath/hg19/encodeDCC/wgEncodeUwDnase/ and http://hgdownload.cse.ucsc.edu/goldenpath/hg19/encodeDCC/wgEncodeOpenChromDnase/.
Additional files
Additional file 1
Supplementary Table S1 contains all jointly overlapping intervals for 141 interval sets of Dnase hypersensitive sites. (XLSX 75 kb)
Additional file 2
Supplementary Material includes Supplementary Tables S2 and S3 for memory usage of JOA ST and ISTF for the semisynthetic datasets. Supplementary Material also provides the pseudocode of the indexed segment tree forest search algorithms. (PDF 135 kb)
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Cite this article
Otlu, B., Can, T. JOA: Joint Overlap Analysis of multiple genomic interval sets. BMC Bioinformatics 20, 121 (2019). https://doi.org/10.1186/s1285901926984
Received:
Accepted:
Published:
Keywords
 Genome analysis
 Joint overlap analysis
 Interval overlap
 Interval intersection
 Genomic interval intersection
 Segment tree
 Indexed segment tree forest
 Space partitioning algorithms