Gene models

Gene models resources describe genes, their transcripts, and the exons that make up each transcript. GAIn provides the GeneModels Python object for accessing these annotations. It reads the resource’s data and represents individual transcripts as TranscriptModel objects, making it possible to find transcripts by gene name or genomic location and examine their coding regions, untranslated regions, and exons.

Gene models can be stored in several formats, including GTF, refFlat, refSeq, and CCDS. GAIn interprets these formats using the format setting in the resource definition, as described in Genomic resources and repositories.

Accessing gene models

A GeneModels object can be built from a resource ID and a repository using build_gene_models_from_resource_id():

import gain.grr as grr

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

If you already have a GenomicResource object, you can use build_gene_models_from_resource() instead:

resource = repo.get_resource("hg38/gene_models/MANE/1.5")
genes = grr.build_gene_models_from_resource(resource)

Both builder functions are available through gain.grr. Calls that resolve to the same resource return the same GeneModels object, which may already have been loaded by another part of the program.

Loading gene models

Unlike a reference genome, which retrieves sequences from an open FASTA file, a GeneModels object loads its annotation data into memory. Call load() before querying genes or transcripts:

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

The load() method returns the GeneModels object itself, allowing it to be chained to the builder as shown above. Use is_loaded() to check whether an object has already been loaded, particularly in code that may receive either a newly built or an already-loaded object.

Because the annotation is kept in memory rather than accessed through an open file handle, close() does nothing. Gene models do not require the open() and close() operations used for some other resource types.

Finding transcripts

Once gene models are loaded, you can look up transcripts by gene name or genomic location.

gene_models_by_gene_name() returns the transcripts associated with a gene symbol. If the gene name is not found, it returns None rather than an empty list. This distinguishes an unknown gene from one that has no transcripts.

gene_models_by_location() returns a list of transcripts overlapping a specified genomic position or interval.

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
)

The first lookup finds transcripts of TP53 and prints each transcript’s identifier, chromosome, transcribed interval, and coding status. The expression transcripts or [] allows the loop to run safely even if the gene name is not found. The second lookup finds transcripts overlapping the specified interval on chromosome 17.

Working with transcripts and exons

A TranscriptModel represents an individual transcript. Its attributes include chrom (chromosome), strand (transcription orientation), tx (transcribed interval), cds (coding interval), and exons (the transcript’s exons). Its methods can identify the coding and untranslated portions of a transcript and calculate their lengths.

For example, after finding the transcripts of a gene, you can inspect the coding regions of one transcript:

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())

The methods cds_regions(), utr5_regions(), utr3_regions(), and all_regions() each return a list of BedRegion objects. The 5’ and 3’ untranslated regions are defined relative to the direction of transcription, so their positions along the chromosome depend on the transcript’s strand.

The exons attribute contains Exon objects. Each exon has a frame value used by the effect-annotation machinery to determine its coding reading frame. An exon with frame = None has not had its frame computed, while frame = -1 identifies a non-coding exon. These values should not be confused. If the source annotation does not provide reading frames, they remain None until update_frames() calculates them. The Exon class documentation in the API reference describes its three attributes, start, stop and frame.

Combining and constructing gene models

join_gene_models() combines several loaded GeneModels objects into one. For example, a primary gene annotation and a supplementary annotation can be merged so that their transcripts can be queried together. The resulting object is already loaded and can be used wherever a loaded GeneModels object is expected, including by effect annotators. Each input must be loaded before it can be joined. An unloaded object is rejected rather than treated as an empty annotation.

from_transcript_models() constructs a loaded, indexed GeneModels object from a resource and a dictionary of transcript models (dict[str, TranscriptModel]). This is useful when transcript annotations have been assembled or modified programmatically and need to be used through the same interface as gene models loaded from a GRR.

API

The classes below are defined in the module gain.genomic_resources.gene_models, and so are the builders build_gene_models_from_resource_id and build_gene_models_from_resource. GeneModels is also exported by gain.grr, as grr.GeneModels. The builders are also exported by gain.grr, as grr.build_gene_models_from_resource_id and grr.build_gene_models_from_resource.

