Gene models

A GeneModels object holds a gene annotation — genes, their transcripts, and the exons of each transcript — parsed out of whatever format the resource stores it in (GTF, refFlat, refSeq, CCDS and several others; Genomic resources and repositories lists the formats and the format key that selects one).

Unlike the other resource types, gene models are loaded wholly into memory rather than read from an open file. The consequence is that load() is expensive and close() does nothing:

from gain.genomic_resources.repository_factory import build_genomic_resource_repository
from gain.genomic_resources.gene_models import build_gene_models_from_resource_id

grr = build_genomic_resource_repository()
genes = build_gene_models_from_resource_id("hg38/gene_models/MANE/1.5", grr).load()

is_loaded() distinguishes a built-but-empty object from a loaded one, which library code that accepts either should check rather than calling load a second time.

Two lookups

gene_models_by_gene_name() goes from a gene symbol to its transcripts, and returns None when the symbol is unknown — not an empty list, so a missing gene and a gene with no transcripts stay distinguishable.

gene_models_by_location() goes the other way, from a position or an interval to the transcripts that overlap it, and returns a list:

transcripts = genes.gene_models_by_gene_name("TP53")
for tm in transcripts or []:
    print(tm.tr_id, tm.chrom, tm.tx, tm.is_coding())

overlapping = genes.gene_models_by_location("chr17", 7_676_000, 7_690_000)

join_gene_models() merges several loaded objects into one — the way to query a primary annotation and a supplementary one as a single set.

Transcripts and exons

A TranscriptModel is the unit most analyses work with. Its coordinates are plain attributes (chrom, strand, tx for the transcribed interval, cds for the coding one, exons), and its methods derive the parts you usually want:

tm = transcripts[0]
if tm.is_coding():
    for region in tm.cds_regions():
        print(region.chrom, region.start, region.stop)
    print("5' UTR:", tm.utr5_regions())
    print("coding length:", tm.cds_len())

cds_regions(), utr5_regions(), utr3_regions() and all_regions() each return a list of BedRegion. Note that utr5/utr3 are named for the transcript’s orientation, so which end of the interval they fall on depends on strand.

The exon list is Exon objects. An exon’s frame is the reading frame the effect-annotation engine needs, and it is not populated by parsing alone: None means “not computed yet”, not “non-coding” (a non-coding exon has -1). A resource whose source format carries no frames holds None until update_frames() fills them in. Exon’s three coordinates are described in its class documentation below rather than as separate entries.

API

gain.genomic_resources.gene_models.build_gene_models_from_resource_id(resource_id: str, grr: GenomicResourceRepo | None = None) → GeneModels

Load gene models from a genomic resource id.

class gain.genomic_resources.gene_models.GeneModels(resource: GenomicResource)

Manage and query gene model data from genomic resources.

This class provides access to gene models loaded from various file formats (GTF, refFlat, refSeq, CCDS, etc.) and offers efficient querying by gene name or genomic location.

The class maintains three internal data structures: - transcript_models: Dict mapping transcript IDs to TranscriptModel objects - gene_models: Dict mapping gene names to lists of TranscriptModel objects - _tx_index: IntervalTree index for fast location-based queries

Variables:

• resource (GenomicResource) – The genomic resource containing gene models.

• config (dict) – Validated configuration from the resource.

• reference_genome_id (str | None) – ID of the reference genome.

• gene_models (dict[str, list[TranscriptModel]]) – Gene name to transcript models mapping.

• transcript_models (dict[str, TranscriptModel]) – Transcript ID to transcript model mapping.

Example

>>> from gain.genomic_resources.gene_models         ...     .gene_models_factory import build_gene_models_from_file
>>> gene_models = build_gene_models_from_file("genes.gtf")
>>> gene_models.load()
>>> # Query by gene name
>>> tp53_transcripts = gene_models.gene_models_by_gene_name("TP53")
>>> # Query by location
>>> transcripts = gene_models.gene_models_by_location("chr17", 7676592)

Note

The gene models must be loaded using the load() method before queries can be performed. The class is thread-safe for concurrent access.

close() → None

Do nothing: gene models hold no open file and stay loaded.

Here so a gene models resource closes like every other.

gene_models_by_gene_name(name: str) → list[TranscriptModel] | None

Retrieve all transcript models for a specific gene.

