FeatureDB

The query database. All methods below are rendered from the live docstrings.

Close it, or use with

DuckDB takes an exclusive lock on the database file for the life of a writable handle. Use a with block, or call close():

with FeatureDB("gencode.duckdb") as db:
    ...

db = FeatureDB("gencode.duckdb", read_only=True)   # shareable by N processes

read_only=True is what lets several worker processes read one annotation database at once. See Connections & concurrency.

Use the _batched methods for bulk work

children(), parents() and region() return Feature objects one at a time — right for exploring, wrong for feeding a model. children_batched(), parents_batched() and region_batched() answer the same question for thousands of anchors in a single query, returning Arrow / pandas / polars without constructing any Feature objects.

gffbase.interface.FeatureDB

FeatureDB(dbfn, default_encoding: str = 'utf-8', keep_order: bool = False, pragmas: dict | None = None, sort_attribute_values: bool = False, text_factory=str, upgrade: str = 'auto', read_only: bool = False, _own_conn: bool | None = None)

Drop-in successor to gffutils.FeatureDB.

Source code in python/gffbase/interface.py

def __init__(
    self,
    dbfn,
    default_encoding: str = "utf-8",
    keep_order: bool = False,
    pragmas: dict | None = None,
    sort_attribute_values: bool = False,
    text_factory=str,
    upgrade: str = "auto",
    read_only: bool = False,
    _own_conn: bool | None = None,
):
    # These four come FIRST, before anything that can raise. `__del__`
    # runs on a half-constructed object too, and reading an attribute that
    # was never assigned would raise a second exception during garbage
    # collection, masking the first.
    self._closed = False
    self._owns_conn = False
    self._seg_cursor: duckdb.DuckDBPyConnection | None = None
    self._read_only = bool(read_only)

    if upgrade not in ("auto", "never", "error"):
        raise ValueError(f"upgrade must be 'auto', 'never' or 'error'; got {upgrade!r}")
    # A read-only handle cannot migrate: the upgrade is DDL. Coerce rather
    # than relying on DuckDB to refuse the write -- `_try_migrate` catches
    # `duckdb.Error` and falls through to v1 compatibility mode, so the
    # attempt would only produce a confusing log line on the way to the
    # same place. `upgrade="error"` is left alone: the caller asked to be
    # told about a v1 database, and they still are.
    if self._read_only and upgrade == "auto":
        upgrade = "never"
    self._upgrade = upgrade
    self.default_encoding = default_encoding
    #: Specification violations tolerated while building this database.
    #: Populated under `mode="compat"`; empty for a reopened database,
    #: which has no record of how it was built.
    self.warnings: list[dict] = []
    self.keep_order = keep_order
    self.sort_attribute_values = sort_attribute_values
    self.text_factory = text_factory
    self._analyzed_flag = False

    # Resolve dbfn → connection.
    #
    # Ownership decides what `close()` is allowed to close. gffbase must
    # never close a connection the caller opened and still holds a
    # reference to -- but it MUST close one it opened itself, or
    # `with create_db(...) as db:` would leak the handle it was written to
    # release. `_own_conn` lets `create_db` transfer ownership explicitly
    # rather than making the rule depend on which branch ran.
    if isinstance(dbfn, duckdb.DuckDBPyConnection):
        self.conn = dbfn
        self.dbfn = ":existing-connection:"
        self._owns_conn = bool(_own_conn)
    elif (
        isinstance(dbfn, tuple)
        and len(dbfn) == 2
        and isinstance(dbfn[0], duckdb.DuckDBPyConnection)
    ):
        # (con, IngestStats) — used internally by `create_db`.
        self.conn = dbfn[0]
        self.dbfn = ":existing-connection:"
        self.warnings = list(getattr(dbfn[1], "warnings", []) or [])
        self._owns_conn = True if _own_conn is None else bool(_own_conn)
    elif isinstance(dbfn, (str, os.PathLike)):
        # `os.PathLike`, not just `str`: the error below says "dbfn must be
        # a path" and this branch used to reject an actual `pathlib.Path`,
        # which is what a caller reaches for first. `create_db` already
        # accepted one, so the two entry points disagreed about their own
        # documented type.
        self.dbfn = os.fspath(dbfn)
        # Reject an embedded NUL before ANY filesystem call. DuckDB is C++
        # and takes the string as a C string, so it stops at the NUL:
        # `FeatureDB("a\0b.duckdb")` created a file called `a` -- a
        # different path than the one asked for. The stray-file cleanup
        # below then called `os.unlink` with the original path, which
        # raises `ValueError` rather than `OSError`, so it escaped the
        # `except` and left the truncated file behind under a confusing
        # error. Anything building a path from untrusted input could write
        # to a location the caller never named.
        if "\x00" in self.dbfn:
            raise ValueError(
                f"database path contains an embedded NUL byte: {self.dbfn!r}. "
                "Paths are passed to DuckDB as C strings, which would "
                "silently truncate at the NUL and open a different file."
            )
        # Whether the file was there BEFORE we connected. DuckDB creates a
        # database on connect, so after the call it always exists and the
        # question can no longer be asked -- which is how a typo'd filename
        # used to end up as an empty database plus a confusing
        # `CatalogException` from the first query against it.
        existed = os.path.exists(self.dbfn)
        if self._read_only and not existed:
            raise FileNotFoundError(
                f"no such database: {self.dbfn}. read_only=True cannot create one; "
                "build it with create_db() first."
            )
        self.conn = duckdb.connect(self.dbfn, read_only=self._read_only)
        self._owns_conn = True if _own_conn is None else bool(_own_conn)
        if not existed:
            # Undo the file DuckDB just made, so a mistyped path does not
            # litter the working directory with empty databases.
            self.conn.close()
            self._closed = True
            for path in (self.dbfn, self.dbfn + ".wal"):
                try:
                    os.unlink(path)
                except (OSError, ValueError):
                    # ValueError too: a path Python refuses to even look at
                    # must not turn best-effort cleanup into the exception
                    # the caller sees instead of the real diagnosis.
                    pass
            raise FileNotFoundError(
                f"no such database: {self.dbfn}. To create one, use "
                f"create_db(source, {self.dbfn!r}); FeatureDB() only opens "
                "databases that already exist."
            )
    else:
        raise TypeError(
            f"dbfn must be a path, DuckDB connection, or (con, stats) tuple; got {type(dbfn)!r}"
        )

    if pragmas:
        self.set_pragmas(pragmas)

    # Recover provenance from `meta`.
    meta = self._read_meta()
    self._apply_schema_version(meta)
    self.dialect = self._parse_dialect(meta.get("dialect"))
    self.fmt = meta.get("fmt", "gff3")
    # How this database was built. Absent from anything written before the
    # key existed, and absent from a v1 database; both were compat, which
    # is what the default says.
    self.mode = meta.get("mode", MODE_COMPAT)
    self.validation = meta.get("validation", VALIDATION_GFFUTILS)
    self.on_error = meta.get("on_error", ON_ERROR_RAISE)
    self._rtree_built = meta.get("rtree_built", "false").lower() == "true"
    self._max_depth = int(meta.get("max_depth", "8"))
    # Defensive: confirm the R-tree index actually exists in this DB.
    if self._rtree_built:
        self._rtree_built = self._has_rtree_index()
    # If the R-tree was built at ingest time, the spatial extension's
    # functions (ST_Intersects, ST_MakeEnvelope, …) must be loaded into
    # the current connection — they are NOT auto-loaded by DuckDB just
    # because the index exists. If the load fails (offline HPC node, etc),
    # gracefully fall back to the multi-column B-tree path.
    if self._rtree_built:
        try:
            self.conn.execute("LOAD spatial")
        except duckdb.Error:
            try:
                self.conn.execute("INSTALL spatial")
                self.conn.execute("LOAD spatial")
            except duckdb.Error:
                self._rtree_built = False

    # Closure-cache vs dynamic-CTE dispatcher. Read the corpus's
    # true hierarchy depth once at open; fall back to a live MAX(depth)
    # query if older DBs don't carry the meta row.
    cmd_meta = meta.get("closure_max_depth")
    if cmd_meta is not None:
        self._closure_max_depth = int(cmd_meta)
    else:
        try:
            row = self.conn.execute("SELECT MAX(depth) FROM closure").fetchone()
            self._closure_max_depth = int(row[0]) if row and row[0] is not None else 0
        except duckdb.Error:
            self._closure_max_depth = 0

    # Load the seqid → y-band map so `_region_sql_rtree` can
    # produce a tightly-bounded ST_MakeEnvelope at query time. Empty when
    # no R-tree was built (we fall back to the B-tree path anyway).
    self._seqid_y_map: dict = {}
    if self._rtree_built:
        try:
            rows = self.conn.execute("SELECT seqid, seqid_y FROM seqid_map").fetchall()
            self._seqid_y_map = {s: int(y) for s, y in rows}
        except duckdb.Error:
            # Older DBs without the map — fall back to B-tree to be safe.
            self._rtree_built = False

    # Directives
    self.directives = [
        row[0]
        for row in self.conn.execute("SELECT directive FROM directives ORDER BY seq").fetchall()
    ]

read_only property

read_only: bool

True if this handle refuses writes.

closed property

closed: bool

True once close() has run.

derived_source property

derived_source: str

source for features this database derives rather than reads.

gffutils_derived under compat so ported scripts that filter on it keep working; gffbase_derived under strict, which reports honest provenance.

schema

schema() -> str

The database schema as SQL text.

A METHOD, not a property: gffutils documents db.schema() and callers write it that way. Exposing it as a property meant the documented call raised TypeError: 'str' object is not callable.

Source code in python/gffbase/interface.py

def schema(self) -> str:
    """The database schema as SQL text.

    A METHOD, not a property: gffutils documents `db.schema()` and callers
    write it that way. Exposing it as a property meant the documented call
    raised `TypeError: 'str' object is not callable`.
    """
    self._require_open("schema")
    rows = self.conn.execute("""
        SELECT sql FROM duckdb_tables() WHERE database_name = current_database()
        UNION ALL
        SELECT sql FROM duckdb_views()  WHERE database_name = current_database()
    """).fetchall()
    return "\n".join(r[0] for r in rows if r[0])

close

close() -> None

Release the DuckDB connection. Idempotent.

The lazily-created segment cursor is closed whether or not the connection is owned -- gffbase created it via conn.cursor(), so gffbase closes it. The connection itself is closed only when this handle opened it: a caller who passed their own connection in still holds it afterwards.