`build_gene_models_from_resource`(resource)	Build the gene models of a genomic resource.
`build_gene_models_from_resource_id`(resource_id)	Build the gene models of the resource with a given id.
`GeneModels`(resource)	Manage and query gene model data from genomic resources.
`GeneModels.close`()	Do nothing: gene models hold no open file and stay loaded.
`GeneModels.from_transcript_models`(resource, ...)	Build a loaded, indexed gene models object over `resource`.
`GeneModels.gene_models_by_gene_name`(name)	Retrieve all transcript models for a specific gene.
`GeneModels.gene_models_by_location`(chrom, ...)	Retrieve transcripts overlapping a genomic position or region.
`GeneModels.gene_names`()	Get list of all gene names in the loaded gene models.
`GeneModels.get_schema`()	The schema a `gene_models` resource's config is checked against.
`GeneModels.has_chromosome`(chrom)	Check if a chromosome has any gene models.
`GeneModels.is_loaded`()	Check whether this object holds usable gene models.
`GeneModels.join_gene_models`(*gene_models)	Merge loaded gene models into a single, loaded GeneModels object.
`GeneModels.load`()	Load gene models from the genomic resource.
`GeneModels.resource_id`	The id of the gene models resource this object wraps.
`TranscriptModel`(gene, tr_id, tr_name, chrom, ...)	Represent a transcript with all its structural features.
`TranscriptModel.all_regions`([ss_extend, prom])	Get all transcript regions with optional extensions.
`TranscriptModel.calc_frames`()	Calculate reading frame for each exon.
`TranscriptModel.cds_len`()	The coding length of the transcript, in base pairs.
`TranscriptModel.cds_regions`([ss_extend])	Compute coding sequence (CDS) regions.
`TranscriptModel.is_coding`()	Check if this transcript is protein-coding.
`TranscriptModel.test_frames`()	Verify that exon frames are correctly set.
`TranscriptModel.total_len`()	The spliced length of the transcript, in base pairs.
`TranscriptModel.update_frames`()	Update the frame attribute of all exons.
`TranscriptModel.utr3_len`()	The length of the 3' untranslated region, in base pairs.
`TranscriptModel.utr3_regions`()	Get 3' untranslated region (3' UTR) segments.
`TranscriptModel.utr5_len`()	The length of the 5' untranslated region, in base pairs.
`TranscriptModel.utr5_regions`()	Get 5' untranslated region (5' UTR) segments.
`Exon`(start, stop[, frame])	Represent a single exon within a transcript.
`Exon.contains`(region)	Check if this exon fully contains a genomic region.

build_gene_models_from_resource(resource: GenomicResource | None) → GeneModels[source]

Build the gene models of a genomic resource.

The returned object may not be loaded. Calls with the same resource return the same object, so an earlier caller may have loaded it.

Parameters:: resource (GenomicResource | None) – A resource of type gene_models.
Returns:: The gene models of the resource. Check is_loaded() or call load() before any query.
Return type:: GeneModels
Raises:: ValueError – If resource is None or is not of type gene_models.

build_gene_models_from_resource_id(resource_id: str, grr: GenomicResourceRepo | None = None) → GeneModels[source]

Build the gene models of the resource with a given id.

The returned object may not be loaded. Calls that find the same resource return the same object, so an earlier caller may have loaded it.

Parameters:

resource_id (str) – The id of a resource of type gene_models in the repository.
grr (GenomicResourceRepo | None) – The repository that resolves resource_id. When it is None, the function builds the default repository with build_genomic_resource_repository().

Returns:

The gene models of the resource. Check is_loaded() or call load() before any query.

Return type:

GeneModels

Raises:

FileNotFoundError – If a protocol-backed repository has no resource with this id.
ValueError – If a group repository has no resource with this id, or if the resource is not of type gene_models.

class GeneModels(resource: GenomicResource)[source]

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 (see is_loaded()) before queries can be performed. The class is thread-safe for concurrent access.

close() → None[source]

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

Here so a gene models resource closes like every other.

classmethod from_transcript_models(resource: GenomicResource, transcript_models: dict[str, TranscriptModel]) → GeneModels[source]

Build a loaded, indexed gene models object over resource.

The transcripts are taken as given, in the order given (the serializers write them in that order); the resource’s chrom mapping is not applied.

Parameters:

resource (GenomicResource) – The gene models resource the result answers to (resource_id, reference genome).
transcript_models (dict[str, TranscriptModel]) – Transcript ID to transcript model.

Returns:

A loaded model over resource.

Return type:

GeneModels

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

Retrieve all transcript models for a specific gene.

Parameters:

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

Returns:

list[TranscriptModel] | None – The transcript models of the gene, or None if the gene models hold no gene with this name. The result is never an empty list.
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][source]

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 (for example “chr1” or “17”).
pos_begin (int) – The start position (1-based, inclusive).
pos_end (int | None, default None) – The end position (1-based, inclusive). If None, the query is the single position pos_begin.

Returns:

list[TranscriptModel] – The transcript models whose transcript region overlaps the query position or region. The list is empty if no transcript overlaps or if the chromosome has no gene models.
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][source]

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

Returns:

list[str] – The gene names (symbols), one for each gene. The list is empty if the gene models hold no genes, for example when they are not loaded yet.
Example – >>> gene_models.load() >>> genes = gene_models.gene_names() >>> print(f”Loaded {len(genes)} genes”)

