pygenogrove 0.3.0

Python bindings for the genogrove C++ library - a specialized B+ tree for genomic intervals

pygenogrove

Python bindings for the genogrove C++ library - a specialized B+ tree data structure optimized for genomic interval storage and querying.

Installation

Building from Source

Requirements:

• C++20 compatible compiler
• CMake 3.15+
• Python 3.8+

# Clone with submodules
git clone --recursive https://github.com/genogrove/pygenogrove.git
cd pygenogrove

# Build using CMake
cmake -S . -B build -DCMAKE_BUILD_TYPE=Release
cmake --build build

# The built module will be in build/pygenogrove.so (or .pyd on Windows)

Using pip

pip install pygenogrove

Using conda/mamba

Quick Start

import pygenogrove as pg

# Create a grove with order 100 (max 99 keys per node)
grove = pg.Grove(100)

# Create intervals — coordinates are closed [start, end] (both inclusive)
interval1 = pg.Interval(100, 200)
interval2 = pg.Interval(150, 250)
interval3 = pg.Interval(300, 400)

# Insert intervals into different chromosomes
grove.insert("chr1", interval1)
grove.insert("chr1", interval2)
grove.insert("chr2", interval3)

print(f"Total intervals: {len(grove)}")  # Output: Total intervals: 3

# Query for overlapping intervals
query = pg.Interval(175, 225)
results = grove.intersect(query, "chr1")

print(f"Found {len(results)} overlapping intervals")
for key in results:
    interval = key.value
    print(f"  {interval.start}-{interval.end}")

Usage Examples

Basic Operations

import pygenogrove as pg

# Create a grove (default order is 3; minimum is 3)
grove = pg.Grove()

# Create and insert intervals
interval = pg.Interval(1000, 2000)
key = grove.insert("chr1", interval)

# Access interval properties (read-only)
print(f"Start: {interval.start}")  # Output: Start: 1000
print(f"End: {interval.end}")      # Output: End: 2000

Important — do not mutate an inserted interval. Interval.start and Interval.end are intentionally read-only, and Interval.set_range(start, end) must only be used on intervals you have NOT yet inserted (e.g. a query interval you want to reuse). Mutating a stored key silently corrupts B+ tree ordering — overlap queries will start returning wrong answers with no error.

Querying Intervals

import pygenogrove as pg

grove = pg.Grove(100)

# Insert some intervals
grove.insert("chr1", pg.Interval(100, 200))
grove.insert("chr1", pg.Interval(300, 400))
grove.insert("chr2", pg.Interval(100, 200))

# Query specific chromosome
query = pg.Interval(150, 350)
results = grove.intersect(query, "chr1")

print(f"Found {len(results)} overlaps in chr1")
for key in results:
    print(f"  Interval: {key.value}")

# Query across all chromosomes
all_results = grove.intersect(query)
print(f"Found {len(all_results)} overlaps across all chromosomes")

Overlap Detection

import pygenogrove as pg

# Static method for checking overlap
interval1 = pg.Interval(100, 200)
interval2 = pg.Interval(150, 250)
interval3 = pg.Interval(300, 400)

print(pg.Interval.overlaps(interval1, interval2))  # True (they overlap)
print(pg.Interval.overlaps(interval1, interval3))  # False (no overlap)

API Reference

Interval

Interval(start: int, end: int)

A genomic interval with closed [start, end] coordinates (0-based, both inclusive).

Attributes (read-only):

• start: Start position (inclusive)
• end: End position (inclusive)

Methods:

• set_range(start, end): Atomically set both endpoints. Only safe on intervals not yet inserted into a Grove (mutating a stored key corrupts B+ tree ordering).
• Interval.overlaps(a, b): Static method to check if two intervals overlap

Grove

Grove(order: int = 3)

A B+ tree container for genomic intervals with multi-index support.

Parameters:

• order: Maximum branching factor (max keys per node = order - 1). Minimum 3.

Methods:

• len(grove) / size() / indexed_vertex_count(): Number of indexed intervals across all indices
• get_order(): Get the order (branching factor) of the tree
• insert(index: str, interval: Interval) -> Key: Insert an interval at the specified index
• intersect(query: Interval) -> QueryResult: Find overlapping intervals across all indices
• intersect(query: Interval, index: str) -> QueryResult: Find overlapping intervals in specific index
• flanking(query: Interval, index: str) -> FlankingResult: Find the nearest non-overlapping keys on either side of the query (predecessor / successor)

FlankingResult (returned by flanking):

• predecessor: the closest key entirely before the query (a Key), or None
• successor: the closest key entirely after the query (a Key), or None