Parameters:

name (str) – The gene name/symbol to search for.

Returns:

List of transcript models for the

gene, or None if the gene is not found.

Return type:

list[TranscriptModel] | None

Example

>>> transcripts = gene_models.gene_models_by_gene_name("BRCA1")
>>> if transcripts:
...     print(f"BRCA1 has {len(transcripts)} transcript variants")

gene_models_by_location(chrom: str, pos_begin: int, pos_end: int | None = None) → list[TranscriptModel]

Retrieve transcripts overlapping a genomic position or region.

This method uses an interval tree index for efficient querying of transcripts by genomic coordinates.

Parameters:

• chrom (str) – The chromosome name (e.g., “chr1”, “17”).

• pos_begin (int) – The start position (1-based, inclusive).

• pos_end (int | None) – The end position (1-based, inclusive). If None, queries a single position.

Returns:

List of TranscriptModel objects whose transcript regions overlap the query position/region. Returns empty list if no overlaps found.

Return type:

list[TranscriptModel]

Example

>>> # Query single position
>>> models = gene_models.gene_models_by_location("chr17", 7676592)
>>> # Query region
>>> models = gene_models.gene_models_by_location(
...     "chr17", 7661779, 7687550
... )
>>> for tm in models:
...     print(f"{tm.gene}: {tm.tr_id}")

Note

Positions are swapped automatically if pos_end < pos_begin.

gene_names() → list[str]

Get list of all gene names in the loaded gene models.

Returns:

List of gene names (symbols).

Return type:

list[str]

Example

>>> gene_models.load()
>>> genes = gene_models.gene_names()
>>> print(f"Loaded {len(genes)} genes")

static get_schema() → dict[str, Any]

The schema a gene_models resource’s config is checked against.

The base resource schema plus filename, format, gene_mapping and chrom_mapping.

has_chromosome(chrom: str) → bool

Check if a chromosome has any gene models.

Parameters:

chrom (str) – The chromosome name to check.

Returns:

True if the chromosome has gene models, False otherwise.

Return type:

bool

Example

>>> if gene_models.has_chromosome("chr1"):
...     print("Chromosome 1 has gene annotations")

is_loaded() → bool

Check if gene models have been loaded.

Returns:

True if load() has been called and completed, False otherwise.

Return type:

bool

Example

>>> if not gene_models.is_loaded():
...     gene_models.load()

static join_gene_models(*gene_models: GeneModels) → GeneModels

Merge multiple gene models into a single GeneModels object.

This combines transcript models from multiple sources into one unified gene models object.

Parameters:

• *gene_models (GeneModels) – Two or more GeneModels objects to

• merge.

Returns:

New GeneModels object containing all transcripts.

Return type:

GeneModels

Raises:

ValueError – If fewer than 2 gene models provided.

Example

>>> gm1 = build_gene_models_from_file("genes1.gtf")
>>> gm2 = build_gene_models_from_file("genes2.gtf")
>>> merged = GeneModels.join_gene_models(gm1, gm2)

Note

Transcript IDs should be unique across all input gene models.

load() → GeneModels

Load gene models from the genomic resource.

This method parses the gene model file and builds internal indexes for efficient querying. It is thread-safe and will only load once.

Returns:

Self, for method chaining.

Return type:

GeneModels

Example

>>> gene_models = build_gene_models_from_file("genes.gtf")
>>> gene_models.load()
>>> num_transcripts = len(gene_models.transcript_models)
>>> print(f"Loaded {num_transcripts} transcripts")

Note

Calling load() multiple times is safe - subsequent calls return immediately if already loaded.

property resource_id: str

The id of the gene models resource this object wraps.

class gain.genomic_resources.gene_models.TranscriptModel(gene: str, tr_id: str, tr_name: str, chrom: str, strand: str, *, tx: tuple[int, int], cds: tuple[int, int], exons: list[Exon] | None = None, attributes: dict[str, Any] | None = None)

Represent a transcript with all its structural features.

A transcript model contains complete information about a gene transcript, including its genomic location, exon structure, coding regions, and additional attributes from the source annotation.

Variables:

• gene (str) – Gene name/symbol (e.g., “TP53”).

• tr_id (str) – Transcript identifier, unique within the gene models.

• tr_name (str) – Original transcript name from source annotation.