static get_schema() → dict[str, Any][source]

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

The base resource schema plus filename, format, gene_mapping and chrom_mapping. filename is required: the loader reads it unconditionally.

Returns:: The Cerberus schema of the resource config.
Return type:: dict[str, Any]

has_chromosome(chrom: str) → bool[source]

Check if a chromosome has any gene models.

Parameters:

chrom (str) – The chromosome name to check.

Returns:

bool – True if at least one transcript model is on the chromosome, False otherwise.
Example – >>> if gene_models.has_chromosome(“chr1”): … print(“Chromosome 1 has gene annotations”)

is_loaded() → bool[source]

Check whether this object holds usable gene models.

True once the model holds indexed transcripts, however it got them: load() parsed its resource, or it was built loaded by join_gene_models or from_transcript_models. load() on a loaded model is a no-op.

Returns:

bool – True if the models are loaded and indexed, False if nobody loaded them yet.
Example – >>> if not gene_models.is_loaded(): … gene_models.load()

static join_gene_models(*gene_models: GeneModels) → GeneModels[source]

Merge loaded gene models into a single, loaded GeneModels object.

The result holds every input’s transcripts, answers to the first input’s resource, and is loaded (see is_loaded()).

Parameters:

*gene_models (GeneModels) – Two or more loaded GeneModels objects to merge.

Returns:

A new loaded GeneModels object that holds all transcripts.

Return type:

GeneModels

Raises:

ValueError – If fewer than 2 gene models are provided, or if any of them was never loaded: merging its empty transcripts would silently drop everything it was meant to bring.

Example: –

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

Note: – Transcript IDs should be unique across all input gene models; a later input’s transcript replaces an earlier one’s under the same ID.

load() → GeneModels[source]

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:

GeneModels – This object, for method chaining.
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 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)[source]

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][source]

Get all transcript regions with optional extensions.

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

Parameters:

ss_extend (int, default 0) – The number of bases to extend into the splice sites at the coding exon boundaries. 0 means no extension.
prom (int, default 0) – The number of bases to extend into the promoter region upstream of the transcription start. 0 means no extension.

Returns:

list[BedRegion] – One region for each exon, extended as the parameters ask.
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][source]

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:

The reading frame of each exon, in the order of exons. A value is 0, 1 or 2 for a coding exon (the bases into the current codon) and -1 for a non-coding exon or a non-coding transcript.

Return type:

list[int]

Raises:

ValueError – If the number of computed frames differs from the number of exons.

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[source]

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.

Returns:: The coding length in base pairs.
Return type:: int

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

Compute coding sequence (CDS) regions.

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

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.

Parameters:

ss_extend (int, default 0) – The number of bases to extend into the splice sites at the exon boundaries. 0 means no extension.

Returns:

list[BedRegion] – The CDS segments, one for each exon the CDS touches. The list is empty for a non-coding transcript.
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[source]

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:

bool – True if the transcript has a coding region (CDS), False for non-coding transcripts.
Example – >>> if transcript.is_coding(): … cds_len = transcript.cds_len() … print(f”CDS length: {cds_len}bp”)

test_frames() → bool[source]

Verify that exon frames are correctly set.

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

Returns:

bool – True if the frame of each exon equals the value from calc_frames(), False otherwise. An exon whose frame is None does not match.
Example – >>> transcript.update_frames() >>> assert transcript.test_frames()

total_len() → int[source]

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.

Returns:: The spliced length in base pairs.
Return type:: int

update_frames() → None[source]

Update the frame attribute of all exons.

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

Returns:

None – The method sets the frame of each exon in place.
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[source]

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.

Returns:: The length of the 3’ UTR in base pairs.
Return type:: int

utr3_regions() → list[BedRegion][source]

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[BedRegion] – The 3’ UTR regions. The list is empty for a non-coding transcript.
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[source]

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.

Returns:: The length of the 5’ UTR in base pairs.
Return type:: int

utr5_regions() → list[BedRegion][source]

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[BedRegion] – The 5’ UTR regions. The list is empty for a non-coding transcript.
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 Exon(start: int, stop: int, frame: int | None = None)[source]

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[source]

Check if this exon fully contains a genomic region.

Parameters:

region (tuple[int, int]) – The (start, stop) positions of the region to check (1-based, inclusive).

Returns:

bool – True if the region lies fully within this exon, False otherwise.
Example – >>> exon = Exon(100, 200) >>> exon.contains((150, 160)) True >>> exon.contains((50, 250)) False

property frame: int | None: The codon reading frame of the exon, or None if not computed.

property start: int: The genomic start position of the exon (1-based, inclusive).

property stop: int: The genomic stop position of the exon (1-based, inclusive).