Keys overlapping the query are excluded; for nested intervals the predecessor is the one with the largest end (smallest gap). Compute the gap distance from the returned key, e.g. query.start - result.predecessor.value.end - 1 (closed coordinates). BedGrove/GffGrove expose flanking too (their results' keys carry .data).

Graph overlay (directed edges between keys):

• add_edge(source: Key, target: Key): Add a directed edge (raises ValueError if a key is None)
• remove_edge(source: Key, target: Key) -> bool: Remove an edge; True if one was removed
• has_edge(source: Key, target: Key) -> bool: Test whether an edge exists
• get_neighbors(source: Key) -> list[Key]: Keys directly reachable from source
• out_degree(source: Key) -> int: Number of outgoing edges from source
• edge_count() -> int: Total number of edges in the overlay
• vertex_count_with_edges() -> int: Number of keys with at least one outgoing edge
• add_external_key(interval: Interval) -> Key: Add a key outside the index that can still participate in the graph (not returned by intersect)

Serialization (zlib-compressed .gg binary):

• serialize(path: str): Write the grove (intervals + graph overlay) to path
• deserialize(path: str) -> Grove (static): Load a grove written by serialize

Key

Wrapper object for intervals stored in the grove. Returned by insert operations.

Attributes:

• value: The interval value of this key

QueryResult

Result object containing matching intervals from a query.

Attributes:

• query: The query interval used for the search
• keys: List of matching keys

Methods:

• __len__(): Number of results
• __iter__(): Iterate over matching keys

BedGrove (interval grove with BED data)

BedGrove is the data-carrying counterpart of Grove: each indexed interval also carries an associated BedEntry payload, so prebuilt .gg files that store BED records can be loaded, queried, and traversed from Python.

import pygenogrove as pg

g = pg.BedGrove(100)

# insert(index, interval, data) — the interval is the key, BedEntry is the payload
entry = pg.BedEntry("chr1", 1000, 2000)   # BED-native coords (0-based, half-open)
entry.name = "BRCA1"
entry.score = 900
entry.strand = "+"
key = g.insert("chr1", pg.Interval(1000, 1999), entry)

# the returned key exposes both the interval value and the BED payload
print(key.value.start, key.data.name)     # 1000 BRCA1

for hit in g.intersect(pg.Interval(1500, 1600), "chr1"):
    print(hit.data.name, hit.data.score)

# serialize/deserialize preserves the BedEntry data
g.serialize("genes.gg")
reloaded = pg.BedGrove.deserialize("genes.gg")

BedGrove exposes the same surface as Grove (multi-index insert/intersect, the graph overlay, and serialize/deserialize), with these differences:

• insert(index: str, interval: Interval, data: BedEntry) -> BedKey takes the BED payload.
• add_external_key(interval: Interval, data: BedEntry) -> BedKey takes the payload too.
• Entry-deriving inserts (no hand-conversion of coordinates):
  • insert(index, entry) -> BedKey — a 2-argument overload: pass a bare BedEntry and the Interval key is derived from its native coordinates (BED's half-open [s, e) → closed [s, e-1]; GFF's 1-based [s, e] → [s-1, e-1]). This is the foolproof way to load records from a reader.
  • insert_bulk(index, entries, presorted=False) -> list[BedKey] — same idea for a whole list of bare entries.
• Fast-path inserts (data-carrying groves only):
  • insert_sorted(index, interval, data) -> BedKey — single insert on the rightmost-append path (skips tree traversal).
  • insert_bulk(index, items, presorted=False) -> list[BedKey] — insert many explicit (Interval, BedEntry) records at once (10–20× faster for large datasets; an empty index is built bottom-up in O(n)). presorted=True assumes the records are already sorted by interval (skips the internal sort).
  • Precondition: sorted/bulk inserts require ascending intervals, and when appending to a non-empty index every new interval must be greater than all existing ones. Violating this corrupts B+ tree ordering — use plain insert if unsure. (GffGrove has all the same methods.)

BedKey is like Key but adds a data attribute:

• value: the interval (returned by copy; do not rely on mutating it)
• data: the associated BedEntry — a live, mutable reference (unlike value, the payload is not part of the B+ tree ordering, so editing it in place is safe)

BedQueryResult is the BedGrove analog of QueryResult (its keys are BedKeys).

BedEntry

A single BED record. Coordinates are BED-native: 0-based, half-open [start, end) (distinct from the closed [start, end] of Interval used as the grove key).

BedEntry(chrom: str, start: int, end: int)

Attributes (read/write):

• chrom (str), start (int), end (int)
• name: Optional[str] (BED4+)
• score: Optional[int] (BED5+)
• strand: Optional[str] — a single character ('+', '-', '.'); assigning an empty or multi-character string raises ValueError, None clears it (BED6+)
• thickness: Optional[ThickInfo] (BED7+)
• item_rgb: Optional[RgbColor] (BED9+)
• blocks: Optional[BlockInfo] (BED12)

ThickInfo(start, end), RgbColor(red, green, blue) (channels 0–255), and BlockInfo(count, sizes, starts) (with list[int] sizes/starts) are the supporting value types. List fields are returned/assigned by copy.

GffGrove (interval grove with GFF/GTF data)

GffGrove is the same data-carrying grove for GFF3/GTF records — identical surface to BedGrove, with a GffEntry payload instead of BedEntry:

import pygenogrove as pg

g = pg.GffGrove(100)

entry = pg.GffEntry("chr1", 1000, 2000, "gene")   # GFF-native coords (1-based, inclusive)
entry.source = "ensembl"
entry.strand = "+"
entry.attributes = {"gene_id": "ENSG1", "gene_name": "BRCA1"}
key = g.insert("chr1", pg.Interval(999, 1999), entry)

print(key.data.type, key.data.get_gene_id())      # gene ENSG1

for hit in g.intersect(pg.Interval(1500, 1600), "chr1"):
    print(hit.data.type, dict(hit.data.attributes))

g.serialize("genes.gg")
reloaded = pg.GffGrove.deserialize("genes.gg")

GffKey mirrors BedKey (value is a copy, data is a live mutable GffEntry reference); GffQueryResult is the GffGrove analog of QueryResult.

GffEntry

A single GFF3/GTF record. Coordinates are GFF-native: 1-based, both endpoints inclusive (distinct from Interval's 0-based closed and BedEntry's 0-based half-open).

GffEntry(seqid: str, start: int, end: int, type: str)

Attributes (read/write):

• seqid (str), source (str), type (str), start (int), end (int)
• score: Optional[float]
• strand: Optional[str] — a single character ('+', '-', '.', '?'); empty or multi-character raises ValueError, None clears it
• phase: Optional[int] (CDS phase 0/1/2)
• attributes: dict[str, str] — the column-9 key/value pairs (returned/assigned by copy)
• format: a GffFormat enum (GFF3 / GTF / UNKNOWN)

Methods: is_gtf(), is_gff3(), get_attribute(key), and the GTF helpers get_gene_id(), get_transcript_id(), get_exon_number(), get_gene_name(), get_gene_biotype() (each returns None when the attribute is absent).

BedReader / GffReader (file iterators)

BedReader and GffReader are single-pass iterators over BED and GFF3/GTF files. Iterate them to get BedEntry / GffEntry records. Plain and gzip/BGZF-compressed (.gz) files are both accepted (auto-detected).

import pygenogrove as pg

# read records one at a time
for entry in pg.BedReader("peaks.bed"):
    print(entry.chrom, entry.start, entry.end, entry.name)

# the common workflow: load a file into a grove. The 2-argument insert derives
# the grove's 0-based closed Interval key from each entry's native coordinates,
# so you don't hand-convert (BED half-open, GFF 1-based) yourself.
g = pg.BedGrove(256)
for e in pg.BedReader("peaks.bed"):
    g.insert(e.chrom, e)

gff = pg.GffGrove(256)
for e in pg.GffReader("genes.gff3"):
    gff.insert(e.seqid, e)

# bulk-load one chromosome at a time (insert_bulk is per-index):
g2 = pg.BedGrove(256)
g2.insert_bulk("chr1", [e for e in pg.BedReader("peaks.bed") if e.chrom == "chr1"])

BedReader(path: str, skip_invalid_lines: bool = False)
GffReader(path: str, skip_invalid_lines: bool = False, validate_gtf: bool = False)

• A missing/unreadable path raises on construction.
• With skip_invalid_lines=False (default) a malformed line raises RuntimeError mid-iteration; with True such lines are skipped. The first data record is validated when the reader is constructed, so a malformed first record raises immediately regardless of this flag.
• GffReader(..., validate_gtf=True) enforces the mandatory GTF2 attributes (gene_id, transcript_id).
• Both expose get_error_message() and get_current_line() for diagnostics.
• The readers are single-pass — they own an htslib file handle and cannot be restarted or iterated twice.

Coordinate systems — Interval is 0-based closed [start, end]; BedEntry is 0-based half-open [start, end); GffEntry is 1-based inclusive [start, end]. Convert deliberately when building grove keys, as shown above.

Current Status

This is an early development version. Currently exposed features:

• Basic grove and interval operations
• Insert and query functionality
• Multi-index support (per chromosome)
• Graph overlay (directed edges, external keys)
• Serialization / deserialization to compressed .gg files
• Nearest-neighbour queries: flanking() (predecessor / successor)
• Associated data: the BedEntry / GffEntry value types and the data-carrying groves BedGrove (grove<interval, bed_entry>) and GffGrove (grove<interval, gff_entry>)
• File readers: BedReader and GffReader (single-pass iterators over BED / GFF3 / GTF files, including .gz)
• Fast-path inserts on data-carrying groves: insert_sorted / insert_bulk, plus entry-deriving insert(index, entry) / insert_bulk(index, entries) overloads that compute the key from a BED/GFF record's native coordinates

Not yet exposed (tracked in #1):

• Genomic coordinates with strand information, and other key types — numeric, kmer (#7)
• BAM/SAM and FASTA readers
• Edge metadata, get_neighbors_if / link_if (require a metadata-carrying grove)

Performance Tips

1. Choose appropriate order: Higher order (e.g., 100-500) reduces tree height for large datasets
2. Separate by chromosome: Use the index parameter to maintain separate trees per chromosome
3. Query specific indices: Query specific chromosomes instead of all indices when possible

License

This project inherits the license from the genogrove C++ library and is therefore licensed under the GPLv3 license.

Related Projects

• genogrove: The underlying C++ library

Citation

If you use pygenogrove in your research, please cite the original genogrove library