Applied - GFF and Genomic Intervals

Core Principle: Genomic feature annotation parsing resolves 1-based closed coordinates in GFF/GTF files, unpacks structured attribute columns, and computes interval intersections against reference genes.

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

Exercise

Answer "what overlaps what" without an interval library. The two mistakes the lesson warns about are both live here: the coordinate systems differ between the two inputs, and overlapping features double-count bases unless you merge them first.

Write overlaps(text), returning a list of output lines. The first input line is the GFF3 feature type to consider. What follows is a GFF3 annotation, then a line containing exactly --, then the query intervals as chrom, start, end in 0-based half-open coordinates, tab separated. Ignore blank lines and lines beginning with #. GFF3 is 1-based closed, so a feature overlaps a query when start - 1 < query_end and end > query_start. Normalise chromosome names on both sides by prefixing chr where it is missing, so that 17 and chr17 are the same sequence.

Print features=<n> queries=<n>, where the feature count is the number of annotation records of the requested type, whatever chromosome they lie on. Then print one line per query in input order, holding chrom:start-end, the number of overlapping features, the number of query bases covered by at least one of them, and the overlapping features' ID attributes sorted ascending and joined with commas, tab separated. Compute the covered bases by clipping each overlapping feature to the query and merging the clipped intervals, so a base under two exons counts once. Use . for the name field when nothing overlaps. If there are no query intervals, print no queries after the first line.