Reference genomes

Genome resources provide the DNA sequences of reference genome assemblies, together with the names and lengths of their chromosomes and other contigs. In a GRR, a genome is a resource of type genome backed by a FASTA file (plain or bgzipped) and its index files. GAIn provides the ReferenceGenome Python object for accessing this information. Unlike a generic GenomicResource, a ReferenceGenome provides methods for retrieving sequences and working with genomic coordinates.

Accessing a genome

A ReferenceGenome object can be built from a resource ID and a repository using build_reference_genome_from_resource_id():

import gain.grr as grr

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

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

resource = repo.get_resource("hg38/genomes/GRCh38-hg38")
genome = grr.build_reference_genome_from_resource(resource)

Both builder functions are available through gain.grr.

Opening and closing a genome

A ReferenceGenome is initially closed. You must call open() before retrieving DNA sequences, because the object uses an open FASTA file handle. The open() method returns the genome object itself, so it can be chained to the builder:

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

For scripts and functions, using a with block ensures that the file handle is released when you finish, even if an error occurs:

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

The with statement closes the genome when the block ends, but it does not open it automatically. The .open() call is therefore necessary. Chromosome names and lengths remain accessible even when the genome is closed, but get_sequence() requires the genome to be open.

You can close an opened genome explicitly with close() and check its state using is_open().

Chromosome names and lengths

The chromosomes attribute lists the contigs in the order recorded in the FASTA index. To obtain the length of a particular chromosome, use get_chrom_length(). To retrieve the complete mapping of contig names to lengths, use get_all_chrom_lengths():

print(genome.chromosomes)
print(genome.get_chrom_length("chr21"))
print(genome.get_all_chrom_lengths())

Retrieving the complete mapping does not require reading the genome sequence. The mapping returned by get_all_chrom_lengths() is the object’s internal dictionary rather than a copy, so it should not be modified.

Different genome assemblies may use chromosome names such as chr1 or 1. The chrom_prefix attribute reports the naming prefix determined from the FASTA index, allowing scripts to accommodate these differences.

For assemblies with defined pseudoautosomal regions, is_pseudoautosomal() uses the PARS information in the resource definition. The configuration of genome resources, including their FASTA files, indexes, and pseudoautosomal regions, is described in Genomic resources and repositories.

Retrieving DNA sequences

get_sequence() retrieves the nucleotides in a specified interval as a string. GAIn uses 1-based coordinates with both endpoints included:

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

For large intervals, fetch() yields the sequence one nucleotide at a time, reading the file in buffered chunks, rather than returning one large string. This can be useful when processing long genomic intervals without holding the entire sequence in memory.

Dividing a genome into regions

split_into_regions() divides a genome into Region objects of a requested size. This is useful when distributing a whole-genome analysis across multiple tasks or processes. You can also restrict the division to a particular chromosome:

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

The last region of each chromosome is represented as Region(chrom, start, None), where None indicates that the region extends to the end of the chromosome. Code processing these regions should therefore not assume that every region has an explicit end coordinate.

API

The class below is defined in the module gain.genomic_resources.reference_genome, and so are the builders build_reference_genome_from_resource_id and build_reference_genome_from_resource. ReferenceGenome is also exported by gain.grr, as grr.ReferenceGenome. The builders are also exported by gain.grr, as grr.build_reference_genome_from_resource_id and grr.build_reference_genome_from_resource.

`build_reference_genome_from_resource`(resource)	Build the reference genome of a genomic resource.
`build_reference_genome_from_resource_id`(...)	Build the reference genome of the resource with a given id.
`ReferenceGenome`(resource)	Query the chromosomes and the nucleotides of a reference genome.
`ReferenceGenome.chrom_prefix`	`"chr"` or `""`.
`ReferenceGenome.chromosomes`	The chromosome names, in the order of the `.fai` index.
`ReferenceGenome.close`()	Close the sequence file and drop the loaded index.
`ReferenceGenome.fetch`(chrom, start, stop[, ...])	Yield the nucleotides of a region, one character at a time.
`ReferenceGenome.get_all_chrom_lengths`()	Return the lengths of all chromosomes.
`ReferenceGenome.get_chrom_length`(chrom)	Return the length of a chromosome.
`ReferenceGenome.get_schema`()	Return the schema that the config of a genome is checked against.
`ReferenceGenome.get_sequence`(chrom, start, stop)	Return the nucleotides of a region as one string.
`ReferenceGenome.is_open`()	Whether `open()` has run and `close()` has not since.
`ReferenceGenome.is_pseudoautosomal`(chrom, pos)	Tell whether a position is in a pseudoautosomal region.
`ReferenceGenome.open`()	Open the sequence file of the genome.
`ReferenceGenome.resource_id`	The id of the genome resource this object wraps.
`ReferenceGenome.split_into_regions`(region_size)	Split the genome into regions of a given size.

