Reference genomes

A ReferenceGenome wraps a genome resource — a FASTA, plain or bgzipped, plus its index — and answers two kinds of question: what chromosomes are there and how long are they, and what are the nucleotides in a given interval.

Opening and closing

The genome holds an open FASTA handle, so it is built closed and must be opened. open() returns the genome itself, so the call chains onto the builder:

from gain.genomic_resources.repository_factory import build_genomic_resource_repository
from gain.genomic_resources.reference_genome import build_reference_genome_from_resource_id

grr = build_genomic_resource_repository()
genome = build_reference_genome_from_resource_id("hg38/genomes/GRCh38-hg38", grr).open()

In anything longer than a script, wrap it in a with block so the handle is released on an exception too. Note that entering the context manager does not open the genome — __enter__ returns it unchanged, and what with adds is the guaranteed close(). Keep the .open():

with build_reference_genome_from_resource_id("hg38/genomes/GRCh38-hg38", grr).open() as genome:
    print(genome.get_sequence("chr21", 5_030_000, 5_030_060))

Dropping the .open() fails less obviously here than it does for a score: the chromosome accessors load the index lazily and so still answer, while get_sequence() asserts. is_open() reports the current state, which is worth checking in library code that may be handed either a fresh or an already-opened genome.

Chromosomes

chromosomes lists the contigs in the order the index gives them, and get_all_chrom_lengths() returns the whole name-to-length mapping in one call. Both read the same lazily-loaded index that get_chrom_length() does, so neither is meaningfully faster than the other; the mapping is worth knowing about because it is the live internal dictionary rather than a copy, so treat it as read-only.

chrom_prefix is the small piece of glue that makes a script portable between builds that spell their contigs chr1 and builds that spell them 1.

Sequence

get_sequence() returns an interval as a string. Coordinates are 1-based and the interval is closed at both ends, matching the convention the rest of GAIn and the underlying formats use:

print(genome.get_sequence("chr21", 5_030_000, 5_030_060))

For intervals large enough that a single string is awkward, fetch() yields the same nucleotides in buffered chunks instead.

split_into_regions() divides the genome into Region pieces of the requested size — the standard way to fan a whole-genome computation out over a task graph or a process pool. The last region of each chromosome is open-ended (Region(chrom, start, None)) rather than padded to the full size, so do not assume every yielded region has an end:

for region in genome.split_into_regions(10_000_000, chromosome="chr21"):
    print(region)

Curator-facing detail — how a genome resource’s genomic_resource.yaml names its FASTA and its index, and what PARS configuration is_pseudoautosomal() reads — is on Genomic resources and repositories.

API

gain.genomic_resources.reference_genome.build_reference_genome_from_resource_id(resource_id: str, grr: GenomicResourceRepo | None = None) → ReferenceGenome

class gain.genomic_resources.reference_genome.ReferenceGenome(resource: GenomicResource)

Provides an interface for quering a reference genome.

property chrom_prefix: str

Return a prefix of all chromosomes of the reference genome.

property chromosomes: list[str]

Return a list of all chromosomes of the reference genome.

close() → None

Close reference genome sequence file-like objects.

fetch(chrom: str, start: int, stop: int | None, buffer_size: int = 512) → Generator[str, None, None]

Yield the nucleotides in a specific region.

While line feed calculation can be inaccurate because not every fetch will start at the start of a line, line feeds add extra characters to read and the output is limited by the amount of nucleotides expected to be read.

get_all_chrom_lengths() → dict[str, int]

Return list of all chromosomes lengths.

get_chrom_length(chrom: str) → int

Return the length of a specified chromosome.

static get_schema() → dict[str, Any]

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

The base resource schema plus filename, index_file, chrom_prefix and PARS.

get_sequence(chrom: str, start: int, stop: int) → str

Return sequence of nucleotides from specified chromosome region.

is_open() → bool

Whether open() has run and close() has not since.

is_pseudoautosomal(chrom: str, pos: int) → bool

Return true if specified position is pseudoautosomal.

open() → ReferenceGenome

Open reference genome resources.

property resource_id: str

The id of the genome resource this object wraps.

split_into_regions(region_size: int, chromosome: str | None = None) → Generator[Region, None, None]

Split the reference genome into regions and yield them.

Can specify a specific chromosome to limit the regions to be in that chromosome only.