Source code in python/gffbase/interface.py

def close(self) -> None:
    """Release the DuckDB connection. Idempotent.

    The lazily-created segment cursor is closed whether or not the
    connection is owned -- gffbase created it via `conn.cursor()`, so
    gffbase closes it. The connection itself is closed only when this
    handle opened it: a caller who passed their own connection in still
    holds it afterwards.
    """
    if self._closed:
        return
    self._closed = True
    if self._seg_cursor is not None:
        try:
            self._seg_cursor.close()
        except duckdb.Error:  # pragma: no cover - already dead
            pass
        self._seg_cursor = None
    if self._owns_conn:
        try:
            self.conn.close()
        except duckdb.Error:  # pragma: no cover - already dead
            pass

validate

validate(level: str = 'fast', **kwargs)

Check this database's structural invariants.

See :func:gffbase.validate.validate_db. Run automatically at the end of a strict-mode ingest; worth running by hand after update(), delete() or coalesce_multipart(), which are the operations that can leave the two halves of a discontinuous feature disagreeing.

Source code in python/gffbase/interface.py

def validate(self, level: str = "fast", **kwargs):
    """Check this database's structural invariants.

    See :func:`gffbase.validate.validate_db`. Run automatically at the end
    of a strict-mode ingest; worth running by hand after `update()`,
    `delete()` or `coalesce_multipart()`, which are the operations that can
    leave the two halves of a discontinuous feature disagreeing.
    """
    from gffbase.validate import validate_db

    return validate_db(self, level=level, **kwargs)

count_features_of_type

count_features_of_type(featuretype: str | None = None) -> int

Count features, optionally of one type.

Parameters:

• featuretype (str | None, default: None ) –

  Restrict the count to this GFF column-3 value ("gene", "exon", …). None counts every feature.

Returns:

• int –

  The number of matching features. A discontinuous feature counts

• int –

  once, however many input lines it was built from.

Example

db.count_features_of_type()          # 6
db.count_features_of_type("exon")    # 3

Source code in python/gffbase/interface.py

def count_features_of_type(self, featuretype: str | None = None) -> int:
    """Count features, optionally of one type.

    Args:
        featuretype: Restrict the count to this GFF column-3 value
            (`"gene"`, `"exon"`, …). `None` counts every feature.

    Returns:
        The number of matching features. A discontinuous feature counts
        once, however many input lines it was built from.

    Example:
        ```python
        db.count_features_of_type()          # 6
        db.count_features_of_type("exon")    # 3
        ```
    """
    self._require_open("count_features_of_type")
    if featuretype is None:
        return scalar(self.conn, "SELECT COUNT(*) FROM features")
    return scalar(
        self.conn, "SELECT COUNT(*) FROM features WHERE featuretype = ?", [featuretype]
    )

featuretypes

featuretypes() -> Iterator[str]

Yield every distinct featuretype in the database, alphabetically.

Yields:

• str –

  Each distinct GFF column-3 value, once.

Example

sorted(db.featuretypes())   # ['CDS', 'exon', 'gene', 'mRNA']

Source code in python/gffbase/interface.py

def featuretypes(self) -> Iterator[str]:
    """Yield every distinct featuretype in the database, alphabetically.

    Yields:
        Each distinct GFF column-3 value, once.

    Example:
        ```python
        sorted(db.featuretypes())   # ['CDS', 'exon', 'gene', 'mRNA']
        ```
    """
    self._require_open("featuretypes")
    for (ft,) in self.conn.execute(
        "SELECT DISTINCT featuretype FROM features ORDER BY featuretype"
    ).fetchall():
        yield ft

seqids

seqids() -> Iterator[str]

Yield every distinct sequence id in the database, alphabetically.

Useful for checking naming convention before a region() query -- chr1, 1 and NC_000001.11 are three different sequences as far as the database is concerned, and GENCODE, Ensembl and RefSeq each pick a different one.

Yields:

• str –

  Each distinct GFF column-1 value, once.

Source code in python/gffbase/interface.py

def seqids(self) -> Iterator[str]:
    """Yield every distinct sequence id in the database, alphabetically.

    Useful for checking naming convention before a `region()` query --
    `chr1`, `1` and `NC_000001.11` are three different sequences as far as
    the database is concerned, and GENCODE, Ensembl and RefSeq each pick a
    different one.

    Yields:
        Each distinct GFF column-1 value, once.
    """
    self._require_open("seqids")
    for (s,) in self.conn.execute(
        "SELECT DISTINCT seqid FROM features ORDER BY seqid"
    ).fetchall():
        yield s

all_features

all_features(limit: RegionLike | None = None, strand: str | None = None, featuretype: str | list[str] | None = None, order_by: str | None = None, reverse: bool = False, completely_within: bool = False) -> Iterator[Feature]

Iterate over every feature in the database.

Parameters:

• limit (RegionLike | None, default: None ) –

  Restrict to a genomic region, as "seqid:start-end" or a (seqid, start, end) tuple. None scans everything.

• strand (str | None, default: None ) –

  Restrict to "+", "-" or ".".

• featuretype (str | list[str] | None, default: None ) –

  One featuretype, or a list of them.

• order_by (str | None, default: None ) –

  Column to sort by. One of id, seqid, source, featuretype, start, end, score, strand, frame, attributes, extra, file_order, length. Anything else raises ValueError -- this is a whitelist, not a SQL fragment.

• reverse (bool, default: False ) –

  Sort descending.

• completely_within (bool, default: False ) –

  With limit, return only features contained entirely inside the region rather than merely overlapping it.

Yields:

• Feature –

  Feature objects in file_order unless order_by says otherwise.

Raises:

• ValueError –

  order_by names something outside the whitelist.

Example

for feature in db.all_features(featuretype="exon", order_by="start"):
    print(feature.id, feature.start)

Source code in python/gffbase/interface.py

def all_features(
    self,
    limit: RegionLike | None = None,
    strand: str | None = None,
    featuretype: str | list[str] | None = None,
    order_by: str | None = None,
    reverse: bool = False,
    completely_within: bool = False,
) -> Iterator[Feature]:
    """Iterate over every feature in the database.

    Args:
        limit: Restrict to a genomic region, as `"seqid:start-end"` or a
            `(seqid, start, end)` tuple. `None` scans everything.
        strand: Restrict to `"+"`, `"-"` or `"."`.
        featuretype: One featuretype, or a list of them.
        order_by: Column to sort by. One of `id`, `seqid`, `source`,
            `featuretype`, `start`, `end`, `score`, `strand`, `frame`,
            `attributes`, `extra`, `file_order`, `length`. Anything else
            raises `ValueError` -- this is a whitelist, not a SQL fragment.
        reverse: Sort descending.
        completely_within: With `limit`, return only features contained
            entirely inside the region rather than merely overlapping it.

    Yields:
        `Feature` objects in `file_order` unless `order_by` says otherwise.

    Raises:
        ValueError: `order_by` names something outside the whitelist.

    Example:
        ```python
        for feature in db.all_features(featuretype="exon", order_by="start"):
            print(feature.id, feature.start)
        ```
    """
    sql, params = self._build_scan_sql(
        base_where=[],
        base_params=[],
        limit=limit,
        strand=strand,
        featuretype=featuretype,
        order_by=order_by,
        reverse=reverse,
        completely_within=completely_within,
    )
    yield from self._yield_features(sql, params)

method

method(*args, **kwargs) -> Iterator[Feature]

Alias for all_features, kept because gffutils defines one.

Source code in python/gffbase/interface.py

def method(self, *args, **kwargs) -> Iterator[Feature]:
    """Alias for `all_features`, kept because gffutils defines one."""
    return self.all_features(*args, **kwargs)

features_of_type

features_of_type(featuretype: str | list[str], limit: RegionLike | None = None, strand: str | None = None, order_by: str | None = None, reverse: bool = False, completely_within: bool = False) -> Iterator[Feature]

Iterate over every feature of one type (or several).

Equivalent to all_features(featuretype=...); both exist because gffutils has both.

Parameters:

• featuretype (str | list[str]) –

  One featuretype ("exon"), or a list of them.

• limit (RegionLike | None, default: None ) –

  Restrict to a genomic region -- see all_features.

• strand (str | None, default: None ) –

  Restrict to "+", "-" or ".".

• order_by (str | None, default: None ) –

  Column to sort by -- see all_features for the permitted names.

• reverse (bool, default: False ) –

  Sort descending.

• completely_within (bool, default: False ) –

  With limit, require full containment.

Yields:

• Feature –

  Matching Feature objects.

Example

genes = list(db.features_of_type("gene"))
both = list(db.features_of_type(["exon", "CDS"]))

Source code in python/gffbase/interface.py

def features_of_type(
    self,
    featuretype: str | list[str],
    limit: RegionLike | None = None,
    strand: str | None = None,
    order_by: str | None = None,
    reverse: bool = False,
    completely_within: bool = False,
) -> Iterator[Feature]:
    """Iterate over every feature of one type (or several).

    Equivalent to `all_features(featuretype=...)`; both exist because
    gffutils has both.

    Args:
        featuretype: One featuretype (`"exon"`), or a list of them.
        limit: Restrict to a genomic region -- see `all_features`.
        strand: Restrict to `"+"`, `"-"` or `"."`.
        order_by: Column to sort by -- see `all_features` for the
            permitted names.
        reverse: Sort descending.
        completely_within: With `limit`, require full containment.

    Yields:
        Matching `Feature` objects.

    Example:
        ```python
        genes = list(db.features_of_type("gene"))
        both = list(db.features_of_type(["exon", "CDS"]))
        ```
    """
    yield from self.all_features(
        limit=limit,
        strand=strand,
        featuretype=featuretype,
        order_by=order_by,
        reverse=reverse,
        completely_within=completely_within,
    )

region

region(region=None, seqid: str | None = None, start: int | None = None, end: int | None = None, strand: str | None = None, featuretype: str | list[str] | None = None, completely_within: bool = False) -> Iterator[Feature]

Features overlapping a genomic interval.

The interval can be given as a string, a tuple, or the three keyword arguments::

db.region("chr1:1000-2000")
db.region(("chr1", 1000, 2000))
db.region(seqid="chr1", start=1000, end=2000)

Parameters

