Applied - GFF and Genomic Intervals

In one line: parsing annotation files and answering "what overlaps what" — where pandas runs out and interval libraries take over.

Uses: Reading and Writing Data · String Accessor · Coordinate Systems · Merging and Joining · SeqFeature and Location

Reading GFF/GTF

Nine tab-separated columns, no header, # comments, 1-based closed:

import pandas as pd

GFF_COLS = ["seqid", "source", "type", "start", "end",
            "score", "strand", "phase", "attributes"]

def read_gff(path):
    df = pd.read_csv(
        path, sep="\t", comment="#", names=GFF_COLS,
        dtype={"seqid": "category", "source": "category", "type": "category",
               "start": "int32", "end": "int32", "strand": "category"},
        na_values=["."],
    )
    df["start0"] = df["start"] - 1        # ← convert to 0-based half-open ONCE
    df["end0"] = df["end"]
    df["length"] = df["end"] - df["start"] + 1
    return df

gff = read_gff("gencode.v45.gtf.gz")
gff["type"].value_counts()

The start0/end0 naming is the discipline. A bare start in a frame that has touched both BED and GFF data is a bug waiting to happen. See Coordinate Systems.

Parsing the attributes column

GFF3 and GTF use different syntax for column 9:

# GFF3:  ID=gene-TP53;Name=TP53;gene_biotype=protein_coding
gff["gene_id"] = gff["attributes"].str.extract(r"(?:^|;)ID=([^;]+)")
gff["gene_name"] = gff["attributes"].str.extract(r"(?:^|;)Name=([^;]+)")
gff["parent"] = gff["attributes"].str.extract(r"(?:^|;)Parent=([^;]+)")

# GTF:  gene_id "ENSG00000141510"; gene_name "TP53";
gff["gene_id"] = gff["attributes"].str.extract(r'gene_id "([^"]+)"')
gff["gene_name"] = gff["attributes"].str.extract(r'gene_name "([^"]+)"')
gff["transcript_id"] = gff["attributes"].str.extract(r'transcript_id "([^"]+)"')

# strip Ensembl version suffixes so joins will match
gff["gene_id_nov"] = gff["gene_id"].str.split(".").str[0]

Always check how many failed to extract:

print(gff["gene_id"].isna().sum(), "rows without a gene_id")

For full GFF3 parsing with parent–child relationships, use gffutils — it builds a SQLite database and lets you query db.children(gene, featuretype="exon") properly. Regex extraction is fine for flat queries and wrong for hierarchy traversal.

Deriving structures

genes = gff.query("type == 'gene'")
exons = gff.query("type == 'exon'")

# exons per transcript
exons.groupby("transcript_id").size().describe()

# transcript span from its exons
tx = exons.groupby("transcript_id").agg(
    seqid=("seqid", "first"), strand=("strand", "first"),
    start0=("start0", "min"), end0=("end0", "max"),
    n_exons=("start0", "size"),
    exonic_bp=("length", "sum"),
)

# introns, derived from consecutive exons
def introns_of(g):
    g = g.sort_values("start0")
    return pd.DataFrame({
        "seqid": g["seqid"].iloc[0],
        "start0": g["end0"].values[:-1],
        "end0": g["start0"].values[1:],
    })

introns = exons.groupby("transcript_id", group_keys=True).apply(
    introns_of, include_groups=False
).reset_index()

Interval overlap — the core problem

pandas cannot do this efficiently. A merge matches on equality; overlap is a range predicate. A naive cross-join of 1 million variants against 60,000 genes is 60 billion comparisons.

The right tools

import pyranges as pr

variants = pr.PyRanges(chromosomes=vcf["CHROM"],
                       starts=vcf["POS"] - 1, ends=vcf["POS"])
genes_pr = pr.PyRanges(chromosomes=genes["seqid"],
                       starts=genes["start0"], ends=genes["end0"])

annotated = variants.join(genes_pr)          # interval join
variants.overlap(genes_pr)                    # boolean filter
variants.nearest(genes_pr)                    # nearest feature + distance
genes_pr.merge()                              # collapse overlapping intervals
a.subtract(b); a.intersect(b)