build_reference_genome_from_resource(resource: GenomicResource) → ReferenceGenome[source]

Build the reference genome of a genomic resource.

Parameters:: resource (GenomicResource) – A resource of type genome.
Returns:: The closed genome. Call open() before any sequence query.
Return type:: ReferenceGenome
Raises:: ValueError – If the resource is not of type genome.

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

Build the reference genome of the resource with a given id.

Parameters:

resource_id (str) – The id of a resource of type genome 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 closed genome. Call open() before any sequence query.

Return type:

ReferenceGenome

Raises:

FileNotFoundError – If the repository has no resource with this id.
ValueError – If the resource is not of type genome.

class ReferenceGenome(resource: GenomicResource)[source]

Query the chromosomes and the nucleotides of a reference genome.

The genome wraps a resource of type genome. The builder returns it closed. Call open() before fetch() or get_sequence(). The chromosome accessors load the index on demand, so they answer without open().

Variables:

resource (GenomicResource) – The genome resource that the object wraps.
config (dict) – The validated configuration of the resource.
pars (dict) – The pseudoautosomal regions from the PARS key of the configuration, as a mapping from the chromosome name to a list of Region. The mapping is empty when the configuration has no PARS key.

property chrom_prefix: str

"chr" or "".

The value is "chr" when the first chromosome name starts with it, and "" otherwise. A chrom_prefix key in the configuration does not change it.

Type:: The prefix of the chromosome names

property chromosomes: list[str]: The chromosome names, in the order of the .fai index.

close() → None[source]

Close the sequence file and drop the loaded index.

Closing a genome that is not open does nothing.

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

Yield the nucleotides of a region, one character at a time.

The coordinates are 1-based and both ends are included. The genome must be open.

Parameters:

chrom (str) – The name of the chromosome.
start (int) – The first position of the region.
stop (int | None) – The last position of the region. A value past the end of the chromosome is cut to the end. When it is None, the region runs to the end of the chromosome.
buffer_size (int) – The number of characters that the reader takes from the file in one step.

Returns:

The nucleotides in upper case. The generator is empty when the genome has no chromosome with this name. In that case a warning is logged.

Return type:

Generator[str, None, None]

get_all_chrom_lengths() → dict[str, int][source]

Return the lengths of all chromosomes.

Returns:: A mapping from the chromosome name to its length, in the order of the .fai index. The mapping is the internal one, not a copy. Do not change it.
Return type:: dict[str, int]

get_chrom_length(chrom: str) → int[source]

Return the length of a chromosome.

Parameters:: chrom (str) – The name of the chromosome.
Returns:: The length of the chromosome in nucleotides.
Return type:: int
Raises:: ValueError – If the genome has no chromosome with this name.

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

Return the schema that the config of a genome is checked against.

Returns:: The base resource schema, with the keys filename, index_file, PARS and chrom_prefix.
Return type:: dict[str, Any]

get_sequence(chrom: str, start: int, stop: int) → str[source]

Return the nucleotides of a region as one string.

The coordinates are 1-based and both ends are included. The genome must be open.

Parameters:

chrom (str) – The name of the chromosome.
start (int) – The first position of the region.
stop (int) – The last position of the region. A value past the end of the chromosome is cut to the end.

Returns:

The nucleotides in upper case. The string is empty when the genome has no chromosome with this name.

Return type:

str

is_open() → bool[source]

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

Returns:: True when the sequence file is open.
Return type:: bool

is_pseudoautosomal(chrom: str, pos: int) → bool[source]

Tell whether a position is in a pseudoautosomal region.

Parameters:

chrom (str) – The name of the chromosome.
pos (int) – The 1-based position.

Returns:

True when the position is in a region of pars. The value is False for a chromosome without such regions, and for every position when the configuration has no PARS key.

Return type:

bool

open() → ReferenceGenome[source]

Open the sequence file of the genome.

Returns:: The genome itself, so that the call chains onto a builder. A genome that is already open is returned unchanged.
Return type:: ReferenceGenome
Raises:: FileNotFoundError – If the index file or the sequence file is missing.

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

Split the genome into regions of a given size.

Parameters:

region_size (int) – The size of each region in nucleotides. The value 0 yields one region for each chromosome.
chromosome (str | None) – The chromosome to split. When it is None, the method splits all chromosomes.

Returns:

The regions, chromosome by chromosome. The last region of each chromosome has stop set to None, and it can be shorter than region_size.

Return type:

Generator[Region, None, None]

Raises:

ValueError – If region_size is negative. The call raises. If chromosome is not a chromosome of the genome, the generator raises on the first step.