region : str or tuple "chrom:start-stop", or (chrom, start, stop). seqid, start, end : optional The same interval, spelled out. Mutually exclusive with region. strand : {"+", "-", "."}, optional Restrict to one orientation. featuretype : str or list of str, optional Restrict to one or several types. completely_within : bool False (default) returns anything that OVERLAPS the interval; True returns only features contained entirely within it.

Which index answers the query is chosen here, not by the caller: an R-tree when one was built and the interval is fully specified, the multi-column B-tree otherwise. The two are semantically identical -- tests/test_spatial_parity.py asserts they return the same features -- so this is a planner decision, not a behavioural one.

A feature with no coordinates (a GFF row carrying . in columns 4 and 5) is not in coordinate space and is never returned.

Source code in python/gffbase/interface.py

def region(
    self,
    region=None,
    seqid: str | None = None,
    start: int | None = None,
    end: int | None = None,
    strand: str | None = None,
    featuretype: str | list[str] | None = None,
    completely_within: bool = False,
) -> Iterator[Feature]:
    """Features overlapping a genomic interval.

    The interval can be given as a string, a tuple, or the three keyword
    arguments::

        db.region("chr1:1000-2000")
        db.region(("chr1", 1000, 2000))
        db.region(seqid="chr1", start=1000, end=2000)

    Parameters
    ----------
    region : str or tuple
        `"chrom:start-stop"`, or `(chrom, start, stop)`.
    seqid, start, end : optional
        The same interval, spelled out. Mutually exclusive with `region`.
    strand : {"+", "-", "."}, optional
        Restrict to one orientation.
    featuretype : str or list of str, optional
        Restrict to one or several types.
    completely_within : bool
        False (default) returns anything that OVERLAPS the interval; True
        returns only features contained entirely within it.

    Which index answers the query is chosen here, not by the caller: an
    R-tree when one was built and the interval is fully specified, the
    multi-column B-tree otherwise. The two are semantically identical --
    `tests/test_spatial_parity.py` asserts they return the same features --
    so this is a planner decision, not a behavioural one.

    A feature with no coordinates (a GFF row carrying `.` in columns 4 and
    5) is not in coordinate space and is never returned.
    """
    rseqid, rstart, rend = self._normalize_region_args(region, seqid, start, end)
    # Decide path.
    use_rtree = (
        self._rtree_built and rseqid is not None and rstart is not None and rend is not None
    )
    if use_rtree:
        sql, params = self._region_sql_rtree(
            rseqid, rstart, rend, strand, featuretype, completely_within
        )
    else:
        sql, params = self._region_sql_btree(
            rseqid, rstart, rend, strand, featuretype, completely_within
        )
    yield from self._yield_features(sql, params)

to_table

to_table(featuretype: str | list[str] | None = None, *, format: str = 'arrow', limit: RegionLike | None = None, strand: str | None = None, order_by: str | None = None, reverse: bool = False, completely_within: bool = False) -> Any

Return the whole database (or a slice of it) as one table.

The columnar counterpart to all_features(): same filters, but the result is a pyarrow.Table / pandas.DataFrame / polars.DataFrame instead of a stream of Feature objects, and no Feature is constructed at any layer.

For "give me every exon as a dataframe", the row-by-row API pays a Python object per row -- millions of them on a whole-genome corpus -- and that allocation dominates everything else. This is one query and one hand-off.

Parameters:

• featuretype (str | list[str] | None, default: None ) –

  One featuretype, or a list of them. None is everything.

• format (str, default: 'arrow' ) –

  "arrow" (default, zero-copy), "df" for pandas, or "polars".

• limit (RegionLike | None, default: None ) –

  Restrict to a genomic region, as "seqid:start-end" or a (seqid, start, end) tuple.

• strand (str | None, default: None ) –

  Restrict to "+", "-" or ".".

• order_by (str | None, default: None ) –

  Column to sort by -- the same whitelist all_features accepts.

• reverse (bool, default: False ) –

  Sort descending.

• completely_within (bool, default: False ) –

  With limit, require full containment.

Returns:

• Any –

  A table in the shape named by format, with one row per feature and the columns of the features table (id, seqid, source, featuretype, start, end, score, strand, frame, file_order).

Raises:

• ValueError –

  order_by names something outside the whitelist, or format is not one of the three.

• ImportError –

  format="df"/"polars" without that package.

Example

exons = db.to_table("exon", format="arrow")
df = db.to_table(["exon", "CDS"], format="df", limit="chr1:1-10000")

Note

Attributes are not included: they are a long-form table, so flattening them would either invent a column per key or collapse multi-valued keys. Query the attributes table with execute() when you need them.

Source code in python/gffbase/interface.py

def to_table(
    self,
    featuretype: str | list[str] | None = None,
    *,
    format: str = "arrow",
    limit: RegionLike | None = None,
    strand: str | None = None,
    order_by: str | None = None,
    reverse: bool = False,
    completely_within: bool = False,
) -> Any:
    """Return the whole database (or a slice of it) as one table.

    The columnar counterpart to `all_features()`: same filters, but the
    result is a `pyarrow.Table` / `pandas.DataFrame` / `polars.DataFrame`
    instead of a stream of `Feature` objects, and no `Feature` is
    constructed at any layer.

    For "give me every exon as a dataframe", the row-by-row API pays a
    Python object per row -- millions of them on a whole-genome corpus --
    and that allocation dominates everything else. This is one query and
    one hand-off.

    Args:
        featuretype: One featuretype, or a list of them. `None` is
            everything.
        format: `"arrow"` (default, zero-copy), `"df"` for pandas, or
            `"polars"`.
        limit: Restrict to a genomic region, as `"seqid:start-end"` or a
            `(seqid, start, end)` tuple.
        strand: Restrict to `"+"`, `"-"` or `"."`.
        order_by: Column to sort by -- the same whitelist `all_features`
            accepts.
        reverse: Sort descending.
        completely_within: With `limit`, require full containment.

    Returns:
        A table in the shape named by `format`, with one row per feature and the columns of the `features` table (id, seqid, source, featuretype, start, end, score, strand, frame, file_order).

    Raises:
        ValueError: `order_by` names something outside the whitelist, or
            `format` is not one of the three.
        ImportError: `format="df"`/`"polars"` without that package.

    Example:
        ```python
        exons = db.to_table("exon", format="arrow")
        df = db.to_table(["exon", "CDS"], format="df", limit="chr1:1-10000")
        ```

    Note:
        Attributes are not included: they are a long-form table, so
        flattening them would either invent a column per key or collapse
        multi-valued keys. Query the `attributes` table with `execute()`
        when you need them.
    """
    self._require_open("to_table")
    sql, params = self._build_scan_sql(
        base_where=[],
        base_params=[],
        limit=limit,
        strand=strand,
        featuretype=featuretype,
        order_by=order_by,
        reverse=reverse,
        completely_within=completely_within,
    )
    return self._materialize_batched(sql, params, format=format)

region_batched

region_batched(regions, featuretype: str | list[str] | None = None, completely_within: bool = False, format: str = 'arrow', explode_segments: bool = False, on_invalid: str = 'raise')

Bulk overlap query. Performs a SINGLE spatial JOIN between every input region and the features table, returning a column-oriented result that maps each query (query_idx) back to its overlapping features.

query_idx indexes regions as you passed it. That is the whole contract of the column -- it is how a caller reassembles per-query groups without re-issuing N queries.

Parameters

regions : iterable of (seqid, start, end) | str | Feature Each item is normalized through _normalize_region_args; the same four input shapes accepted by region() are accepted here. featuretype : str | list[str] | None Optional features.featuretype filter applied to all regions. completely_within : bool If True, only features fully contained in the region are returned (default False — overlap is sufficient). format : "arrow" | "df" | "polars" Return shape (default "arrow"). on_invalid : {"raise", "skip"} What to do with an item that does not normalize to a region. "raise" (default) raises ValueError naming the offending position and value. "skip" drops it while leaving every other item's query_idx at its position in regions, so the gap is visible rather than silently closed up.

bool

Emit one row per physical INPUT LINE rather than one per logical feature, adding a seg_idx column. A discontinuous feature then contributes a row per segment, each with its own coordinates, score and phase -- which is what a caller writing a coverage track or exporting to a line-oriented format actually needs.

Offered on the tabular APIs only, never on region() / children() / all_features(): those must keep yielding Feature objects, and letting FeatureSegment rows leak into them would corrupt legacy consumers.

Result columns are query_idx, query_seqid, query_start, query_end, id, seqid, source, featuretype, start, end, score, strand, frame, file_order -- plus seg_idx when explode_segments.

Source code in python/gffbase/interface.py