• chrom (str) – Chromosome name (e.g., “chr17”, “17”).

• strand (str) – Strand orientation (“+” or “-“).

• tx (tuple[int, int]) – Transcript start and end positions (1-based, closed interval).

• cds (tuple[int, int]) – Coding sequence start and end positions (1-based, closed interval). For non-coding transcripts, cds[0] >= cds[1].

• exons (list[Exon]) – List of Exon objects in genomic order.

• attributes (dict[str, Any]) – Additional annotation attributes (e.g., gene_biotype, gene_version).

Example

>>> from gain.genomic_resources.gene_models.transcript_models import \
...     TranscriptModel, Exon
>>> tm = TranscriptModel(
...     gene="TP53",
...     tr_id="ENST00000269305",
...     tr_name="TP53-201",
...     chrom="17",
...     strand="-",
...     tx=(7661779, 7687550),
...     cds=(7668402, 7687490),
...     exons=[Exon(7661779, 7661822), Exon(7668402, 7669690)],
...     attributes={"gene_biotype": "protein_coding"},
... )
>>> print(f"Coding: {tm.is_coding()}")
Coding: True
>>> regions = tm.cds_regions()
>>> print(f"CDS has {len(regions)} regions")

Note

• All coordinates use 1-based, closed intervals

• CDS includes both start and stop codons

• Exons should be in genomic order (not necessarily 5’ to 3’)

all_regions(ss_extend: int = 0, prom: int = 0) → list[BedRegion]

Get all transcript regions with optional extensions.

Returns all exonic regions, optionally extending into splice sites and promoter regions.

Parameters:

• ss_extend (int) – Number of bases to extend into splice sites at coding exon boundaries. Default is 0.

• prom (int) – Number of bases to extend into promoter region upstream of transcription start. Default is 0.

Returns:

List of all transcript regions, potentially

extended based on parameters.

Return type:

list[BedRegion]

Example

>>> # Basic exonic regions
>>> regions = transcript.all_regions()
>>> # With splice site extension
>>> regions = transcript.all_regions(ss_extend=3)
>>> # With promoter region
>>> regions = transcript.all_regions(prom=2000)

Note

Promoter extension is strand-aware: extends upstream of the transcription start (before first exon for +, after last for -).

calc_frames() → list[int]

Calculate reading frame for each exon.

Computes the codon reading frame (0, 1, or 2) for each exon based on the CDS coordinates and strand orientation.

Frames are computed from the CDS intersected with the exon set, so bases annotated outside the exons – an incomplete stop codon, say – shift no frame.

Returns:

Reading frame for each exon. Values are:

• 0, 1, or 2 for coding exons (bases into current codon)

• -1 for non-coding exons or non-coding transcripts

Return type:

list[int]

Example

>>> frames = transcript.calc_frames()
>>> for exon, frame in zip(transcript.exons, frames):
...     if frame >= 0:
...         print(f"Exon {exon.start}-{exon.stop}: frame {frame}")

Note

Frame calculation is strand-aware and considers exon order. Use update_frames() to set frame attribute on Exon objects.

cds_len() → int

The coding length of the transcript, in base pairs.

The sum of cds_regions(), so it counts coding sequence only and is 0 for a non-coding transcript.

cds_regions(ss_extend: int = 0) → list[BedRegion]

Compute coding sequence (CDS) regions.

Extracts the portions of exons that contain coding sequence, optionally extending into splice sites.

Parameters:

ss_extend (int) – Number of bases to extend into splice sites at exon boundaries. Default is 0 (no extension).

The CDS is intersected with the exon set first, so the result never leaves an exon. An empty list is returned for a non-coding transcript and for one whose CDS misses its exons entirely.

Returns:

List of BedRegion objects representing CDS

segments. Returns empty list for non-coding transcripts.

Return type:

list[BedRegion]

Example

>>> cds_regions = transcript.cds_regions()
>>> for region in cds_regions:
...     print(f"{region.chrom}:{region.start}-{region.stop}")
>>> # With splice site extension
>>> extended = transcript.cds_regions(ss_extend=3)

Note

CDS regions include both start and stop codons where the source annotates them – an incomplete transcript may have neither. Use collect_gtf_cds_regions() from serialization module to exclude the stop codon for GTF format.

is_coding() → bool