pyranges wraps nested containment lists — O(log n) per query. bioframe offers a similar API that stays closer to pandas DataFrames. pybedtools wraps BEDTools, which is the gold standard if you already have it installed.

Pure Python, for smaller problems

from intervaltree import IntervalTree
from collections import defaultdict

trees = defaultdict(IntervalTree)
for r in genes.itertuples():
    trees[r.seqid][r.start0:r.end0] = r.gene_name       # half-open

def genes_at(chrom, pos0):
    return [iv.data for iv in trees[chrom][pos0]]

vcf["genes"] = [genes_at(c, p - 1) for c, p in zip(vcf["CHROM"], vcf["POS"])]
vcf = vcf.explode("genes")

Build once, query many times. Fine for tens of thousands of intervals.

Sorted-array approach

For a one-off overlap of two sorted position lists, np.searchsorted is fast and dependency-free:

import numpy as np

starts = genes_chr["start0"].to_numpy()
ends = genes_chr["end0"].to_numpy()
order = np.argsort(starts)
starts, ends = starts[order], ends[order]

idx = np.searchsorted(starts, positions, side="right") - 1
hit = (idx >= 0) & (positions < ends[idx])

Correct only for non-overlapping intervals — which genes are not. Useful for exons within a transcript, or for binned data.

The classic mistakes

Off-by-one at the format boundary:

# GFF is 1-based closed; a Python slice is 0-based half-open
seq[gff_start - 1 : gff_end]           # ✓
seq[gff_start : gff_end]                # ✗ — drops the first base

The overlap predicate. Two half-open intervals overlap when:

a_start < b_end and b_start < a_end

Note strict < on both. Using <= makes book-ended intervals (where one ends exactly where the next begins) count as overlapping, which is wrong for half-open coordinates and is exactly how adjacent exons get merged.

Strand. Overlap is usually strand-agnostic, but "the nearest upstream gene" is not — upstream means lower coordinates on the plus strand and higher on the minus strand.

Practical recipes

# gene body coverage: how much of each gene is covered by an interval set?
cov = pr.PyRanges(cov_df).merge()
genes_pr.coverage(cov)

# promoters (2 kb upstream, strand-aware)
genes["prom_start0"] = np.where(genes["strand"] == "+",
                                genes["start0"] - 2000, genes["end0"])
genes["prom_end0"] = np.where(genes["strand"] == "+",
                              genes["start0"], genes["end0"] + 2000)
genes["prom_start0"] = genes["prom_start0"].clip(lower=0)

# fixed-width genomic bins
BIN = 1_000_000
variants["bin"] = variants["POS"] // BIN
density = variants.groupby(["CHROM", "bin"], observed=True).size()

# is a variant in an exon, and if so which?
exon_tree = build_trees(exons)
vcf["in_exon"] = [bool(exon_tree[c][p - 1]) for c, p in zip(vcf.CHROM, vcf.POS)]

# distance to the nearest splice site — feeds splice-effect prediction
boundaries = np.sort(np.concatenate([exons["start0"], exons["end0"]]))
i = np.searchsorted(boundaries, positions)
dist = np.minimum(np.abs(positions - boundaries[np.clip(i-1, 0, None)]),
                  np.abs(positions - boundaries[np.clip(i, None, len(boundaries)-1)]))

The promoter definition is worth noting as an example of strand-awareness done explicitly — the naive start - 2000 is wrong for half your genes.

Chromosome naming, again

set(gff["seqid"]) & set(vcf["CHROM"])     # ← if this is empty, that is your bug

GENCODE uses chr1, Ensembl uses 1. Normalise at load, assert the intersection is non-empty. This check has saved more debugging time than any other single line in this vault.

Common mistakes

• Off-by-one converting GFF to 0-based.
• <= in the overlap predicate, merging book-ended intervals.
• A pandas cross-join for overlap.
• chr1 vs 1.
• Ignoring strand for directional queries.
• Regex-parsing GFF3 hierarchy instead of using gffutils.
• Version suffixes on gene IDs breaking joins.
• Not clipping negative coordinates after extending intervals.
• Assuming one transcript per gene.

See also

Coordinate Systems · Applied - Genomic File Formats · SeqFeature and Location · String Accessor · Merging and Joining · Applied - Annotating a Variant Table