def region_batched(
    self,
    regions,
    featuretype: str | list[str] | None = None,
    completely_within: bool = False,
    format: str = "arrow",
    explode_segments: bool = False,
    on_invalid: str = "raise",
):
    """Bulk overlap query. Performs a SINGLE spatial JOIN between every
    input region and the features table, returning a column-oriented
    result that maps each query (`query_idx`) back to its overlapping
    features.

    `query_idx` indexes `regions` as you passed it. That is the whole
    contract of the column -- it is how a caller reassembles per-query
    groups without re-issuing N queries.

    Parameters
    ----------
    regions : iterable of (seqid, start, end) | str | Feature
        Each item is normalized through `_normalize_region_args`; the
        same four input shapes accepted by `region()` are accepted here.
    featuretype : str | list[str] | None
        Optional `features.featuretype` filter applied to all regions.
    completely_within : bool
        If True, only features fully contained in the region are returned
        (default False — overlap is sufficient).
    format : "arrow" | "df" | "polars"
        Return shape (default `"arrow"`).
    on_invalid : {"raise", "skip"}
        What to do with an item that does not normalize to a region.
        ``"raise"`` (default) raises `ValueError` naming the offending
        position and value. ``"skip"`` drops it while leaving every other
        item's `query_idx` at its position in `regions`, so the gap is
        visible rather than silently closed up.

    explode_segments : bool
        Emit one row per physical INPUT LINE rather than one per logical
        feature, adding a `seg_idx` column. A discontinuous feature then
        contributes a row per segment, each with its own coordinates,
        score and phase -- which is what a caller writing a coverage track
        or exporting to a line-oriented format actually needs.

        Offered on the tabular APIs only, never on `region()` /
        `children()` / `all_features()`: those must keep yielding
        `Feature` objects, and letting `FeatureSegment` rows leak into
        them would corrupt legacy consumers.

    Result columns are query_idx, query_seqid, query_start,
    query_end, id, seqid, source, featuretype, start, end, score,
    strand, frame, file_order -- plus seg_idx when `explode_segments`.
    """
    self._require_open("region_batched")
    if on_invalid not in {"raise", "skip"}:
        raise ValueError(f"on_invalid must be 'raise' or 'skip'; got {on_invalid!r}")

    # `query_idx` carries the item's position in `regions`, NOT its
    # position among the ones that survived normalization. Those two used
    # to be the same expression -- `range(len(rows))` over the filtered
    # list -- so a single unparseable region silently shifted every later
    # query's index by one, and the caller mapped whole result groups onto
    # the wrong input. Nothing raised, and the answer stayed plausible.
    rows = []
    for idx, r in enumerate(regions):
        if isinstance(r, tuple) and len(r) == 3 and all(v is not None for v in r):
            seqid, rs, re_ = r[0], int(r[1]), int(r[2])
        else:
            seqid, rs, re_ = self._normalize_region_args(r, None, None, None)
        if seqid is None or rs is None or re_ is None:
            if on_invalid == "raise":
                raise ValueError(
                    f"regions[{idx}] = {r!r} does not describe a region "
                    "(need a 'seqid:start-end' string, a (seqid, start, end) "
                    "tuple, or a Feature). Pass on_invalid='skip' to drop it "
                    "and keep the remaining query_idx values aligned to the input."
                )
            continue
        rows.append((idx, seqid, int(rs), int(re_)))
    if not rows:
        return self._empty_region_batched(format, explode_segments)

    import pyarrow as pa

    regions_table = pa.table(
        {
            "query_idx": [r[0] for r in rows],
            "query_seqid": [r[1] for r in rows],
            "query_start": [r[2] for r in rows],
            "query_end": [r[3] for r in rows],
        }
    )
    self.conn.register("__staging_regions", regions_table)
    try:
        # The R-tree path uses the seqid_y-encoded envelope so that
        # DuckDB's spatial index segregates chromosomes.
        # B-tree fallback is identical SQL minus the ST_Intersects
        # predicate.
        ft_where, ft_params = self._featuretype_filter(featuretype, qualifier="f")
        ft_clause = (" AND " + " AND ".join(ft_where)) if ft_where else ""
        within_clause = (
            ' AND f.start >= q.query_start AND f."end" <= q.query_end'
            if completely_within
            else ""
        )
        # Same envelope-superset recheck as `_segment_overlap`, but
        # correlated against the staged query columns rather than bound
        # parameters. Empty -- and so byte-identical to v1 -- unless this
        # database actually holds a multipart feature. `sg`, not `s`: the
        # R-tree arm already binds `s` to the seqid lookup.
        # Exploding replaces the logical projection with `segments_all`,
        # which yields exactly one row per physical input line. The join to
        # `features` stays -- it is what carries the R-tree predicate -- so
        # the spatial index is still doing the selection and the extra join
        # only expands the rows it found.
        src = "sa" if explode_segments else "f"
        seg_cols = ", sa.seg_idx" if explode_segments else ""
        # `(query_idx, start)` alone is not a total order: two features
        # sharing a start came back in whatever order the join happened to
        # produce, which differed between the logical and exploded forms of
        # the same query. `id` makes it deterministic, and `seg_idx` keeps a
        # feature's own lines in file order.
        seg_order = ", sa.seg_idx" if explode_segments else ""
        explode_join = ""
        if explode_segments:
            bounds = (
                'sa.start >= q.query_start AND sa."end" <= q.query_end'
                if completely_within
                else 'sa.start <= q.query_end AND sa."end" >= q.query_start'
            )
            explode_join = " JOIN segments_all sa ON sa.feature_id = f.id AND " + bounds

        seg_clause = ""
        if self._n_multipart and not completely_within and not explode_segments:
            seg_clause = (
                " AND (f.n_segments = 1 OR EXISTS ("
                "SELECT 1 FROM segments sg WHERE sg.feature_id = f.id "
                'AND sg.start <= q.query_end AND sg."end" >= q.query_start))'
            )

        if self._rtree_built and self._seqid_y_map:
            # Inline the seqid → seqid_y map as a small VALUES table so
            # the JOIN can prune by chromosome inside the R-tree.
            values_pairs = ",".join(
                f"(?, {y})" for y in (self._seqid_y_map[s] for s in self._seqid_y_map)
            )
            seqid_y_params = list(self._seqid_y_map.keys())
            sql = f"""
                WITH seqid_lookup(seqid, seqid_y) AS (
                    VALUES {values_pairs}
                )
                SELECT
                    q.query_idx, q.query_seqid, q.query_start, q.query_end,
                    f.id, {src}.seqid, {src}.source, {src}.featuretype,
                    {src}.start, {src}."end" AS "end",
                    {src}.score, {src}.strand, {src}.frame, {src}.file_order{seg_cols}
                FROM __staging_regions q
                JOIN seqid_lookup s ON s.seqid = q.query_seqid
                JOIN features f
                  ON f.seqid = q.query_seqid
                  AND ST_Intersects(
                      f.bbox,
                      ST_MakeEnvelope(q.query_start, s.seqid_y,
                                      q.query_end,   s.seqid_y + 1)){explode_join}
                WHERE 1=1{within_clause}{seg_clause}{ft_clause}
                ORDER BY q.query_idx, {src}.start, f.id{seg_order}
            """
            params = list(seqid_y_params) + ft_params
        else:
            sql = f"""
                SELECT
                    q.query_idx, q.query_seqid, q.query_start, q.query_end,
                    f.id, {src}.seqid, {src}.source, {src}.featuretype,
                    {src}.start, {src}."end" AS "end",
                    {src}.score, {src}.strand, {src}.frame, {src}.file_order{seg_cols}
                FROM __staging_regions q
                JOIN features f
                  ON f.seqid = q.query_seqid
                  AND f.start <= q.query_end
                  AND f."end"  >= q.query_start{explode_join}
                WHERE 1=1{within_clause}{seg_clause}{ft_clause}
                ORDER BY q.query_idx, {src}.start, f.id{seg_order}
            """
            params = list(ft_params)

        return self._materialize_batched(sql, params, format=format)
    finally:
        try:
            self.conn.unregister("__staging_regions")
        except Exception:
            pass

children

children(id, level: int | None = None, featuretype: str | list[str] | None = None, order_by=None, reverse: bool = False, limit=None, completely_within: bool = False) -> Iterator[Feature]

Descendants of a feature, in hierarchy order.

Parameters

id : str or Feature The anchor. A Feature built by hand has id is None and is rejected rather than silently matching nothing. level : int, optional Only descendants exactly this many steps down -- 1 is direct children. None (default) walks the whole subtree. featuretype : str or list of str, optional order_by : str or sequence of str, optional A whitelisted sort key; anything else raises ValueError naming the accepted set. See docs/security/2026-sql-injection.md for why this is not a free-text field. reverse : bool Applies to every key, not just the last. limit : str or tuple, optional Restrict to a genomic interval, as region() accepts it. completely_within : bool With limit, require containment rather than overlap.

Reads from the materialized closure table when level is within the database's recorded depth, and falls back to a recursive CTE when a deeper walk is asked for. Both are cycle-safe: a Parent graph containing a loop is traversed once, not until the depth budget runs out.

Source code in python/gffbase/interface.py

def children(
    self,
    id,
    level: int | None = None,
    featuretype: str | list[str] | None = None,
    order_by=None,
    reverse: bool = False,
    limit=None,
    completely_within: bool = False,
) -> Iterator[Feature]:
    """Descendants of a feature, in hierarchy order.

    Parameters
    ----------
    id : str or Feature
        The anchor. A `Feature` built by hand has `id is None` and is
        rejected rather than silently matching nothing.
    level : int, optional
        Only descendants exactly this many steps down -- 1 is direct
        children. None (default) walks the whole subtree.
    featuretype : str or list of str, optional
    order_by : str or sequence of str, optional
        A whitelisted sort key; anything else raises `ValueError` naming
        the accepted set. See `docs/security/2026-sql-injection.md` for
        why this is not a free-text field.
    reverse : bool
        Applies to every key, not just the last.
    limit : str or tuple, optional
        Restrict to a genomic interval, as `region()` accepts it.
    completely_within : bool
        With `limit`, require containment rather than overlap.

    Reads from the materialized closure table when `level` is within the
    database's recorded depth, and falls back to a recursive CTE when a
    deeper walk is asked for. Both are cycle-safe: a `Parent` graph
    containing a loop is traversed once, not until the depth budget runs
    out.
    """
    target_id = _require_feature_id(id)
    yield from self._relation_query(
        target_id,
        level,
        featuretype,
        order_by,
        reverse,
        limit,
        completely_within,
        direction="children",
    )

parents

parents(id, level: int | None = None, featuretype: str | list[str] | None = None, order_by=None, reverse: bool = False, completely_within: bool = False, limit=None) -> Iterator[Feature]

Ancestors of a feature, in hierarchy order.

The mirror of :meth:children, taking the same arguments and making the same routing decision. level=1 is direct parents.

GFF3 permits a feature to name several Parents, so the hierarchy is a DAG rather than a tree and one ancestor can be reachable by more than one path. Each is returned once.

Source code in python/gffbase/interface.py

def parents(
    self,
    id,
    level: int | None = None,
    featuretype: str | list[str] | None = None,
    order_by=None,
    reverse: bool = False,
    completely_within: bool = False,
    limit=None,
) -> Iterator[Feature]:
    """Ancestors of a feature, in hierarchy order.

    The mirror of :meth:`children`, taking the same arguments and making
    the same routing decision. `level=1` is direct parents.

    GFF3 permits a feature to name several `Parent`s, so the hierarchy is
    a DAG rather than a tree and one ancestor can be reachable by more
    than one path. Each is returned once.
    """
    target_id = _require_feature_id(id)
    yield from self._relation_query(
        target_id,
        level,
        featuretype,
        order_by,
        reverse,
        limit,
        completely_within,
        direction="parents",
    )

children_batched

children_batched(feature_ids, level: int | None = None, featuretype: str | list[str] | None = None, format: str = 'arrow', explode_segments: bool = False)

Bulk children lookup. Returns the descendants of ALL feature_ids in a single vectorized DuckDB query.

Parameters

