Almost every genomics analysis ends up asking a question about intervals: which ChIP-seq peaks fall in promoters, how many reads overlap each gene, which variants lie inside exons, or which regions are shared between two experiments. This guide explains what genomic intervals are, the file formats that store them, the coordinate conventions that cause so many off-by-one errors, the core operations (intersect, merge, subtract, closest, coverage), the main tools compared (bedtools, GenomicRanges, PyRanges, bioframe), and how to troubleshoot the problems you are most likely to meet.
A genomic interval (also called a region, range or feature) is a stretch of a genome defined by three pieces of information: the chromosome (or contig or scaffold) name, a start position and an end position. A fourth piece, the strand (plus, minus or unknown), is often added, and extra columns carry a name, a score or other annotations.
Examples of things we describe as intervals:
A gene, an exon, a CDS or a UTR in an annotation file.
A ChIP-seq, ATAC-seq or CUT&Tag peak.
A CpG island, an enhancer, a promoter, a repeat or a transposable element.
A structural variant, a deletion or a copy-number segment.
The positions covered by a sequencing read (a read is also an interval on the reference).
A window of fixed size used for binning the genome.
Because an interval is just coordinates, you can compare very different data types with the same operations. That is what makes interval arithmetic so useful and so common in bioinformatics.
2. Why intervals matter: typical questions they answer
Annotation: which gene is closest to each peak, and how far away is it?
Overlap and intersection: how many peaks from experiment A are also found in experiment B?
Counting: how many reads overlap each gene or each bin?
Filtering: remove regions that overlap blacklisted, low-mappability or repeat regions.
Defining new regions: promoters (for example transcription start site plus or minus 2 kb), introns (gene minus exons), intergenic regions (the complement of genes).
Merging: combine nearby or overlapping peaks into consensus regions.
Extraction: obtain the DNA sequence of a list of regions, or the variants inside them.
Enrichment testing: are my regions overlapping an annotation more often than by chance?
Coordinate conversion: move regions between genome assemblies (for example hg19 to hg38).
3. File formats and their coordinate systems
Different formats describe intervals in different ways. Mixing them without checking the coordinate system is the number one cause of silent errors.
Format
Typical content
Columns (minimum)
Coordinate system
BED
Generic intervals
chrom, start, end (BED3), optionally name, score, strand (BED6) up to BED12
Not applicable; required by many bedtools commands
FASTA index (.fai)
Index of a FASTA file
name, length, offset, ...
Not applicable
4. The most important concept: 0-based versus 1-based coordinates
There are two conventions to number bases, and both are in daily use.
**0-based, half-open [start, end)**: the first base of a chromosome is position 0, the start is included, and the end is excluded. The length of the interval is simply end - start. This is used by BED, bedGraph, bigWig, BAM internals, Python and, internally, PyRanges and bedtools.
**1-based, closed [start, end]**: the first base is position 1, and both ends are included. The length is end - start + 1. This is used by GFF, GTF, VCF, SAM text, UCSC browser display, samtools region strings (chr1:1000-2000) and R's IRanges/GRanges.
The same three-base feature is therefore written differently depending on the format:
The same feature written in two conventions: BED uses 0-based, half-open coordinates, while GFF, VCF and GRanges use 1-based, closed coordinates.
Conversion rules
BED to 1-based: start + 1, end unchanged.
1-based to BED: start - 1, end unchanged.
A single base at 1-based position 100 is chr1 99 100 in BED.
A zero-length interval (insertion site) is allowed in BED (start == end) but not in GFF.
💡 Typical bug: converting a GTF to BED without subtracting 1 from the start. Every gene is then shifted by one base and the first base of each feature is lost. For overlaps of large regions this goes unnoticed, but it breaks exon-level analyses, variant annotation and motif extraction.
Other coordinate pitfalls
Chromosome names: chr1 (UCSC) versus 1 (Ensembl and NCBI). Files from different sources must use the same names, or every overlap comes back empty.
Genome assembly version: positions are meaningless unless both files use the same assembly (hg19 versus hg38, mm10 versus mm39). Convert with liftOver or CrossMap, and expect some regions not to convert.
Sorting: many tools need intervals sorted by chromosome and then by start. Lexicographic order (chr1, chr10, chr11, ..., chr2) is different from natural order (chr1, chr2, ..., chr10). Use the same order in every file.
Strand: operations may or may not respect strand. Check the flag in each tool.
Line endings and spaces: BED files must be tab-separated. Windows line endings (^M) and spaces instead of tabs break parsers.
5. Core interval operations
Most analyses are built from a small set of operations.
Operation
Question answered
bedtools command
Intersect
Which parts of A overlap B?
intersect
Subtract
What remains of A after removing B?
subtract
Merge
Combine overlapping or nearby intervals into one
merge
Closest
What is the nearest B feature to each A, and how far?
closest
Window
Which B features are within N bases of each A?
window
Slop and flank
Extend an interval, or take its upstream or downstream neighbors
slop, flank
Complement
Which regions of the genome are not covered?
complement
Coverage
How much of each A is covered by B, or how many reads fall in it?
coverage, multicov
Shuffle
Randomize positions, to build a null distribution
shuffle
Jaccard
How similar are two sets of intervals overall?
jaccard
Makewindows
Split the genome into bins of a chosen size
makewindows
Getfasta
Extract the DNA sequence of intervals
getfasta
Three basic operations on two interval sets A and B: intersect keeps shared bases, merge fuses overlapping intervals, subtract removes B from A (schematic made for Bioeduc).
6. Tools compared
Tool
Language and interface
Strengths
Limits
Best for
bedtools
Command line (C++)
Standard, fast, many operations, works with BED/GFF/VCF/BAM, integrates in pipelines
Needs sorted files for best performance, text-based
Pandas-like, fast (Rust or Cython backend), good API
API changed in versions 1.x; check your version
Python workflows
bioframe
Python (pandas)
Works directly on DataFrames, clear semantics, 3D-genome community
Less feature-rich than bedtools
Pandas users
pybedtools
Python wrapper
Calls bedtools from Python, convenient
Depends on the bedtools binary
Scripts around bedtools
tabix, bcftools, samtools
Command line
Random access to region in compressed, indexed files
Retrieval only, no complex operations
Querying large VCF, BED or BAM files by region
deepTools, IGV, UCSC tools
CLI and viewers
Signal over regions, heatmaps, visual inspection
Specific tasks
ChIP-seq/ATAC-seq analysis and QC
There is no single winner. A practical rule: use bedtools for command-line pipelines, GenomicRanges/plyranges when the rest of the analysis is in R, and PyRanges or bioframe when it is in Python.
7. bedtools in practice
Prepare the files
# Sort everything in the same order (the one in your genome file)
sort -k1,1 -k2,2n peaks.bed > peaks.sorted.bed
bedtools sort -i genes.bed -g genome.chrom.sizes > genes.sorted.bed
# Convert a GTF gene feature to BED (gawk): subtract 1 from the start!
gawk 'BEGIN{OFS="\t"} $3=="gene"{match($0,/gene_id "([^"]+)"/,a); print $1,$4-1,$5,a[1],0,$7}' genes.gtf \
| sort -k1,1 -k2,2n > genes.bed
# Chromosome sizes from a FASTA index
samtools faidx genome.fa && cut -f1,2 genome.fa.fai > genome.chrom.sizes
Intersect: all variants of the basic overlap query
# Peaks that overlap at least one gene, report each A record once
bedtools intersect -a peaks.sorted.bed -b genes.sorted.bed -u
# Report both records of every overlap pair (a peak can appear several times)
bedtools intersect -a peaks.sorted.bed -b genes.sorted.bed -wa -wb
# Add the number of overlapping bases
bedtools intersect -a peaks.sorted.bed -b genes.sorted.bed -wo
# Peaks with NO overlap (for example outside blacklist regions)
bedtools intersect -a peaks.sorted.bed -b blacklist.bed -v
# Number of overlapping genes per peak (0 if none)
bedtools intersect -a peaks.sorted.bed -b genes.sorted.bed -c
# Keep all peaks, filling with "." when there is no overlap (left outer join)
bedtools intersect -a peaks.sorted.bed -b genes.sorted.bed -loj
# Require 50 percent overlap of both features (reciprocal), same strand only
bedtools intersect -a A.bed -b B.bed -f 0.5 -r -s
Useful options: -f and -F set the minimum fraction of A or B that must overlap, -r makes it reciprocal, -s requires the same strand and -S the opposite strand, and -sorted streams pre-sorted files using little memory (add -g genome.chrom.sizes to verify chromosome order).
Annotate peaks with the nearest gene
# Nearest gene, with signed distance relative to the gene's orientation (-D ref)
bedtools closest -a peaks.sorted.bed -b genes.sorted.bed -D ref -t first > peaks_nearest_gene.tsv
Build promoters, introns and intergenic regions
# 2 kb upstream of each gene, strand-aware (promoter-like region)
bedtools flank -i genes.sorted.bed -g genome.chrom.sizes -l 2000 -r 0 -s > upstream2kb.bed
# Extend each gene by 1 kb on both sides
bedtools slop -i genes.sorted.bed -g genome.chrom.sizes -b 1000 > genes_plus1kb.bed
# Introns-like regions: gene minus exons
bedtools subtract -a genes.sorted.bed -b exons.sorted.bed > introns.bed
# Intergenic regions: complement of genes
bedtools merge -i genes.sorted.bed | bedtools complement -i - -g genome.chrom.sizes > intergenic.bed
Merge, window and bin
# Merge peaks closer than 100 bases and keep the names of merged features
bedtools merge -i peaks.sorted.bed -d 100 -c 4 -o collapse > merged_peaks.bed
# Genome bins of 10 kb
bedtools makewindows -g genome.chrom.sizes -w 10000 > bins_10kb.bed
Count reads in regions
# Reads per region for one BAM
bedtools coverage -a regions.bed -b sample.bam -counts > region_counts.tsv
# Several BAM files at once
bedtools multicov -bams s1.bam s2.bam s3.bam -bed regions.bed > counts_matrix.tsv
For gene counting in RNA-seq, specialized counters such as featureCounts handle paired reads, strandedness and multi-mapping more carefully than a generic overlap.
Sequences, statistics and conversions
# Extract sequences (uses the BED 0-based convention; -s reverse-complements minus strand)
bedtools getfasta -fi genome.fa -bed regions.bed -name -s > regions.fa
# Similarity of two peak sets
bedtools jaccard -a peaksA.sorted.bed -b peaksB.sorted.bed
# Move to another assembly
liftOver peaks_hg19.bed hg19ToHg38.over.chain.gz peaks_hg38.bed unmapped.bed
8. GenomicRanges in R
In Bioconductor, intervals are stored in GRanges objects (1-based, closed) with strand and metadata columns. rtracklayer::import() reads BED files and automatically converts to the 1-based convention, so you do not subtract or add manually.
library(GenomicRanges); library(rtracklayer)
peaks <- import("peaks.bed") # BED -> GRanges (coordinate conversion handled)
genes <- import("genes.gtf") # keeps gene_id, gene_name as metadata
genes <- genes[genes$type == "gene"]
# Overlaps
ov <- findOverlaps(peaks, genes)
hits <- peaks[queryHits(ov)]
mcols(hits)$gene <- genes$gene_id[subjectHits(ov)]
peaks_in_genes <- subsetByOverlaps(peaks, genes)
n_overlaps <- countOverlaps(peaks, genes)
# Nearest gene and distance
nn <- distanceToNearest(peaks, genes)
mcols(peaks)$nearest <- genes$gene_id[subjectHits(nn)]
mcols(peaks)$dist <- mcols(nn)$distance
# Promoters (strand aware): 2 kb upstream and 200 bp downstream of TSS
prom <- promoters(genes, upstream = 2000, downstream = 200)
# Set operations
merged <- reduce(peaks, min.gapwidth = 100)
remaining <- setdiff(genes, blacklist)
# Export back to BED
export(merged, "merged_peaks.bed")
plyranges offers tidy verbs such as find_overlaps(), join_overlap_inner(), group_by() and summarise(), which make complex overlap summaries easy to read. For peak annotation in ChIP-seq and ATAC-seq, ChIPseeker assigns peaks to promoters, UTRs, exons, introns and intergenic regions and plots the distribution.
9. PyRanges and bioframe in Python
PyRanges stores intervals in a DataFrame with Chromosome, Start, End (and optionally Strand) columns in the BED convention.
import pyranges as pr
peaks = pr.read_bed("peaks.bed")
genes = pr.read_gtf("genes.gtf")
genes = genes[genes.Feature == "gene"]
overlaps = peaks.join(genes) # one row per overlapping pair
in_genes = peaks.overlap(genes) # peaks overlapping any gene
n_overlap = peaks.count_overlaps(genes) # number of overlaps per peak
nearest = peaks.nearest(genes) # nearest gene and distance
merged = peaks.merge(slack=100) # merge nearby peaks
remaining = genes.subtract(blacklist) # remove blacklisted bases
💡 PyRanges 1.x renamed several methods compared with 0.x versions (for example the join functions). Check the documentation of the version you installed.
bioframe works directly on pandas DataFrames with chrom, start, end columns, which is convenient if your data are already tables:
import bioframe as bf
df_peaks = bf.read_table("peaks.bed", schema="bed3")
df_genes = bf.read_table("genes.bed", schema="bed6")
inter = bf.overlap(df_peaks, df_genes, how="inner", suffixes=("_peak", "_gene"))
merged = bf.merge(df_peaks, min_dist=100)
close = bf.closest(df_peaks, df_genes)
10. Region-based access to variants, reads and signal
Large files should not be loaded fully into memory when you only need a few regions.
# VCF: variants inside regions (indexed VCF: -R; streaming without index: -T)
bcftools view -R regions.bed.gz variants.vcf.gz -Oz -o variants_in_regions.vcf.gz
# BAM: reads overlapping regions
samtools view -L regions.bed -b sample.bam > sample.in_regions.bam
# Compress and index a BED file for random access
bgzip regions.sorted.bed && tabix -p bed regions.sorted.bed.gz
tabix regions.sorted.bed.gz chr2:1,000,000-1,050,000 # query uses 1-based coordinates
11. Testing whether an overlap is significant
Finding that 40 percent of your peaks overlap promoters does not mean much unless you know what to expect by chance. Build a null model that preserves what matters, such as the interval lengths and the chromosome distribution, then compare.
bedtools shuffle randomizes intervals while respecting chromosome sizes. Use -excl to avoid gaps and blacklisted regions, and -chrom to keep intervals on their original chromosomes.
bedtools fisher computes a Fisher exact test on the overlap of two sets.
regioneR (R) and GAT (Python/command line) run permutation tests with flexible randomization strategies.
Choose the background carefully: a random draw from the whole genome is a poor null when your regions can only occur in accessible, mappable or gene-rich areas.
Sort once, then stream: sorted inputs let bedtools use -sorted, which reads files in a single pass with very little memory instead of building an index in RAM.
Compress and index with bgzip and tabix for repeated region queries.
Filter early: restrict to the chromosomes or regions you need before expensive joins.
Avoid many-to-many blowups: if a peak overlaps hundreds of features, -wa -wb produces hundreds of rows. Use -u, -c or closest -t first when one row per peak is enough.
In R and Python, keep intervals in native objects (GRanges, PyRanges), not plain strings, and avoid looping over rows: use vectorized overlaps instead.
For very large single-cell fragment files (scATAC-seq), use tabix-indexed fragments and tools designed for them (for example Signac, ArchR, SnapATAC2).
13. Troubleshooting: symptoms, causes and fixes
Symptom or error
Likely cause
What to do
Intersect returns nothing
Different chromosome names (chr1 versus 1) or different assemblies
List the chromosome names of each file with cut -f1 and sort -u, compare them, rename with sed 's/^chr//', and verify the assembly
ERROR: chromosome sort ordering or unsorted warning
Files sorted differently, or not sorted
Sort both with the same command or bedtools sort -g genome.chrom.sizes
Error: malformed BED entry
Spaces instead of tabs, header lines, missing columns, Windows line endings
inspect hidden characters with cat -A, convert with dos2unix, remove track or browser lines
Everything is shifted by one base
Mixed 0-based and 1-based coordinates
Check the source format, convert GTF/VCF starts with start - 1 for BED
Gene counts differ between tools
Open versus closed intervals, strandedness, how multi-mappers or overlapping genes are treated
Use a dedicated counter for RNA-seq and document the options
Duplicate peaks after intersect
Peak overlaps several features with -wa -wb
Use -u, -c, or collapse with groupBy; use closest -t first for one hit
-s or -S gives no results
Strand column missing or .
Provide a BED6 file with strand, or drop the strand option
liftOver loses many regions
Regions fall in non-syntenic or rearranged sequence
Inspect the unmapped.bed, adjust -minMatch, or remap the original data
tabix complains about compression or index
File not compressed with bgzip or not sorted
Run bgzip (not gzip), sort first, then tabix -p bed
getfasta returns the wrong bases
Wrong coordinate convention or wrong FASTA assembly
Verify with a feature of known sequence and ensure the FASTA index matches the assembly
Out of memory in R or Python
Millions of intervals with large metadata, or accidental many-to-many join
Subset columns, use sorted streaming tools, or process by chromosome
Zero-length intervals dropped
Insertions or BED records with start == end
Pad by one base or use tools that support empty intervals
GRanges end smaller than start
Importing BED without rtracklayer conversion
Use import() instead of reading the table manually
A fast sanity check for any new overlap pipeline: take one feature of which you know the position, run the pipeline on it, and confirm with a genome browser (IGV or UCSC) that the output coordinates match what you see.
14. Best-practice checklist
Know the coordinate system of every input file before combining them.
Use the same assembly version and chromosome naming in all files.
Sort and compress large files, and keep a chrom.sizes file for the genome.
Decide explicitly whether operations should be strand-aware and whether partial overlaps count.
Choose a meaningful null model before interpreting overlap enrichment.
Save the exact commands and tool versions that produced each region file.
Inspect a few results manually in a genome browser.
15. Frequently asked questions
Is BED 0-based or 1-based? BED start is 0-based and the end is exclusive. The first base of a chromosome is 0 1.
How do I convert between BED and GFF coordinates? BED start equals GFF start minus 1, and the end is the same. Reverse the operation going back.
Should I use bedtools or an R/Python package? Use bedtools for pipelines and simple operations on files, and R or Python packages when overlaps are one step of a larger analysis in those languages.
How do I handle overlapping genes? Decide whether to keep all matches, assign by priority (for example protein-coding first), or use strand information to separate them. Document the rule.
What is a blacklist? A set of genomic regions with artefactual signal in sequencing assays (for example the ENCODE blacklist). Peaks overlapping them are normally removed.
16. References
Quinlan AR, Hall IM (2010). BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics 26:841-842.
Lawrence M et al. (2013). Software for computing and annotating genomic ranges. PLoS Computational Biology 9:e1003118.
Stovner EB, Saetrom P (2020). PyRanges: efficient comparison of genomic intervals in Python. Bioinformatics 36:918-919.
Open2C et al. (2024). Bioframe: operations on genomic intervals in pandas dataframes. Bioinformatics 40:btae088.
Lee S, Cook D, Lawrence M (2019). plyranges: a grammar of genomic data transformation. Genome Biology 20:4.
Yu G, Wang LG, He QY (2015). ChIPseeker: an R/Bioconductor package for ChIP peak annotation, comparison and visualization. Bioinformatics 31:2382-2383.
Amemiya HM, Kundaje A, Boyle AP (2019). The ENCODE blacklist: identification of problematic regions of the genome. Scientific Reports 9:9354.
Gel B et al. (2016). regioneR: an R/Bioconductor package for the association analysis of genomic regions based on permutation tests. Bioinformatics 32:289-291.
💬 Commentaires