Pybedtools

pybedtools

Use this skill for genome interval set operations from Python. It wraps BEDTools and is appropriate for ATAC-seq peaks, ChIP-seq peaks, blacklist filtering, enhancer/gene overlap, GWAS locus annotation, VCF region filtering, BAM coverage, and closest-feature queries.

Setup

pybedtools requires both the Python package and the bedtools binary.

conda install -c bioconda bedtools pybedtools

Use pip install pybedtools only when bedtools is already available on PATH.

Common Patterns

import pybedtools

peaks = pybedtools.BedTool("peaks.bed").sort()
genes = pybedtools.BedTool("genes.gtf").sort()

overlap = peaks.intersect(genes, wa=True, wb=True)
merged = peaks.merge()
nearest = peaks.closest(genes, d=True)
coverage = genes.coverage(peaks)

DataFrame integration:

bt = pybedtools.BedTool.from_dataframe(df[["chrom", "start", "end", "name"]])
out = bt.intersect("targets.bed", wa=True, wb=True).to_dataframe()

Decision Rules

• Use pysam for low-level BAM/VCF reading and writing.
• Use pybedtools for set algebra over genomic intervals.
• Use duckdb or polars for large tabular aggregation after intervals are computed.
• Sort inputs before operations where BEDTools expects sorted data.

Pitfalls

• BED is 0-based half-open; GFF/GTF/VCF are usually 1-based. Convert deliberately.
• Chromosome naming must match (chr1 vs 1).
• Temporary files can accumulate in long sessions; call pybedtools.cleanup() when appropriate.
• For reproducible randomization, set genome files and seeds explicitly in shuffle workflows.