feature_ids : iterable of str | Feature Anchors. May contain Feature objects or raw ID strings. level : int | None None → all descendants (closure cache when materialized, otherwise dynamic CTE). Integer → exact-depth point lookup. featuretype : str | list[str] | None Optional filter on features.featuretype. format : "arrow" | "df" | "polars" Return shape (default "arrow" — pyarrow.Table).

bool

Emit one row per physical INPUT LINE rather than one per logical feature, adding a seg_idx column. Tabular APIs only -- see region_batched.

Result columns are anchor (the parent ID supplied), descendant_id, seqid, source, featuretype, start, end, score, strand, frame, file_order, depth -- plus seg_idx when explode_segments.

Source code in python/gffbase/interface.py

def children_batched(
    self,
    feature_ids,
    level: int | None = None,
    featuretype: str | list[str] | None = None,
    format: str = "arrow",
    explode_segments: bool = False,
):
    """Bulk children lookup. Returns the descendants of ALL `feature_ids`
    in a single vectorized DuckDB query.

    Parameters
    ----------
    feature_ids : iterable of str | Feature
        Anchors. May contain `Feature` objects or raw ID strings.
    level : int | None
        None → all descendants (closure cache when materialized,
        otherwise dynamic CTE). Integer → exact-depth point lookup.
    featuretype : str | list[str] | None
        Optional filter on `features.featuretype`.
    format : "arrow" | "df" | "polars"
        Return shape (default `"arrow"` — `pyarrow.Table`).

    explode_segments : bool
        Emit one row per physical INPUT LINE rather than one per logical
        feature, adding a `seg_idx` column. Tabular APIs only -- see
        `region_batched`.

    Result columns are anchor (the parent ID supplied),
    descendant_id, seqid, source, featuretype, start, end, score,
    strand, frame, file_order, depth -- plus seg_idx when
    `explode_segments`.
    """
    return self._batched_relation(
        feature_ids,
        level=level,
        featuretype=featuretype,
        direction="children",
        format=format,
        explode_segments=explode_segments,
    )

parents_batched

parents_batched(feature_ids, level: int | None = None, featuretype: str | list[str] | None = None, format: str = 'arrow', explode_segments: bool = False)

Bulk parents lookup. Mirrors children_batched but walks the closure / edges in the reverse direction.

Source code in python/gffbase/interface.py

def parents_batched(
    self,
    feature_ids,
    level: int | None = None,
    featuretype: str | list[str] | None = None,
    format: str = "arrow",
    explode_segments: bool = False,
):
    """Bulk parents lookup. Mirrors `children_batched` but walks the
    closure / edges in the reverse direction."""
    return self._batched_relation(
        feature_ids,
        level=level,
        featuretype=featuretype,
        direction="parents",
        format=format,
        explode_segments=explode_segments,
    )

delete

delete(features: FeatureLike | Iterable[FeatureLike], make_backup: bool = True, **kwargs) -> FeatureDB

Delete features, and everything that referenced them.

Removes the rows from features, attributes, segments and edges, then rebuilds the transitive closure so no path through a deleted node survives -- deleting a transcript really does remove its exons from its gene's descendants.

Parameters:

• features (FeatureLike | Iterable[FeatureLike]) –

  A feature id, a Feature, or an iterable of either.

• make_backup (bool, default: True ) –

  Accepted for gffutils compatibility and currently ignored; no .bak is written.

Returns:

• FeatureDB –

  self, so calls can be chained.

Raises:

• ReadOnlyError –

  The database was opened with read_only=True.

Example

db.delete("transcript_1")
db.delete([f.id for f in db.features_of_type("CDS")])

Source code in python/gffbase/interface.py

def delete(
    self, features: FeatureLike | Iterable[FeatureLike], make_backup: bool = True, **kwargs
) -> FeatureDB:
    """Delete features, and everything that referenced them.

    Removes the rows from `features`, `attributes`, `segments` and `edges`,
    then rebuilds the transitive closure so no path through a deleted node
    survives -- deleting a transcript really does remove its exons from its
    gene's descendants.

    Args:
        features: A feature id, a `Feature`, or an iterable of either.
        make_backup: Accepted for gffutils compatibility and currently
            ignored; no `.bak` is written.

    Returns:
        `self`, so calls can be chained.

    Raises:
        ReadOnlyError: The database was opened with `read_only=True`.

    Example:
        ```python
        db.delete("transcript_1")
        db.delete([f.id for f in db.features_of_type("CDS")])
        ```
    """
    self._require_writable("delete")
    ids = self._coerce_ids(features)
    if not ids:
        return self
    placeholders = ",".join("?" * len(ids))
    self.conn.execute(f"DELETE FROM features WHERE id IN ({placeholders})", ids)
    self.conn.execute(f"DELETE FROM attributes WHERE feature_id IN ({placeholders})", ids)
    # `segments` too. Without this the segment rows outlive their feature,
    # and `segments_all` joins them straight back -- so every physical-level
    # consumer (`export_sqlite`, the validator, `explode_segments`) would
    # keep reporting lines of a feature the caller deleted.
    self.conn.execute(f"DELETE FROM segments WHERE feature_id IN ({placeholders})", ids)
    self.conn.execute(
        f"DELETE FROM edges WHERE parent IN ({placeholders}) OR child IN ({placeholders})",
        ids + ids,
    )
    # Rebuild the closure rather than deleting the rows that mention these
    # ids. Deleting only those rows leaves behind every *transitive* row
    # that merely ROUTED THROUGH a deleted node: remove the mRNA from
    # gene -> mRNA -> exon and the depth-2 `gene -> exon` row survives,
    # because it names neither the mRNA as ancestor nor as descendant. So
    # `children(gene, level=None)` kept yielding exons of a transcript
    # that no longer existed. `update()` already rebuilds for the same
    # reason; deletion needs it at least as much.
    self._rebuild_closure()
    self._refresh_depth_meta()
    return self

update

update(data: Iterable, make_backup: bool = True, **kwargs) -> FeatureDB

Add features to an existing database.

Appends to features, attributes and edges, then rebuilds the transitive closure and re-reads the corpus statistics the relational dispatcher routes on.

Parameters:

• data (Iterable) –

  An iterable of Feature or ParsedFeature objects, or another FeatureDB whose features are copied in.

• make_backup (bool, default: True ) –

  Accepted for gffutils compatibility and currently ignored; no .bak is written.

Returns:

• FeatureDB –

  self, so calls can be chained.

Raises:

• ReadOnlyError –

  The database was opened with read_only=True.

• TypeError –

  An item is neither a Feature nor a ParsedFeature.

Example

introns = list(db.create_introns())
db.update(introns)

Source code in python/gffbase/interface.py

def update(self, data: Iterable, make_backup: bool = True, **kwargs) -> FeatureDB:
    """Add features to an existing database.

    Appends to `features`, `attributes` and `edges`, then rebuilds the
    transitive closure and re-reads the corpus statistics the relational
    dispatcher routes on.

    Args:
        data: An iterable of `Feature` or `ParsedFeature` objects, or
            another `FeatureDB` whose features are copied in.
        make_backup: Accepted for gffutils compatibility and currently
            ignored; no `.bak` is written.

    Returns:
        `self`, so calls can be chained.

    Raises:
        ReadOnlyError: The database was opened with `read_only=True`.
        TypeError: An item is neither a `Feature` nor a `ParsedFeature`.

    Example:
        ```python
        introns = list(db.create_introns())
        db.update(introns)
        ```
    """
    self._require_writable("update")
    # Minimal update: accept an iterable of Feature objects and
    # append them to features + attributes + edges, then refresh closure.
    from gffbase.feature import ParsedFeature
    from gffbase.ingest import _ArrowBatchBuilder

    # The builder needs the seqid_to_y dict so it can stamp
    # seqid_y (and bbox, when the R-tree is live) inline. Reuse the map
    # the FeatureDB already loaded from `seqid_map`.
    builder = _ArrowBatchBuilder(
        self._seqid_y_map,
        has_spatial=bool(self._rtree_built),
    )
    order = scalar_or(self.conn, "SELECT COALESCE(MAX(file_order), 0) FROM features", 0)

    if isinstance(data, FeatureDB):
        data = list(data.all_features())
    for feat in data:
        order += 1
        if isinstance(feat, Feature):
            blob = (
                feat._attributes_blob
                if feat._attributes_blob is not None
                else feat._format_attributes().encode("utf-8")
            )
            pairs = [(k, v, i) for k, vs in feat.attributes.items() for i, v in enumerate(vs)]
            pf = ParsedFeature(
                seqid=feat.seqid,
                source=feat.source,
                featuretype=feat.featuretype,
                start=feat.start,
                end=feat.end,
                score=feat.score,
                strand=feat.strand,
                frame=feat.frame,
                attributes_blob=blob,
                attributes_pairs=pairs,
                extra=list(feat.extra),
            )
            fid = feat.id or f"{feat.featuretype}_{order}"
        elif isinstance(feat, ParsedFeature):
            pf = feat
            fid = next(
                (v for k, v, _ in pf.attributes_pairs if k == "ID"),
                f"{pf.featuretype}_{order}",
            )
        else:
            raise TypeError(f"update() does not accept {type(feat)!r}")
        builder.append(fid, pf, order)
    builder.flush_into(self.conn)
    # Rebuild from edges (cheap on small updates), then re-read the
    # statistics the dispatcher routes on -- an update can deepen the
    # hierarchy or introduce the first multipart feature, and both were
    # previously left at whatever they were when the handle opened.
    self._rebuild_closure()
    self._refresh_depth_meta()
    return self

add_relation

add_relation(parent, child, level: int = 1, parent_func=None, child_func=None) -> FeatureDB

Link one parent to one child.

parent_func / child_func receive (parent, child), and whatever they RETURN is written back to the database. That is upstream's contract and the reason assign_child returns the child: the callback exists to edit an attribute (Parent=) that then has to be persisted. Previously the return value was discarded and nothing was written, so child_func=assign_child set an attribute on a throwaway object.

A string id is resolved to a Feature first, so callbacks fire whether the caller passed objects or ids -- they used to be skipped silently for ids.

Source code in python/gffbase/interface.py

