bio-genome-intervals-gtf-gff-handling

What this skill does

## Version Compatibility

Reference examples tested with: gffutils 0.13+, gffread 0.12+, gtfparse 2.x, pyranges 0.1+ (or 1.0+ - see note), AGAT 1.4+.

Before using code patterns, verify installed versions match. If versions differ:
- CLI: `<tool> --version` then `<tool> --help` to confirm flags
- Python: `pip show <package>` then `help(module.function)` to check signatures

Two version landmines specific to this skill: (1) **gtfparse changed its return type** - older releases returned a pandas DataFrame, gtfparse >=2.x returns a **polars** DataFrame by default; pass `result_type='pandas'` before chaining pandas idioms (`.copy()`, boolean masks). (2) **pyranges has a major-version API split** - pyranges 0.x and the 1.0 rewrite differ in method names and attribute access; check `import pyranges; pyranges.__version__` before pasting code. gffutils stores **1-based** coordinates while pyranges stores **0-based** - their `start` fields differ by one by design. If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt rather than retrying.

# GTF/GFF Handling

**"Pull these features (or their sequences) out of my annotation, convert it, or find out why my counts are wrong."** -> Treat the file as a serialized gene-model tree: walk gene->transcript->exon/CDS, derive implicit features, and reconcile coordinate and namespace conventions before trusting any number.
- CLI: `gffread in.gtf -T -o out.gtf` (convert), `gffread -w tx.fa -g genome.fa in.gtf` (FASTA), `agat_convert_sp_gxf2gxf.pl` (sanitize)
- Python: `gffutils.create_db(...)` then `db.children(gene, featuretype='exon')` (tree query); `gtfparse.read_gtf(..., result_type='pandas')` / `pyranges.read_gtf(...)` (dataframe)

## The Single Most Important Modern Insight -- A GTF/GFF3 Is a Serialized Gene-Model Tree, Not a Table of Intervals

Almost every painful bug here comes from the file looking like a CSV while behaving like a tree, or from a coordinate/provenance mismatch the tools never warn about - and **every one of these failures is silent**: nothing throws, the wrong answer just propagates. Three load-bearing facts the tutorials skip:

1. **The coordinate conversion is asymmetric.** GTF/GFF3 are 1-based fully inclusive `[start, end]`; BED (and pyranges-internal) are 0-based half-open `[start-1, end)`. Convert to BED by **subtracting 1 from the start only - the end is unchanged** (the inclusive 1-based end and the exclusive 0-based end are the same integer). Doing `start-1` AND `end-1` shifts the feature one base left and is the classic over-correction: invisible in coverage/overlap, catastrophic in CDS translation (a one-base frameshift garbles the protein). gffutils keeps 1-based, pyranges stores 0-based, so their `start` fields differ by one *correctly* - never "fix" that discrepancy.

2. **The all-zero count matrix.** featureCounts/htseq-count match a read to a feature by **string equality on the chromosome name**, so `chr1` != `1` != `NC_000001.11` produces a perfectly well-formed matrix of **zeros with no error or warning** - the only signal is `~0%` assigned in the summary. Same bug one altitude up: gene-ID version suffixes (`ENSG00000223972.5` vs `ENSG00000223972`) silently drop rows on an annotation join. Audit every cross-file key (chromosomes between BAM/GTF/FASTA, gene IDs between GTF/count-matrix/annotation) by **set intersection, never by eye**, before any count or join.

3. **phase is not frame, and the stop codon is a 3-bp ghost.** Phase (column 8) is the strand-aware count of bases to trim from the segment's transcriptional 5' end to reach the next codon (0/1/2) - **recompute it on any CDS edit** (AGAT/gffread do; hand-editing coordinates without fixing phase frameshifts the translation). GTF (Ensembl/GENCODE) **excludes** the stop codon from CDS; GenBank/GFF3 often **include** it - so a CDS length off by exactly 3 nt (or a protein +/-1 stop) between two sources is a convention mismatch, not a bug.

## Tool Taxonomy

