gain.genomic_resources.gene_models package

Submodules

gain.genomic_resources.gene_models.default_attributes module

Encoding of the atts column of the default gene models format.

The column packs a transcript’s attributes into one field as key:value pairs joined by ;. Both delimiters, and the escape character itself, are backslash-escaped inside keys and values, so that an attribute value holding either delimiter survives a save/load round trip.

escape_default_attribute(value: str) → str

Escape the attribute delimiters in a key or a value.

format_default_attributes(attributes: dict) → str

Pack a transcript’s attributes into the atts column.

parse_default_attributes(atts: str) → dict[str, str]

Unpack the atts column into the attributes it holds.

unescape_default_attribute(value: str) → str

Reverse escape_default_attribute.

A backslash not followed by a delimiter or another backslash is literal, so free text that happens to carry one – as NCBI RefSeq notes do – reads back unchanged. A backslash directly in front of one of those characters is always taken as an escape: the column records nothing that would tell an escape apart from a literal backslash there.

gain.genomic_resources.gene_models.gene_models module

class GeneModels(resource: GenomicResource)

Bases: ResourceConfigValidationMixin

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

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

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

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]

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]

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]

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

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

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

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

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.

create_regions_from_genes(gene_models: GeneModels, genes: list[str], regions: list[Region] | None, gene_regions_heuristic_cutoff: int = 20, gene_regions_heuristic_extend: int = 20000) → list[Region] | None

Produce a list of regions from given gene symbols.

If given a list of regions, will merge the newly-created regions from the genes with the provided ones.

gain.genomic_resources.gene_models.gene_models_factory module

build_gene_models_from_file(file_name: str, file_format: str | None = None, gene_mapping_file_name: str | None = None, chrom_mapping_file_name: str | None = None) → GeneModels

Load gene models from local filesystem.

build_gene_models_from_resource(resource: GenomicResource | None) → GeneModels

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

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.

gain.genomic_resources.gene_models.parsers module

class ColumnarLayout(accepted_columns: tuple[tuple[str, ...], ...], gene_columns: tuple[str, ...], attribute_columns: tuple[str, ...] = ())

Bases: object

One UCSC-derived columnar gene-models layout.