def add_relation(
    self, parent, child, level: int = 1, parent_func=None, child_func=None
) -> FeatureDB:
    """Link one parent to one child.

    `parent_func` / `child_func` receive `(parent, child)`, and whatever
    they RETURN is written back to the database. That is upstream's
    contract and the reason `assign_child` returns the child: the callback
    exists to edit an attribute (`Parent=`) that then has to be persisted.
    Previously the return value was discarded and nothing was written, so
    `child_func=assign_child` set an attribute on a throwaway object.

    A string id is resolved to a `Feature` first, so callbacks fire whether
    the caller passed objects or ids -- they used to be skipped silently
    for ids.
    """
    return self.add_relations(
        [(parent, child)], level=level, parent_func=parent_func, child_func=child_func
    )

add_relations

add_relations(pairs, level: int = 1, parent_func=None, child_func=None) -> FeatureDB

add_relation for many pairs, with ONE closure rebuild.

Deriving the closure costs a recursive CTE over every edge, so doing it per pair makes a bulk operation quadratic. merge_all links every component of every merged feature and is the caller that made this necessary.

Source code in python/gffbase/interface.py

def add_relations(self, pairs, level: int = 1, parent_func=None, child_func=None) -> FeatureDB:
    """`add_relation` for many pairs, with ONE closure rebuild.

    Deriving the closure costs a recursive CTE over every edge, so doing it
    per pair makes a bulk operation quadratic. `merge_all` links every
    component of every merged feature and is the caller that made this
    necessary.
    """
    self._require_writable("add_relations")
    pairs = list(pairs)
    if not pairs:
        return self

    edges: list[list[str]] = []
    touched: dict[str, Feature] = {}
    for parent, child in pairs:
        parent = self[parent] if isinstance(parent, str) else parent
        child = self[child] if isinstance(child, str) else child
        edges.append([_require_feature_id(parent), _require_feature_id(child)])
        if parent_func is not None:
            updated = parent_func(parent, child)
            if updated is not None:
                touched[_require_feature_id(updated)] = updated
        if child_func is not None:
            updated = child_func(parent, child)
            if updated is not None:
                touched[_require_feature_id(updated)] = updated

    self.conn.executemany("INSERT INTO edges(parent, child) VALUES (?, ?)", edges)
    for feature in touched.values():
        self._write_back(feature)

    self._rebuild_closure()
    self._refresh_depth_meta()
    return self

interfeatures

interfeatures(features, new_featuretype=None, merge_attributes: bool = True, numeric_sort: bool = False, dialect=None, attribute_func=None, update_attributes=None)

The gaps between consecutive features.

attribute_func takes ONE argument -- an attribute mapping -- and returns one, matching the oracle. It is applied to each neighbour's attributes before they are merged, not to the merged result. The previous three-argument (prev, cur, attrs) form meant any gffutils caller passing a callback got a TypeError.

Source code in python/gffbase/interface.py

def interfeatures(
    self,
    features,
    new_featuretype=None,
    merge_attributes: bool = True,
    numeric_sort: bool = False,
    dialect=None,
    attribute_func=None,
    update_attributes=None,
):
    """The gaps between consecutive features.

    `attribute_func` takes ONE argument -- an attribute mapping -- and
    returns one, matching the oracle. It is applied to each neighbour's
    attributes *before* they are merged, not to the merged result. The
    previous three-argument `(prev, cur, attrs)` form meant any gffutils
    caller passing a callback got a `TypeError`.
    """
    from gffbase.helpers import merge_attributes as _merge_attrs

    if attribute_func is None:

        def attribute_func(a):
            return a

    feats = _with_coordinates(features)
    last = None
    for cur in feats:
        if last is None:
            last = cur
            continue
        if cur.seqid != last.seqid:
            # A gap between two different sequences is not a gap. Without
            # this the pair produced a feature spanning nothing, stamped
            # with the previous sequence's name.
            last = cur
            continue

        new_start = last.end + 1
        new_end = cur.start - 1
        if new_end < new_start:
            last = cur
            continue

        if merge_attributes:
            attrs = _merge_attrs(
                attribute_func(dict(last.attributes)),
                attribute_func(dict(cur.attributes)),
                numeric_sort=numeric_sort,
            )
        else:
            attrs = {}
        if update_attributes:
            attrs.update(update_attributes)

        # A feature may not carry several IDs, so a merged pair's two IDs
        # become one hyphenated id rather than a multi-valued attribute.
        if len(attrs.get("ID", [])) > 1:
            attrs["ID"] = ["-".join(attrs["ID"])]

        yield Feature(
            seqid=last.seqid,
            source=self.derived_source,
            featuretype=(
                new_featuretype
                if new_featuretype is not None
                else f"inter_{last.featuretype}_{cur.featuretype}"
            ),
            start=new_start,
            end=new_end,
            score=".",
            # Where the flanks disagree the gap has no orientation. This
            # used to inherit the left flank's strand unconditionally.
            strand=cur.strand if last.strand == cur.strand else ".",
            attributes=attrs,
            dialect=dialect or self.dialect,
        )
        last = cur

merge

merge(features, merge_criteria=None, multiline: bool = False)

Collapse runs of features that satisfy every criterion.

Consumes features in the order given. That is the oracle's contract and it matters: merge_all supplies a specific merge_order, and re-sorting here (which this used to do) silently discarded it.

A feature that merged with nothing is yielded unchanged with children set to no_children, so a caller can tell a real merge from a pass-through by truthiness. merge_all depends on exactly that.

Source code in python/gffbase/interface.py

def merge(self, features, merge_criteria=None, multiline: bool = False):
    """Collapse runs of features that satisfy every criterion.

    Consumes `features` **in the order given**. That is the oracle's
    contract and it matters: `merge_all` supplies a specific
    `merge_order`, and re-sorting here (which this used to do) silently
    discarded it.

    A feature that merged with nothing is yielded unchanged with
    `children` set to `no_children`, so a caller can tell a real merge from
    a pass-through by truthiness. `merge_all` depends on exactly that.
    """
    from gffbase import merge_criteria as mc

    if merge_criteria is None:
        merge_criteria = (mc.seqid, mc.overlap_end_inclusive, mc.strand, mc.feature_type)
    elif not isinstance(merge_criteria, (list, tuple)):
        merge_criteria = [merge_criteria]

    accum: Feature | None = None
    components: list[Feature] = []
    last_id: str | None = None

    for f in _with_coordinates(features):
        if accum is None:
            # A feature that fails its own criteria can never merge with
            # anything, so pass it straight through rather than opening a
            # run with it. Without this pre-pass such a feature was
            # silently accumulated into the next run.
            if all(pred(f, f, components) for pred in merge_criteria):
                accum, components, last_id = f, [f], None
            else:
                yield _finalize_merge(f, no_children)
            continue

        if not components:
            # `accum` came from a previous run's tail and has not been
            # checked against its own criteria yet.
            if all(pred(accum, accum, components) for pred in merge_criteria):
                components.append(accum)
            else:
                yield _finalize_merge(accum, no_children)
                accum, last_id = f, None
                continue

        if not all(pred(accum, f, components) for pred in merge_criteria):
            yield _finalize_merge(accum, components)
            accum, components, last_id = f, [], None
            continue

        if len(components) == 1:
            # About to merge for real, so stop mutating the caller's
            # feature and take a copy with an identity of its own.
            accum = self._clone_for_merge(accum)
            if not last_id:
                last_id = self._next_autoincrement_id(accum.featuretype)
            accum.id = last_id
            accum.attributes["ID"] = last_id
        components.append(f)

        # Ambiguity flags: where the components disagree, say so rather
        # than silently keeping the first one's value.
        if f.seqid not in accum.seqid.split(","):
            accum.seqid += "," + f.seqid
        if f.strand != accum.strand:
            accum.strand = "."
        if f.frame != accum.frame:
            accum.frame = "."
        if f.featuretype != accum.featuretype:
            accum.featuretype = "sequence_feature"
        # Both ends, not just the far one: with a caller-chosen
        # `merge_order` the run is not necessarily start-sorted, and only
        # extending `end` quietly truncated the merged feature.
        if f.start < accum.start:
            accum.start = f.start
        if f.end > accum.end:
            accum.end = f.end

    if accum is not None:
        yield _finalize_merge(accum, components)

merge_all

merge_all(merge_order=('seqid', 'featuretype', 'strand', 'start'), merge_criteria=None, featuretypes_groups=(None,), exclude_components: bool = False) -> list[Feature]

Merge everything in the database and write the results back.

Three things were wrong here and all three were silent. The method returned every input feature including ones that merged with nothing, so a caller could not tell what had actually been merged; it persisted nothing, despite documenting that "the resulting records are added to the database"; and it accepted exclude_components and ignored it, so asking for the components to be removed did nothing at all.

Source code in python/gffbase/interface.py

def merge_all(
    self,
    merge_order=("seqid", "featuretype", "strand", "start"),
    merge_criteria=None,
    featuretypes_groups=(None,),
    exclude_components: bool = False,
) -> list[Feature]:
    """Merge everything in the database and **write the results back**.

    Three things were wrong here and all three were silent. The method
    returned every input feature including ones that merged with nothing,
    so a caller could not tell what had actually been merged; it persisted
    nothing, despite documenting that "the resulting records are added to
    the database"; and it accepted `exclude_components` and ignored it, so
    asking for the components to be removed did nothing at all.
    """
    if not len(featuretypes_groups):
        # An empty tuple used to mean "merge nothing" and return []. The
        # oracle reads it as "no featuretype filter".
        featuretypes_groups = (None,)

    result: list[Feature] = []
    for group in featuretypes_groups:
        merged_in_group = [
            merged
            for merged in self.merge(
                self.all_features(featuretype=group, order_by=merge_order),
                merge_criteria=merge_criteria,
            )
            if merged.children
        ]
        for merged in merged_in_group:
            self._insert(merged)
            result.append(merged)

        if exclude_components:
            self.delete(
                [c for merged in merged_in_group for c in merged.children],
                make_backup=False,
            )
        else:
            # One batched call, not one per child: `add_relation` re-derives
            # the closure each time, so the oracle's per-child loop would
            # make this O(children x full rebuild).
            self.add_relations(
                [(merged, child) for merged in merged_in_group for child in merged.children],
                level=1,
                child_func=assign_child,
            )
    return result

create_introns

create_introns(exon_featuretype: str = 'exon', grandparent_featuretype: str | None = 'gene', parent_featuretype: str | None = None, new_featuretype: str = 'intron', merge_attributes: bool = True, numeric_sort: bool = False) -> Iterator[Feature]

Introns, computed per transcript.

