Analysis gallery¶
Complete analyses, each one a real question with the code that answers it and
the output it actually produces. Every Python block on this page runs in CI
against the vendored tests/data/hierarchy.gff3 fixture, so the numbers below
are what the code prints — not what it printed once.
For per-method snippets see the Usage gallery; for corpus-specific workflows see the Cookbooks.
The setup every example below assumes:
from gffbase import create_db
db = create_db(HIERARCHY_PATH, str(TMP / "gallery.duckdb"), force=True)
Extract splice sites¶
Every donor and acceptor in the annotation — the input a splice-site model like Splam or OpenSpliceAI trains on. An intron is the gap between consecutive exons of the same transcript, so this is one sort and one pairwise walk.
for tx in db.features_of_type("mRNA"):
exons = sorted(db.children(tx.id, featuretype="exon"), key=lambda e: e.start)
for left, right in zip(exons, exons[1:]):
print(
f"{tx.id} donor {left.end + 1} acceptor {right.start - 1} "
f"intron {right.start - left.end - 1} bp"
)
At whole-genome scale, batch it
On a real annotation, fetch every exon in one query instead of looping:
db.children_batched(transcript_ids, featuretype="exon", format="arrow")
returns an Arrow table with an anchor column naming each row's
transcript, so you can group without re-querying. See
Machine learning workflows.
create_introns() gives you the intron features directly if you want them as
Feature objects rather than coordinates.
Pick the longest isoform per gene¶
The standard reduction to one transcript per gene — by exonic length, not
genomic span, which is what children_bp measures.
for gene in db.features_of_type("gene"):
isoforms = list(db.children(gene.id, featuretype="mRNA"))
longest = max(isoforms, key=lambda t: db.children_bp(t, child_featuretype="exon"))
span = db.children_bp(longest, child_featuretype="exon")
print(f"{gene.id}: {longest.id} ({span} exonic bp of {len(isoforms)} isoforms)")
For a curated answer rather than a computed one, MANE ships exactly one representative transcript per gene — see MANE.
Exon-count distribution¶
Single-exon genes behave differently from multi-exon ones in almost every analysis, so this histogram is worth looking at before trusting a result.
from collections import Counter
counts = Counter(
len(list(db.children(t.id, featuretype="exon")))
for t in db.features_of_type("mRNA")
)
for n_exons, n_tx in sorted(counts.items()):
print(f"{n_exons} exon(s): {n_tx} transcript(s)")
Coding fraction per transcript¶
How much of each transcript is protein-coding. A transcript with no CDS is a non-coding isoform, which this makes visible rather than silently averaging away.
for tx in db.features_of_type("mRNA"):
coding = db.children_bp(tx, child_featuretype="CDS")
exonic = db.children_bp(tx, child_featuretype="exon")
pct = 100 * coding / exonic if exonic else 0
print(f"{tx.id}: {coding}/{exonic} bp coding ({pct:.0f}%)")
Hand the whole annotation to pandas¶
When the analysis is columnar, skip Feature objects entirely.
to_table() takes the same filters as all_features() and
returns Arrow, pandas or polars.
Restrict to a locus with limit=, and get a dataframe instead:
region = db.to_table("exon", limit=("chr1", 100, 300), format="df")
print(len(region), "exons in chr1:100-300")
What is in an unfamiliar annotation?¶
Before any of the above, from the shell:
Feature counts by type, sequence count and names, and how the database was built. See Command line.
Where to go next¶
- Cookbooks — the same techniques against GENCODE, RefSeq and MANE, at whole-genome scale.
- Machine learning workflows — bulk extraction into PyTorch and Hugging Face with no per-row Python.
- Usage gallery — one snippet per public method.