| Tool | Role | Mechanism | When |
|------|------|-----------|------|
| gffutils | Queryable gene-tree DB (Python) | builds a SQLite DB; `children`/`parents`/`region` traverse the hierarchy; keeps 1-based coords | walk gene->transcript->exon/CDS, derive introns, query by ID/coordinate |
| gffread | Converter + sequence extractor (CLI) | fast C++, genome-aware; one binary | GTF<->GFF3, extract transcript/CDS/protein FASTA, region filter |
| pyranges | Vectorized interval engine (Python) | PyRanges/pandas-like; stores 0-based half-open | overlap joins, set ops, dataframe-native interval work |
| gtfparse | GTF -> dataframe (Python) | one call explodes column 9 into attribute columns | quick column/filter work; NOT hierarchy-aware (flat table) |
| AGAT | GFF/GTF sanitizer (Perl CLI) | reconstructs the full tree; adds missing features, fixes IDs/phase, deflates attributes | a malformed/non-standard file - run FIRST, before parsing |

## Decision Tree by Scenario

| Scenario | Recommended | Why |
|----------|-------------|-----|
| Walk the gene/transcript/exon hierarchy, derive introns | gffutils `create_db` + `children`/`parents` | the hierarchy is the point; flat parsers lose it |
| Convert GTF<->GFF3 or extract transcript/CDS/protein FASTA | gffread (`-T`, `-w`/`-x`/`-y -g`) | genome-aware, knows the stop-codon convention |
| Quick column/filter on a clean modern GTF | gtfparse (`result_type='pandas'`) | one-call dataframe; verify return type first |
| Overlap/set ops, large in-memory interval joins | pyranges | vectorized; route arithmetic -> interval-arithmetic |
| File malformed: no `##gff-version`, missing gene/exon lines, dup IDs, mixed conventions | AGAT `agat_convert_sp_gxf2gxf.pl` first | sanitize once vs writing a brittle parser around it |
| GTF -> BED for bedtools | `start-1`, end unchanged (-> bed-file-basics) | the off-by-one boundary is where it bites |
| Counts came out all-zero or DE join dropped rows | intersect seqid / gene-ID namespaces | string-equality match; no error is emitted |
| Counting reads per gene/feature | -> rna-quantification/featurecounts-counting | the seqid/strand landmines live there; set `-s` from chemistry |
| Judge whether the annotation itself is sound | -> genome-annotation/annotation-qc | this skill operates on the file, not its quality |

## Walk the Gene Tree and Derive Introns (gffutils)

**Goal:** Traverse gene -> transcript -> exon and reconstruct features (introns) that the file does not store explicitly.

**Approach:** Build a SQLite DB once (disabling gene/transcript inference when those lines already exist, for a ~100x speedup), then query children ordered by position and synthesize introns from the exon gaps.

```python
import gffutils

# disable_infer_* is GTF-only and applies when gene/transcript lines ALREADY exist (modern GENCODE/Ensembl) -> ~100x faster
db = gffutils.create_db('annotation.gtf', 'annotation.db', force=True,
                        disable_infer_genes=True, disable_infer_transcripts=True,
                        merge_strategy='create_unique')

gene = db['ENSG00000141510']                                    # gffutils returns 1-based coords (raw record)
for tx in db.children(gene, featuretype=['mRNA', 'transcript'], order_by='start'):
    exons = list(db.children(tx, featuretype='exon', order_by='start'))
    introns = list(db.interfeatures(exons, new_featuretype='intron'))   # introns are not stored - derived from exon gaps
    print(tx.id, len(exons), 'exons', len(introns), 'introns')
```

A modern GTF without `disable_infer_*` triggers the slow inference/merge machinery; an *older* minimal GTF lacking gene/transcript lines needs inference ON so gffutils reconstructs the envelopes. Match the flag to the file. introns, UTRs (`exon - CDS`), and TSS are derived, not stored - never infer biological absence from a missing feature line.

## Convert Formats and Extract Sequences (gffread)

gffread is genome-aware and respects the stop-codon convention, so it is