grandparent_featuretype="gene" descends one level first and computes the gaps within each transcript separately. Treating the gene as the direct anchor -- which this used to do -- pools the exons of every isoform into one sorted list, so the "introns" of a multi-isoform gene were computed across transcript boundaries and were not introns of anything.

Source code in python/gffbase/interface.py

def create_introns(
    self,
    exon_featuretype: str = "exon",
    grandparent_featuretype: str | None = "gene",
    parent_featuretype: str | None = None,
    new_featuretype: str = "intron",
    merge_attributes: bool = True,
    numeric_sort: bool = False,
) -> Iterator[Feature]:
    """Introns, computed **per transcript**.

    `grandparent_featuretype="gene"` descends one level first and computes
    the gaps within each transcript separately. Treating the gene as the
    direct anchor -- which this used to do -- pools the exons of every
    isoform into one sorted list, so the "introns" of a multi-isoform gene
    were computed across transcript boundaries and were not introns of
    anything.
    """
    for anchor in self._exon_anchors(grandparent_featuretype, parent_featuretype):
        exons = self.children(anchor, level=1, featuretype=exon_featuretype, order_by="start")
        yield from self.interfeatures(
            exons,
            new_featuretype=new_featuretype,
            merge_attributes=merge_attributes,
            numeric_sort=numeric_sort,
            dialect=self.dialect,
        )

create_splice_sites

create_splice_sites(exon_featuretype: str = 'exon', grandparent_featuretype: str | None = 'gene', parent_featuretype: str | None = None, merge_attributes: bool = True, numeric_sort: bool = False) -> Iterator[Feature]

The two-base splice sites flanking each intron.

A splice site is a dinucleotide -- GT at the donor, AG at the acceptor -- so these are 2 bp features, not the 1 bp ones this used to emit. They are typed by their position in the transcript rather than in the genome, so the left site of a minus-strand transcript is its 3' site.

Emission order is every left site, then every right site, matching the oracle; the intron's merged attributes are carried through with the ID prefixed by the featuretype so the two sites of one intron differ.

Source code in python/gffbase/interface.py

def create_splice_sites(
    self,
    exon_featuretype: str = "exon",
    grandparent_featuretype: str | None = "gene",
    parent_featuretype: str | None = None,
    merge_attributes: bool = True,
    numeric_sort: bool = False,
) -> Iterator[Feature]:
    """The two-base splice sites flanking each intron.

    A splice site is a dinucleotide -- GT at the donor, AG at the acceptor
    -- so these are 2 bp features, not the 1 bp ones this used to emit.
    They are typed by their position in the transcript rather than in the
    genome, so the left site of a minus-strand transcript is its 3' site.

    Emission order is every left site, then every right site, matching the
    oracle; the intron's merged attributes are carried through with the ID
    prefixed by the featuretype so the two sites of one intron differ.
    """
    for side in ("left", "right"):
        for anchor in self._exon_anchors(grandparent_featuretype, parent_featuretype):
            exons = self.children(
                anchor, level=1, featuretype=exon_featuretype, order_by="start"
            )
            if anchor.strand == "+":
                featuretype = (
                    "five_prime_cis_splice_site"
                    if side == "left"
                    else "three_prime_cis_splice_site"
                )
            elif anchor.strand == "-":
                featuretype = (
                    "three_prime_cis_splice_site"
                    if side == "left"
                    else "five_prime_cis_splice_site"
                )
            else:
                # No orientation, so neither end is 5' or 3'.
                featuretype = "splice_site"

            for site in self.interfeatures(
                exons,
                new_featuretype=featuretype,
                merge_attributes=merge_attributes,
                numeric_sort=numeric_sort,
                dialect=self.dialect,
            ):
                if side == "left":
                    site.end = site.start + 1
                else:
                    site.start = site.end - 1
                if site.attributes.get("ID"):
                    site.attributes["ID"] = [f"{featuretype}_{site.attributes['ID'][0]}"]
                yield site

children_bp

children_bp(feature: FeatureLike, child_featuretype: str = 'exon', merge: bool = False, merge_criteria: Sequence | None = None, **kwargs) -> int

Total base pairs covered by a feature's children.

Parameters:

• feature (FeatureLike) –

  The parent, as an id or a Feature.

• child_featuretype (str, default: 'exon' ) –

  Which children to measure.

• merge (bool, default: False ) –

  Merge overlapping children first, so shared bases are counted once. Without it, overlapping children double-count.

• merge_criteria (Sequence | None, default: None ) –

  Predicates controlling what may merge; see gffbase.merge_criteria. Defaults to same seqid, strand and featuretype with inclusive overlap.

Returns:

• int –

  The summed length in base pairs.

Raises:

• ValueError –

  The removed ignore_strand argument was passed.

• TypeError –

  Any other unexpected keyword argument.

Example

db.children_bp("transcript_1", child_featuretype="exon", merge=True)

Source code in python/gffbase/interface.py

def children_bp(
    self,
    feature: FeatureLike,
    child_featuretype: str = "exon",
    merge: bool = False,
    merge_criteria: Sequence | None = None,
    **kwargs,
) -> int:
    """Total base pairs covered by a feature's children.

    Args:
        feature: The parent, as an id or a `Feature`.
        child_featuretype: Which children to measure.
        merge: Merge overlapping children first, so shared bases are
            counted once. Without it, overlapping children double-count.
        merge_criteria: Predicates controlling what may merge; see
            `gffbase.merge_criteria`. Defaults to same seqid, strand and
            featuretype with inclusive overlap.

    Returns:
        The summed length in base pairs.

    Raises:
        ValueError: The removed `ignore_strand` argument was passed.
        TypeError: Any other unexpected keyword argument.

    Example:
        ```python
        db.children_bp("transcript_1", child_featuretype="exon", merge=True)
        ```
    """
    if kwargs:
        # Accepting and ignoring these was worse than refusing them:
        # `ignore_strand` was removed upstream precisely because it gave
        # the wrong answer, and silently dropping it returns a number that
        # looks right.
        if "ignore_strand" in kwargs:
            raise ValueError(
                "'ignore_strand' has been deprecated; please use merge_criteria to "
                "control how features should be merged. E.g., leave out the mc.strand "
                "criteria to ignore strand."
            )
        raise TypeError(f"children_bp() got unexpected keyword arguments {list(kwargs)}")

    kids = self.children(feature, featuretype=child_featuretype, order_by="start")
    if merge:
        kids = self.merge(kids, merge_criteria=merge_criteria)
    total = 0
    for k in kids:
        if k.start is not None and k.end is not None:
            total += k.end - k.start + 1
    return total

bed12

bed12(feature: FeatureLike, block_featuretype: Sequence[str] = ('exon',), thick_featuretype: Sequence[str] = ('CDS',), thin_featuretype: Sequence[str] | None = None, name_field: str = 'ID', color: str | None = None) -> str

Render a feature and its children as one BED12 line.

Parameters:

• feature (FeatureLike) –

  The parent, as an id or a Feature.

• block_featuretype (Sequence[str], default: ('exon',) ) –

  Child types that become BED blocks (exons).

• thick_featuretype (Sequence[str], default: ('CDS',) ) –

  Child types that define the thick region (coding sequence).

• thin_featuretype (Sequence[str] | None, default: None ) –

  Child types that define the thin region. When given, it is honoured rather than inferred.

• name_field (str, default: 'ID' ) –

  Attribute used for BED column 4. Falls back to the feature id when absent.

• color (str | None, default: None ) –

  RGB string for column 9, e.g. "255,0,0".

Returns:

• str –

  A tab-separated BED12 line, without a trailing newline.

Note

A feature with no thick children is emitted as entirely thick. blockSizes and blockStarts carry no trailing comma.

Example

print(db.bed12("transcript_1"))

Source code in python/gffbase/interface.py

def bed12(
    self,
    feature: FeatureLike,
    block_featuretype: Sequence[str] = ("exon",),
    thick_featuretype: Sequence[str] = ("CDS",),
    thin_featuretype: Sequence[str] | None = None,
    name_field: str = "ID",
    color: str | None = None,
) -> str:
    """Render a feature and its children as one BED12 line.

    Args:
        feature: The parent, as an id or a `Feature`.
        block_featuretype: Child types that become BED blocks (exons).
        thick_featuretype: Child types that define the thick region
            (coding sequence).
        thin_featuretype: Child types that define the thin region. When
            given, it is honoured rather than inferred.
        name_field: Attribute used for BED column 4. Falls back to the
            feature id when absent.
        color: RGB string for column 9, e.g. `"255,0,0"`.

    Returns:
        A tab-separated BED12 line, without a trailing newline.

    Note:
        A feature with no thick children is emitted as entirely thick.
        `blockSizes` and `blockStarts` carry no trailing comma.

    Example:
        ```python
        print(db.bed12("transcript_1"))
        ```
    """
    if thick_featuretype and thin_featuretype:
        raise ValueError("Can only specify one of `thick_featuretype` or `thin_featuretype`")
    if isinstance(feature, str):
        feature = self[feature]
    blocks = sorted(
        _with_coordinates(self.children(feature, featuretype=list(block_featuretype))),
        key=lambda f: (f.start, f.end),
    )
    if feature.start is None or feature.end is None:
        raise ValueError(
            f"cannot build a BED12 record for {feature.id!r}: "
            "feature has no start/end coordinates"
        )
    chrom_start = feature.start - 1
    chrom_end = feature.end

    if thin_featuretype:
        # The complement of `thick`: the caller names the UNtranslated
        # parts, and the thick span is what lies between them. Accepted and
        # silently ignored before, so `thin_featuretype=["UTR"]` returned a
        # record with the thick span covering the whole feature.
        thin = sorted(
            _with_coordinates(self.children(feature, featuretype=list(thin_featuretype))),
            key=lambda f: (f.start, f.end),
        )
        if thin:
            thick_start, thick_end = thin[0].end, thin[-1].start - 1
        else:
            thick_start, thick_end = feature.start, feature.end
    else:
        thick = sorted(
            _with_coordinates(self.children(feature, featuretype=list(thick_featuretype))),
            key=lambda f: (f.start, f.end),
        )
        if thick:
            thick_start, thick_end = thick[0].start - 1, thick[-1].end
        else:
            # No CDS: the oracle marks the whole feature thick, using its
            # 1-based start. Collapsing both to `chrom_start` -- which this
            # did -- renders an entirely thin feature, the opposite claim.
            thick_start, thick_end = feature.start, feature.end

    try:
        name_value = feature.attributes[name_field][0]
    except (KeyError, IndexError):
        name_value = "."
    score = feature.score if feature.score not in (".", "") else "0"
    # `.` is a legal BED strand and means "unstranded". Rewriting it to `+`
    # asserts an orientation the source did not have.
    strand = feature.strand if feature.strand in ("+", "-") else "."
    rgb = (color or "0,0,0").replace(" ", "").strip()
    # BED12 requires blockCount to equal the number of entries in
    # blockSizes and blockStarts, so a block with missing coordinates has
    # to drop out of all three together, not just the two lists.
    sized = [(b.start, b.end) for b in blocks if b.start is not None and b.end is not None]
    if not sized:
        # A feature with no block children is one block: itself.
        sized = [(feature.start, feature.end)]
    # The blocks must span the feature. BED12's blockStarts are offsets
    # from chromStart and the last block has to reach chromEnd, so a
    # feature whose children do not span it produces a structurally
    # invalid line -- one that names a range it does not cover.
    #
    # gffutils refuses this ("Start of first exon (%s) does not match start
    # of feature (%s)"), and it is right to: the usual cause is asking for
    # a `block_featuretype` the feature does not have all of, and silently
    # emitting a wrong line sends it downstream into a genome browser.
    first_start, last_end = sized[0][0], sized[-1][1]
    if first_start != feature.start:
        raise ValueError(
            f"Start of first exon ({first_start}) does not match start of "
            f"feature ({feature.start})"
        )
    if last_end != feature.end:
        raise ValueError(
            f"End of last exon ({last_end}) does not match end of feature ({feature.end})"
        )

    block_count = len(sized)
    # No trailing comma. UCSC tolerates one, but the oracle emits none and
    # a differential comparison sees every line as different.
    block_sizes = ",".join(str(end - start + 1) for start, end in sized)
    block_starts = ",".join(str((start - 1) - chrom_start) for start, _end in sized)
    return "\t".join(
        str(x)
        for x in (
            feature.seqid,
            chrom_start,
            chrom_end,
            name_value,
            score,
            strand,
            thick_start,
            thick_end,
            rgb,
            block_count,
            block_sizes,
            block_starts,
        )
    )