Check if this transcript is protein-coding.

The CDS is measured against the exon set, so a transcript whose annotated CDS misses its exons entirely translates nothing and is not coding. This keeps is_coding() in step with cds_regions(), which callers index into without checking.

Returns:

True if the transcript has a coding region (CDS),

False for non-coding transcripts.

Return type:

bool

Example

>>> if transcript.is_coding():
...     cds_len = transcript.cds_len()
...     print(f"CDS length: {cds_len}bp")

test_frames() → bool

Verify that exon frames are correctly set.

Compares the frame attribute of each exon with the calculated frame to ensure consistency.

Returns:

True if all exon frames match calculated values,

False otherwise.

Return type:

bool

Example

>>> transcript.update_frames()
>>> assert transcript.test_frames()

total_len() → int

The spliced length of the transcript, in base pairs.

The sum of the exon lengths – NOT the genomic span tx, which includes the introns between them.

update_frames() → None

Update the frame attribute of all exons.

Calculates reading frames using calc_frames() and updates the frame attribute of each Exon object.

Example

>>> transcript.update_frames()
>>> for exon in transcript.exons:
...     print(f"Exon frame: {exon.frame}")

Note

This modifies the Exon objects in place.

utr3_len() → int

The length of the 3’ untranslated region, in base pairs.

The sum of utr3_regions(). “3’” is relative to the transcript’s own orientation, so on a - strand transcript this is the lower-coordinate end.

utr3_regions() → list[BedRegion]

Get 3’ untranslated region (3’ UTR) segments.

The 3’ UTR extends from the stop codon (translation end) to the transcription end. Strand orientation is considered.

The CDS is intersected with the exon set first, so a CDS overrunning its terminal exon splits the transcript exactly as a flush one would.

Returns:

List of 3’ UTR regions. Returns empty list for

non-coding transcripts.

Return type:

list[BedRegion]

Example

>>> utr3 = transcript.utr3_regions()
>>> utr3_length = sum(r.stop - r.start + 1 for r in utr3)
>>> print(f"3' UTR: {utr3_length}bp")

Note

For positive strand: regions after CDS end. For negative strand: regions before CDS start.

utr5_len() → int

The length of the 5’ untranslated region, in base pairs.

The sum of utr5_regions(). “5’” is relative to the transcript’s own orientation, so on a - strand transcript this is the higher-coordinate end.

utr5_regions() → list[BedRegion]

Get 5’ untranslated region (5’ UTR) segments.

The 5’ UTR extends from the transcription start to the start codon (translation start). Strand orientation is considered.

The CDS is intersected with the exon set first, so a CDS overrunning its terminal exon splits the transcript exactly as a flush one would.

Returns:

List of 5’ UTR regions. Returns empty list for

non-coding transcripts.

Return type:

list[BedRegion]

Example

>>> utr5 = transcript.utr5_regions()
>>> utr5_length = sum(r.stop - r.start + 1 for r in utr5)
>>> print(f"5' UTR: {utr5_length}bp")

Note

For positive strand: regions before CDS start. For negative strand: regions after CDS end.

class gain.genomic_resources.gene_models.Exon(start: int, stop: int, frame: int | None = None)

Represent a single exon within a transcript.

An exon is a segment of a transcript that is retained in the mature RNA after splicing. This class stores the genomic coordinates and codon reading frame of an exon.

Variables:

• start (int) – Genomic start position (1-based, inclusive).

• stop (int) – Genomic end position (1-based, inclusive).

• frame (int | None) – Codon reading frame (0, 1, or 2) for coding exons, and -1 for a non-coding one – see calc_frames. None means the frame has not been computed yet, not that the exon is non-coding: it is what an exon built without one holds until update_frames fills it in. Serializing a model still holding None is refused rather than guessed.

Example

>>> exon = Exon(start=100, stop=200, frame=0)
>>> print(exon.start, exon.stop)
100 200
>>> exon.contains((150, 160))
True

contains(region: tuple[int, int]) → bool

Check if this exon fully contains a genomic region.

Parameters:

region (tuple[int, int]) – A (start, stop) position tuple to check.

Returns:

True if the region is fully contained within this exon.

Return type:

bool

Example

>>> exon = Exon(100, 200)
>>> exon.contains((150, 160))
True
>>> exon.contains((50, 250))
False