The five layouts gain reads – refFlat, refSeq, CCDS, knownGene and UCSC genePred – are one record loop over one raw read. They differ on three axes and no others, which are the three fields below (gain#941). Everything else the loop does is the same for all of them and lives in parse_columnar_format: the half-open-to-inclusive coordinate shift, suffixing a transcript name into a unique id, the exon read, and update_frames().

gain’s own output format is deliberately not one of these. It is read by column name rather than by position, is already in gain’s coordinates, carries a third exon column, and builds its attributes by parsing a dedicated column rather than by copying whole cells.

accepted_columns: tuple[tuple[str, ...], ...]

The column lists this format accepts, tried in order. Only genePred has more than one – the ten-column genePred core and the fifteen-column genePredExt. Neither attempt can consume the other’s file, though the two branches of parse_raw rule that out differently: probe_header matches the header against these names and probe_columns counts them.

attribute_columns: tuple[str, ...] = ()

The columns copied into TranscriptModel.attributes, in this order – attributes are written back out in iteration order, so the order is part of the layout. This is the union over the accepted column lists, not a subset any one file carries: it is how genePred’s two widths share one row, the narrow one carrying none of these five. parse_columnar_format narrows it to the width that actually matched, once, before reading any record.

These are pinned to text on the headerless read as well (gain#973), so naming a column here decides its dtype and not only its presence.

gene_columns: tuple[str, ...]

Where the gene label comes from, best candidate first. The first column carrying a non-blank cell wins; when none does, the last column named here supplies its cell anyway, blank or absent or not. A one-column rule is therefore just that column, read unconditionally, which is what four of the five layouts want.

This list is also part of what parse_columnar_format hands parse_raw to pin to text, so adding a column here does more than reorder the fallback: it changes the dtype that column is read at (gain#963).

class FormatInference(matched: tuple[str, ...], rejected: tuple[tuple[str, str], ...], sampled_rows: int, tie_break: str | None = None)

Bases: object

What trying every supported format against a file prefix established.

A format rejects a file through one of two channels: it raises, or it quietly returns no transcript models. Both end up in rejected, so the ledger has no holes – a reason the reader cannot see is the whole of gain#856.

property file_format: str | None

The inferred format, or None when the file stays ambiguous.

Exactly one matching format is an inference; a multi-format collision resolved by a content tie-break is one too.

matched: tuple[str, ...]

rejected: tuple[tuple[str, str], ...]

report() → str

Render why inference did not settle on a single format.

sampled_rows: int

tie_break: str | None = None

GTF_CDS_FEATURES = frozenset({'CDS'})

Features that state the coding sequence itself, one record per coding stretch of an exon. Widened into cds exactly as the codon records are, and wherever a codon record is missing they are the only statement of the extent there is. For a complete transcript they add nothing: GENCODE and Ensembl exclude the stop codon from their CDS records, so the codon span already covers them. NCBI includes it – folding both sources together answers the same under either convention, so this carries no assumption about flavour.

GTF_CODON_FEATURES = frozenset({'start_codon', 'stop_codon'})

Features that delimit the coding sequence. Each record widens its transcript’s cds interval to cover the codon’s span.

GTF_EXONLESS_TRANSCRIPT_FEATURES = frozenset({'miRNA', 'pre_miRNA'})

Transcript-level features FlyBase emits with no exon records at all, so admitting them would add hundreds of transcript models carrying no sequence. This is a policy about these two spellings: it skips them up front, before their children are read, which is what lets a child record be reported against a named skipped transcript. An accepted feature that turns out to have no exon child is a separate matter – it is dropped after the whole file is read (gain#965), so no transcript reaches the models with an empty exon list by either route. Moving one of these into GTF_TRANSCRIPT_FEATURES should be deliberate, not a silent behaviour change.

GTF_EXON_FEATURES = frozenset({'exon'})

Features that append an exon to their transcript’s model.

GTF_IGNORED_FEATURES = frozenset({'3UTR', '5UTR', 'UTR', 'five_prime_utr', 'gene', 'three_prime_utr'})

Features whose records contribute nothing to the models and are skipped outright, before attribute parsing – so an ignored record is not required to carry a transcript_id (Ensembl gene records genuinely lack one). gene restates what every transcript-level record already carries, and the UTR spellings are implied by the exons. The exonless biotypes are deliberately not here: their skip runs after attribute parsing, so their children’s errors can name the skipped transcript.

GTF_SELENOCYSTEINE_FEATURES = frozenset({'Selenocysteine'})

Features that mark a site within their transcript and contribute nothing to the model. GENCODE emits one Selenocysteine record per recoded UGA codon of a selenoprotein. Taking no measurement from them is a redundancy, not a policy: every such site falls inside a CDS record of the same transcript, so cds already covers it – all 130 records across the 88 selenoproteins of GENCODE v49 comprehensive. Dispatched as child records, so that a record with no parent transcript is reported rather than silently turned into a transcript of its own.

GTF_TRANSCRIPT_FEATURES = frozenset({'mRNA', 'ncRNA', 'pseudogene', 'rRNA', 'snRNA', 'snoRNA', 'tRNA', 'transcript'})

Features that introduce a transcript. Ensembl and RefSeq emit the literal transcript; FlyBase instead names the transcript by its biotype. Every entry here is handled identically – it creates a transcript model keyed by transcript_id. Supporting a flavour usually takes more than this set: FlyBase also relies on the 5UTR/3UTR spellings in GTF_IGNORED_FEATURES and on gene_symbol as its gene label. Check a new file’s cut -f3 | sort -u against the module’s GTF_* constants, which between them spell out the whole vocabulary the loop dispatches on. Flavour is an intake concern only – serialization.py normalises back out, always writing transcript and gene_name.

UCSC_GENEPRED_LAYOUT = ColumnarLayout(accepted_columns=(('name', 'chrom', 'strand', 'txStart', 'txEnd', 'cdsStart', 'cdsEnd', 'exonCount', 'exonStarts', 'exonEnds'), ('name', 'chrom', 'strand', 'txStart', 'txEnd', 'cdsStart', 'cdsEnd', 'exonCount', 'exonStarts', 'exonEnds', 'score', 'name2', 'cdsStartStat', 'cdsEndStat', 'exonFrames')), gene_columns=('name2', 'name'), attribute_columns=('score', 'name2', 'cdsStartStat', 'cdsEndStat', 'exonFrames'))

The only layout accepting two widths, and the only one whose gene label has a fallback: the narrow form has no alternate-name column at all, and the wide form may carry a blank one. UCSC’s own genePred and genePredExt table definitions – the sole specification either width has – are quoted in parse_ucscgenepred_models_format.

get_parser(fileformat: str) → Callable[[IO, dict[str, str] | None, int | None], dict[str, TranscriptModel] | None] | None

Get gene models parser based on file format.

infer_gene_models_format(infile: IO) → FormatInference

Try every supported format against a prefix of infile.

load_gene_mapping(resource: GenomicResource) → dict[str, str]

Load alternative names for genes.

Assume that its first line has two column names

load_transcript_models(resource: GenomicResource) → dict[str, TranscriptModel]

Load gene models.

parse_ccds_gene_models_format(infile: IO, gene_mapping: dict[str, str] | None = None, nrows: int | None = None) → dict[str, TranscriptModel] | None

Parse CCDS gene models file format.

A CCDS model’s gene and transcript name are the same string; see _REFSEQ_COLUMNS for why this format and refSeq are two entries.

parse_columnar_format(layout: ColumnarLayout, infile: IO, gene_mapping: dict[str, str] | None = None, nrows: int | None = None) → dict[str, TranscriptModel] | None

Parse a columnar gene-models file against one layout.

Returns None when the file matches none of the layout’s accepted column lists – the GeneModelsParser rejection convention, which format inference reads back to say why a format lost.

Both of the layout’s per-record rules are resolved against the width that matched before any record is read, rather than re-derived per record. That is not only for speed, though it is worth about 250ns a record over files that run to the hundreds of thousands: which columns a record carries is settled by the match, because parse_raw asserts the frame’s columns are exactly the ones it was given.

parse_default_gene_models_format(infile: IO, gene_mapping: dict[str, str] | None = None, nrows: int | None = None) → dict[str, TranscriptModel] | None

Parse default gene models file format.

parse_gtf_gene_models_format(infile: IO, gene_mapping: dict[str, str] | None = None, nrows: int | None = None) → dict[str, TranscriptModel] | None

Parse GTF gene models file format.

parse_known_gene_models_format(infile: IO, gene_mapping: dict[str, str] | None = None, nrows: int | None = None) → dict[str, TranscriptModel] | None

Parse known gene models file format.

parse_raw(infile: IO, expected_columns: list[str], nrows: int | None = None, comment: str | None = None, text_columns: tuple[str, ...] = ()) → DataFrame | None

Parse raw gene models data based on expected columns.

Both branches keep a blank cell as the empty string it was. Letting pandas filter it instead made what a cell became a property of its whole column rather than of itself: one blank re-typed the column, so a well-formed record serialized differently depending on whether some other row was blank, and the blank itself reached serialization as the fabricated token nan (gain#931).

What each branch pins differs, and only because of what it costs. The headered branch has always pinned every column to text. The headerless branch pins _IDENTIFYING_COLUMNS, which is what settles the typing gain#929 left here – a chromosome column of bare digits was handed over as the int 17, and a transcript index keyed by that is unreachable by a query for “17”. The two read the same values either way; the pin decides only how much of the frame is object.

text_columns is how a caller names the rest of what has to reach the model as the file spelled it. A gene label is the second such column, and for the same reason: it keys the gene index, so a gene labelled by a bare digit was unreachable by a lookup for its own name, and which label a record got depended on which branch below recognised the file (gain#963). Only the caller knows which column that is – it differs per layout, and for one of them it is the alternate name with the transcript name behind it. gain’s own output format, which does not come through here, pins its gene column the same way and always has.

A layout’s attribute columns are the third, and the first named for something other than a key: they are copied into the model whole and written back out, so leaving them inferred made a record’s attribute a property of which branch read the file. It is the one axis of this where two published resources actually disagreed – refSeq’s #bin, score and exonCount arrive as ints from the headerless file and as text from the headered one (gain#973).

Pinning them rewrites no published resource: both refSeq files serialize to the same bytes either way, over their full length. It is not output-neutral in general though, and must not be – a score of 007 inferred as 7 serializes as 7, and gain’s own format, having only text to write an attribute as, leaves nothing downstream able to say what the file first recorded.

The GTF reader shares this and names none of these columns, so it keeps the inference its own arithmetic depends on. It does not escape na_filter: a blank cell reaches it as '' too, which is what its blank-attributes guard now decides on.

parse_ref_flat_gene_models_format(infile: IO, gene_mapping: dict[str, str] | None = None, nrows: int | None = None) → dict[str, TranscriptModel] | None

Parse refFlat gene models file format.

parse_ref_seq_gene_models_format(infile: IO, gene_mapping: dict[str, str] | None = None, nrows: int | None = None) → dict[str, TranscriptModel] | None

Parse refSeq gene models file format.

parse_ucscgenepred_models_format(infile: IO, gene_mapping: dict[str, str] | None = None, nrows: int | None = None) → dict[str, TranscriptModel] | None

Parse UCSC gene prediction models file fomrat.

table genePred “A gene prediction.”

(
string  name;               "Name of gene"
string  chrom;              "Chromosome name"
char[1] strand;             "+ or - for strand"
uint    txStart;            "Transcription start position"
uint    txEnd;              "Transcription end position"
uint    cdsStart;           "Coding region start"
uint    cdsEnd;             "Coding region end"
uint    exonCount;          "Number of exons"
uint[exonCount] exonStarts; "Exon start positions"
uint[exonCount] exonEnds;   "Exon end positions"
)

table genePredExt “A gene prediction with some additional info.”

(
string name;            "Name of gene (usually transcript_id from
                            GTF)"
string chrom;           "Chromosome name"
char[1] strand;         "+ or - for strand"
uint txStart;           "Transcription start position"
uint txEnd;             "Transcription end position"
uint cdsStart;          "Coding region start"
uint cdsEnd;            "Coding region end"
uint exonCount;         "Number of exons"
uint[exonCount] exonStarts; "Exon start positions"
uint[exonCount] exonEnds;   "Exon end positions"
int score;              "Score"
string name2;           "Alternate name (e.g. gene_id from GTF)"
string cdsStartStat;    "Status of CDS start annotation (none,
                            unknown, incomplete, or complete)"
string cdsEndStat;      "Status of CDS end annotation
                            (none, unknown,
                            incomplete, or complete)"
lstring exonFrames;     "Exon frame offsets {0,1,2}"
)

probe_columns(infile: IO, expected_columns: list[str], comment: str | None = None) → bool

Probe gene models file based on expected columns.

probe_header(infile: IO, expected_columns: list[str], comment: str | None = None) → bool

Probe gene models file header based on expected columns.

read_gene_models_tsv(infile: IO, **kwargs: Any) → DataFrame

Read a gene-models table, keeping every cell as its own text.

Every read that builds records goes through here, so that na_filter is off in one place rather than five. pandas otherwise reads a blank cell – and several spellings that are not blank, NA and NULL among them – as a float NaN, which reaches serialization as the fabricated token nan and re-types the column around it (gain#931).

probe_header and probe_columns do not: they read one row to recognise a layout and never look at a value, so what a blank cell becomes there cannot reach a record.

That the setting had to be repeated per call site is how a read got missed: the gene mapping kept filtering long after the two model reads stopped, and wrote nan into the gene column.

What each caller pins with dtype still differs, and is theirs to decide – the layouts do not agree on which columns are text.

gain.genomic_resources.gene_models.record_cells module

Reading one cell of a columnar gene-models record, and refusing it.

The gene-models parsers all face the same problem: a cell arrives as whatever pandas made of it, and a record built from a cell that cannot be read is worse than no record at all. What a bad cell should produce is one thing – a ValueError naming the record and the column – and it is gathered here rather than repeated in each parser (gain#907, gain#929).

Two shapes of message live here, because a record is named by two of its own columns:

• Once the transcript name and chromosome are known, everything else is reported against them: transcript NM_000546 at chr17 has ....

• Those two cannot be named that way themselves, so each falls back to the record’s position in the file, plus whichever of the pair is readable: gene models record 2 at chr17 has a blank name column.

The layouts that use these are near-duplicates of one another, and driving them from a table instead is gain#941; this module is what such a table would call.

QUOTED_TEXT_LIMIT = 60

How much of a cell to quote back when reporting it. Shared with _scan_gtf_attributes, so that the two messages truncate alike.

cell_text(value: Any) → str

Render what pandas made of a cell as the text the file held.

Since gain#931 the reads keep a cell as its own text, so what the file said is what arrives: a blank cell is '', and the spellings pandas would otherwise have taken for missing values – NA, NULL and nan among them – are the words they are, and the messages built from this can say which one it read.

No read in this module produces a NaN any more – a row that stops short of a column yields '' as well, on both the columnar and the GTF paths (measured on pandas 3.0.2). The pd.isna branch is kept anyway: it costs one comparison on a path that is already building an error message, parse_coordinate relies on it to tell a missing number from one it should convert, and the supported pandas range is wider than the version this was measured on.

value is annotated Any rather than object because pd.isna has no overload for the latter.

parse_coordinate(value: Any, column: str, tr_name: object, chrom: object) → int

Read a single coordinate column, naming its record.

The columnar layouts already wrapped these in int(), which does reject a blank cell – but as cannot convert float NaN to integer, naming neither the record nor the column, and the gain#856 ledger then offers that to the reader as the reason a format was rejected. The default format did not convert at all, so a blank cell became a transcript bound of NaN (gain#929).

A coordinate spelled 100.0 is read as 100 whether it arrives as text or as a number. It used to be only the latter: a column spelled that way throughout was inferred as float on the headerless path and int(100.0) kept it parsing, while the headered path pinned it to text and int("100.0") did not, so the same file parsed or failed depending on which branch recognised it. Since gain#931 reads every columnar cell as text, the text conversion is the one that has to accept both spellings.

OverflowError is caught alongside the rest because inf is a coordinate pandas accepts and int() will not take. It reaches here only on the read path that infers a float column, so without it the two paths report the same file differently – and the one that escaped named neither the record nor the column.

This runs four times per record on files reaching into the hundreds of thousands, so neither path builds a message until there is one to build. Since gain#931 the text path is the common one – the five UCSC-derived layouts pin every column to a string dtype – and the numeric path is reached by the default format’s four bound columns, which are the ones left to inference.

parse_exon_bounds(rec: dict, tr_name: object, chrom: object) → tuple[list[int], list[int]]

Read the paired exon-position columns of a columnar record.

Every columnar layout but the default one spells the pair the same way, so they share this rather than repeating the pair of reads and the length check between them.

parse_exon_positions(value: Any, column: str, tr_name: object, chrom: object) → list[int]

Read a comma-separated coordinate column, naming its record.

pandas delivers a blank cell as a float NaN, which used to reach str.strip and escape as an AttributeError naming a float (gain#907). Text that is simply not a coordinate list fails int the same way, and leaves the reader just as stuck, so both are reported here as one thing: this record’s column could not be read.

Where a GTF record has to be placed by feature and position – its transcript_id being what tends to be missing – a columnar record is named by the transcript name and chromosome every columnar layout carries in columns of their own.

The quoted cell is what pandas made of the column, not the file’s own bytes – see cell_text. The int failure stays on the chain, so the offending token survives the truncation.

parse_transcript_bounds(rec: dict, tr_name: object, chrom: object) → tuple[tuple[int, int], tuple[int, int]]

Read the transcript and coding bounds of a columnar record.

The five UCSC-derived layouts spell these four columns the same way and share the half-open convention that shifts each start by one, so they share this rather than repeating it between them (gain#941).

record_identity(rec: dict, record: int, name_column: str, chrom_column: str) → tuple[Any, Any]

Read the two columns that say which record a columnar row is.

Both used to be taken as they came. A blank one became a float NaN in the model: a NaN chromosome keys the transcript index all by itself, so the record is unreachable by every location query, and a NaN transcript name reaches serialization as the literal token nan – and is suffixed into a transcript id of nan_1, an identifier no file ever carried (gain#929).

These two are what every other message here names a record by, so they cannot be named that way themselves. Each falls back to the record’s position in the file, plus whichever of the pair is still readable.

Blankness is decided on the text, but what is returned is the cell pandas handed over, untouched – a guard meant only to reject must not re-type a cell that parses. It used to matter which read path handed it over: the headerless one inferred a dtype, so a chromosome column of bare digits came back as the int 17, and a transcript index keyed by that is unreachable by a query for “17”. gain#931 settled that at the read boundary, where it belonged, so both paths now hand over text.

require_cell(value: Any, column: str, tr_name: object, chrom: object) → Any

Read a load-bearing text column of an already-identified record.

Blank is refused rather than carried: a strand reaches update_frames(), so a record without one does not merely hold an odd value – its exon frames come out as if it had a strand, and the output is quietly wrong rather than missing (gain#929).

Blankness is decided on the text, but what is returned is the cell pandas handed over, untouched – see record_identity.

require_equal_exon_counts(tr_name: object, chrom: object, **columns: list[int]) → None

Refuse a record whose exon columns disagree on how many exons.

This was a bare assert, which carries no message; the gain#856 ledger renders whatever a parser raised, so what reached the reader as the reason a format was rejected was AssertionError (no message). Naming the record and the counts costs nothing and is the whole of what the reader needed.

The counts themselves are only tallied once they disagree: this runs once per record, and the records that reach it agree.

unparsable(column: str, tr_name: object, chrom: object, text: str) → ValueError

Report a cell that a record cannot be built from.

gain.genomic_resources.gene_models.serialization module

DEFAULT_FORMAT_COLUMNS = ('chr', 'trID', 'trOrigId', 'gene', 'strand', 'tsBeg', 'txEnd', 'cdsStart', 'cdsEnd', 'exonStarts', 'exonEnds', 'exonFrames', 'atts')

The header row of the default format, in column order.

build_gtf_record(transcript: TranscriptModel, feature: str, start: int, stop: int, attrs: str) → tuple[tuple[str, int, int, int], str]

Build an indexed GTF format record for a feature.

calc_frame_for_gtf_cds_feature(transcript: TranscriptModel, region: BedRegion) → int

Calculate frame for the given feature.

collect_gtf_cds_regions(strand: str, cds_regions: list[BedRegion]) → list[BedRegion]

Returns list of all regions that represent the CDS.

collect_gtf_start_codon_regions(strand: str, cds_regions: list[BedRegion]) → list[BedRegion]

Returns list of all regions that represent the start codon.

collect_gtf_stop_codon_regions(strand: str, cds_regions: list[BedRegion]) → list[BedRegion]

Returns list of all regions that represent the stop codon.

find_exon_cds_region_for_gtf_cds_feature(transcript: TranscriptModel, region: BedRegion) → tuple[Exon, BedRegion]

Find exon and CDS region that contains the given feature.

gene_models_to_gtf(gene_models: GeneModels, *, sort_by_position: bool = True) → StringIO

Output a GTF format string representation.

Refuses a model nobody loaded; a loaded, empty model serializes to an empty string, with a warning.

get_exon_number_for(transcript: TranscriptModel, start: int, stop: int) → int

Get the exon number for a genomic region.

Returns the exon number (in transcript order) that overlaps the given genomic coordinates.

Parameters:

• start (int) – Start position (1-based).

• stop (int) – End position (1-based).

Returns:

Exon number (1-based) in transcript orientation.

Returns 0 if no overlapping exon found.

Return type:

int

Example

>>> # For a region within the second exon of a + strand transcript
>>> exon_num = transcript.get_exon_number_for(1000, 1050)
>>> print(f"Region is in exon {exon_num}")

Note

Exon numbering is strand-aware: - Positive strand: numbered 5’ to 3’ (exon 1 is first) - Negative strand: numbered 5’ to 3’ (exon 1 is last in genome)

gtf_canonical_index(index: tuple[str, int, int, int]) → tuple

save_as_default_gene_models(gene_models: GeneModels, output_filename: str, *, gzipped: bool = True) → None

Save gene models in a file in default file format.

Refuses a model nobody loaded, and two transcript shapes the format has no spelling for, before the output file is opened. A loaded, empty model is written as the header row alone, with a warning.

transcript_to_gtf(transcript: TranscriptModel) → list[tuple[tuple[str, int, int, int], str]]

Output an indexed list of GTF-formatted features of a transcript.

gain.genomic_resources.gene_models.to_gpf_gene_models_format module

main(argv: list[str] | None = None) → None

Convert gene models to default GPF gene models format.

gain.genomic_resources.gene_models.transcript_models module

class Exon(start: int, stop: int, frame: int | None = None)

Bases: object

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]) – 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).

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)

Bases: object

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, 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]

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

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]

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

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

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

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

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

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]

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

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]

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.

Module contents