iter_by_parent_childs

iter_by_parent_childs(featuretype: str = 'gene', level: int | None = None, order_by: str | None = None, reverse: bool = False, completely_within: bool = False) -> Iterator[list[Feature]]

Group the database by parent, yielding one list per parent.

Parameters:

• featuretype (str, default: 'gene' ) –

  The parent featuretype to group by.

• level (int | None, default: None ) –

  How deep to collect children. None takes the whole subtree; 1 takes direct children only.

• order_by (str | None, default: None ) –

  Column to sort parents by -- see all_features.

• reverse (bool, default: False ) –

  Sort parents descending.

• completely_within (bool, default: False ) –

  Passed through to the child query.

Yields:

• list[Feature] –

  A list per parent, the parent first followed by its children.

Example

for group in db.iter_by_parent_childs("gene"):
    gene, children = group[0], group[1:]

Source code in python/gffbase/interface.py

def iter_by_parent_childs(
    self,
    featuretype: str = "gene",
    level: int | None = None,
    order_by: str | None = None,
    reverse: bool = False,
    completely_within: bool = False,
) -> Iterator[list[Feature]]:
    """Group the database by parent, yielding one list per parent.

    Args:
        featuretype: The parent featuretype to group by.
        level: How deep to collect children. `None` takes the whole
            subtree; `1` takes direct children only.
        order_by: Column to sort parents by -- see `all_features`.
        reverse: Sort parents descending.
        completely_within: Passed through to the child query.

    Yields:
        A list per parent, the parent first followed by its children.

    Example:
        ```python
        for group in db.iter_by_parent_childs("gene"):
            gene, children = group[0], group[1:]
        ```
    """
    for parent in self.features_of_type(featuretype, order_by=order_by, reverse=reverse):
        kids = list(
            self.children(
                parent,
                level=level,
                order_by=order_by,
                reverse=reverse,
                completely_within=completely_within,
            )
        )
        yield [parent, *kids]

attribute_search

attribute_search(text: str, featuretype=None) -> Iterator[Feature]

Features with an attribute VALUE matching text, case-insensitively.

text is SQL LIKE syntax, so % and _ are wildcards; a bare string matches as a substring, which is what a caller searching for an accession expects.

gffutils' CLI calls db.attribute_search(...), but no such method exists anywhere in gffutils -- only in that call site and in an obsolete test file -- so gffutils-cli search raises AttributeError on every invocation. This is a working implementation rather than a port of one.

The search is over the long-form attributes table, not over the raw column-9 blob, so it matches DECODED values: searching for a;b finds a feature whose source said a%3Bb.

Source code in python/gffbase/interface.py

def attribute_search(self, text: str, featuretype=None) -> Iterator[Feature]:
    """Features with an attribute VALUE matching `text`, case-insensitively.

    `text` is SQL `LIKE` syntax, so `%` and `_` are wildcards; a bare
    string matches as a substring, which is what a caller searching for an
    accession expects.

    gffutils' CLI calls `db.attribute_search(...)`, but no such method
    exists anywhere in gffutils -- only in that call site and in an
    obsolete test file -- so `gffutils-cli search` raises `AttributeError`
    on every invocation. This is a working implementation rather than a
    port of one.

    The search is over the long-form `attributes` table, not over the raw
    column-9 blob, so it matches DECODED values: searching for `a;b` finds
    a feature whose source said `a%3Bb`.
    """
    pattern = text if any(ch in text for ch in "%_") else f"%{text}%"
    params: list = [pattern]
    clause = ""
    if featuretype:
        if isinstance(featuretype, str):
            clause = " AND f.featuretype = ?"
            params.append(featuretype)
        else:
            types = list(featuretype)
            clause = f" AND f.featuretype IN ({','.join('?' * len(types))})"
            params.extend(types)

    sql = (
        f"SELECT {self._select_feature_aliased('f')} FROM features f "
        "WHERE EXISTS (SELECT 1 FROM attributes a WHERE a.feature_id = f.id "
        f"AND lower(a.value) LIKE lower(?)){clause} ORDER BY f.file_order"
    )
    yield from self._yield_features(sql, params)

execute

execute(query: str)

Execute arbitrary SQL. Returns DuckDB's relation cursor. SQLite-style queries against features_compat and relations_compat views are supported; see compat_views.sql.

Deliberately NOT guarded against read_only: this is the escape hatch, the SQL is the caller's, and DuckDB's own refusal names the statement it rejected, which is more use here than a generic message from us.

Source code in python/gffbase/interface.py

def execute(self, query: str):
    """Execute arbitrary SQL. Returns DuckDB's relation cursor.
    SQLite-style queries against ``features_compat`` and ``relations_compat``
    views are supported; see ``compat_views.sql``.

    Deliberately NOT guarded against ``read_only``: this is the escape
    hatch, the SQL is the caller's, and DuckDB's own refusal names the
    statement it rejected, which is more use here than a generic message
    from us.
    """
    self._require_open("execute")
    return self.conn.execute(query.rstrip(";"))

analyze

analyze() -> None

Refresh DuckDB's planner statistics for this database.

Worth running once after a large update(); the planner otherwise keeps costing queries against the shape the database had at ingest.

Raises:

• ReadOnlyError –

  The database was opened with read_only=True.

Source code in python/gffbase/interface.py

def analyze(self) -> None:
    """Refresh DuckDB's planner statistics for this database.

    Worth running once after a large `update()`; the planner otherwise
    keeps costing queries against the shape the database had at ingest.

    Raises:
        ReadOnlyError: The database was opened with `read_only=True`.
    """
    # ANALYZE writes statistics into the database, so it is a mutation.
    self._require_writable("analyze")
    self.conn.execute("ANALYZE")
    self._analyzed_flag = True

set_pragmas

set_pragmas(pragmas: dict) -> None

Apply DuckDB settings, ignoring pragmas that only SQLite has.

Legacy callers pass constants.default_pragmas -- synchronous, journal_mode, main.page_size, main.cache_size -- none of which DuckDB has. Those are skipped, which is what makes a gffutils script run here unchanged.

Both the name and the value used to be interpolated straight into the statement, and the whole loop body sat inside except duckdb.Error: continue. DuckDB executes trailing statements, so

db.set_pragmas({"threads": "1; DROP TABLE attributes"})

dropped the table -- and because the exception was swallowed, a payload that failed was silent too. Names are now matched against DuckDB's own settings catalog and values rendered as SQL literals, so nothing a caller supplies reaches the parser as syntax.

Matching the live catalog rather than a hardcoded list means the check tracks whatever DuckDB build is installed, instead of going stale and rejecting settings a newer version added.

Source code in python/gffbase/interface.py

def set_pragmas(self, pragmas: dict) -> None:
    """Apply DuckDB settings, ignoring pragmas that only SQLite has.

    Legacy callers pass `constants.default_pragmas` -- `synchronous`,
    `journal_mode`, `main.page_size`, `main.cache_size` -- none of which
    DuckDB has. Those are skipped, which is what makes a gffutils script
    run here unchanged.

    Both the name and the value used to be interpolated straight into the
    statement, and the whole loop body sat inside `except duckdb.Error:
    continue`. DuckDB executes trailing statements, so

        db.set_pragmas({"threads": "1; DROP TABLE attributes"})

    dropped the table -- and because the exception was swallowed, a
    payload that *failed* was silent too. Names are now matched against
    DuckDB's own settings catalog and values rendered as SQL literals, so
    nothing a caller supplies reaches the parser as syntax.

    Matching the live catalog rather than a hardcoded list means the check
    tracks whatever DuckDB build is installed, instead of going stale and
    rejecting settings a newer version added.
    """
    known = {
        row[0] for row in self.conn.execute("SELECT name FROM duckdb_settings()").fetchall()
    }
    for k, v in pragmas.items():
        name = str(k)
        if name not in known:
            # Not a DuckDB setting. Previously this was indistinguishable
            # from "DuckDB rejected the value"; now it is a decision.
            _log.debug("set_pragmas: skipping %r, not a DuckDB setting", name)
            continue
        # `name` is echoed from the catalog, so it cannot carry syntax.
        self.conn.execute(f"SET {name} = {_sql_literal(v)}")