Add optional chromosome FASTA support to Genome (third FASTA tier)

[iskandr]

Motivation

pyensembl currently ships transcript and protein FASTAs but not the chromosome FASTA, so any consumer that needs raw bases at arbitrary genomic positions (intronic, intergenic, flanking) has no upstream path.

Downstream OpenVax tools have repeatedly hit this — varcode in particular has genome.genome.sequence(...) calls inside try/except blocks that silently bail out (varcode/cryptic_exons.py:326), because the prior code expected the API would exist. The "wishful API call" is itself a signal: the natural shape lives here.

Use cases blocked today:

• varcode left-alignment of indels in intronic positions (openvax/varcode#369).
• varcode cryptic-exon scoring for breakpoints outside annotated exons (openvax/varcode#372).
• varcode sequence-aware splice prediction (openvax/varcode#297).
• isovar / vaxrank intronic-retention protein computation (the splice-outcome stubs originally tracked in openvax/varcode#296, now superseded).
• Any downstream consumer that wants Genome.sequence(contig, start, end).

Proposed API

Mirrors the existing transcript/protein FASTA pattern:

release = EnsemblRelease(81, download_genome_fasta=True)
release.sequence("7", 117_480_000, 117_480_100)  # + strand bases, 1-based inclusive

Construction-time opt-in keeps the default install lean (the chromosome FASTA is large; most pyensembl users don't need it). When omitted or False, sequence() either raises a clear MissingGenomeFastaError or returns "" — design choice flagged below.

Parallels to existing slots:

Existing	New
_get_transcript_fasta_paths	_get_genome_fasta_paths
transcript_fasta_urls	genome_fasta_urls
requires_transcript_fasta	requires_genome_fasta
transcript_sequence(transcript_id)	sequence(contig, start, end)

Where to download

Ensembl ships chromosome FASTAs at:

http://ftp.ensembl.org/pub/release-{N}/fasta/{species}/dna/

Three flavors per release / species:

File	Contents	Approx size (human)
{Species}.{Assembly}.dna.primary_assembly.fa.gz	Primary chromosomes only (chr1-22, X, Y, MT). Recommended default.	~830MB compressed, ~3GB uncompressed
{Species}.{Assembly}.dna.toplevel.fa.gz	Primary chromosomes + alt scaffolds + decoys.	~880MB compressed
{Species}.{Assembly}.dna.chromosome.{N}.fa.gz	Per-chromosome split. Useful for partial fetches.	varies

Plus masking variants: dna (unmasked), dna_sm (soft-masked, lowercase for repeats), dna_rm (hard-masked, N's for repeats). Default should be dna unmasked — consumers can soft-mask themselves if needed, and effect-prediction tools usually want the raw bases.

Recommend: download dna.primary_assembly.fa.gz by default; provide --toplevel and --masked={none,soft,hard} flags for opt-in alternatives. Per-chromosome split is implementation detail — pyfaidx can index the merged FASTA fine.

Sharing sequence across Ensembl releases

The chromosome sequence is stable within a major assembly. GRCh38.p13 and GRCh38.p14 have identical bases on chr1-22, X, Y, MT — patches add alt scaffolds and fix annotation but don't change primary chromosome bases. Releases 81 through ~109 all use GRCh38.p13. Re-downloading the same ~3GB for each release is wasteful.

Proposed: cache keyed on (assembly, patch), not release number

Storage:

~/.pyensembl/dna_cache/
    GRCh38.p13/
        Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz
        Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz.fai
    GRCh37.p13/
        Homo_sapiens.GRCh37.dna.primary_assembly.fa.gz
        ...

Each EnsemblRelease symlinks (or refers via a metadata file) into the shared cache. Reinstalling release 82 when release 81 is already installed for GRCh38.p13 reuses the existing FASTA.

Mapping release → patch level

Ensembl's release announcements list the patch level (e.g., release 81 uses GRCh38.p3, release 95 uses GRCh38.p12, release 110 uses GRCh38.p14). Two options:

1. Hard-coded table in pyensembl.species or a new pyensembl.assembly_patches module:

   release_to_assembly_with_patch = {
       (81, "homo_sapiens"): "GRCh38.p3",
       (82, "homo_sapiens"): "GRCh38.p3",
       ...
       (110, "homo_sapiens"): "GRCh38.p14",
   }

   Pros: deterministic, offline, easy to audit. Cons: needs updating every release.

2. Read from FTP at install time: parse the assembly version from README or CHECKSUMS on the release directory. Pros: always current. Cons: needs network even for cached operations.

Recommend (1) with periodic auto-update: ship the static mapping in the package, run a CI job nightly that opens a PR if Ensembl has released a new version. Most users get fast offline behavior; the table stays fresh.

What about alt scaffolds and patches?

Patches do affect alt scaffolds and decoys (a new alt added in p14 isn't in p13). For consumers that only care about primary chromosomes, sharing-by-major-assembly is safe. For consumers that touch alt contigs, sharing should be by full (assembly, patch) tuple.

Recommend: default cache key is the full (assembly, patch) string. Within one patch, all releases share. Across patches, separate FASTAs (you can have GRCh38.p13 and GRCh38.p14 side by side).

A --shared-major-only flag can opt in to coarser sharing for users who never touch alt contigs.

Local versions of the genome

Users with custom FASTAs (curated decoys, organism-specific tweaks, BSgenome dumps, refgenie assets, Galaxy genomes):

release = EnsemblRelease(81, genome_fasta_path="/path/to/custom_GRCh38.fa")

Or via CLI:

pyensembl install --release 81 --genome-fasta-path /path/to/custom_GRCh38.fa

When a path is provided:

1. Skip the download — point pyensembl's _get_genome_fasta_paths at the user-supplied path.
2. Verify contigs — read contig names from the user FASTA and compare against the names pyensembl expects from the GTF; warn on mismatch (helps catch GRCh37 attached to GRCh38, or chr-prefixed-vs-not contigs).
3. Store metadata alongside the path so subsequent pyensembl list shows it as user-attached, not downloaded.
4. Don't cache under the shared (assembly, patch) tree — a user-attached FASTA might be patched / customized and shouldn't pollute the canonical cache.

Data management

CLI additions:

# Opt-in download alongside transcript/protein FASTAs.
pyensembl install --release 81 --with-genome-fasta

# Just the genome FASTA, assuming the rest is already installed.
pyensembl install --release 81 --only-genome-fasta

# Local path attachment (no download).
pyensembl install --release 81 --genome-fasta-path /path/to/custom.fa

# Inspect what's installed.
pyensembl list                    # shows whether genome FASTA is present
pyensembl list --check-genome-fasta  # verifies file exists, fai index valid

# Disk cleanup.
pyensembl prune --orphan-genome-fastas   # remove FASTAs no installed release references

Python API mirrors:

release.download(genome_fasta=True)      # download missing parts
release.install(genome_fasta=True)       # download + index
release.requires_genome_fasta            # False by default
release.genome_fasta_path                # property; None if not installed

Open design questions

1. Default for download_genome_fasta: I'd keep False (opt-in). ~3GB is a lot to surprise people with.
2. What sequence() returns when no FASTA is installed: raise MissingGenomeFastaError (subclass of ValueError so consumers catching ValueError keep working). Returning "" silently is the kind of silent-degradation bug that originally motivated this work.
3. pyfaidx vs samtools faidx vs in-tree FASTA reader: pyfaidx is already a transitive dep through several OpenVax tools and is pure Python. Recommend pyfaidx for the reader; the FASTA is opened lazily on first sequence() call.
4. Soft-masking handling: download unmasked by default. Provide release.sequence(contig, start, end, mask="upper") to uppercase soft-masked bases automatically, with "raw" for verbatim bytes.

What I can offer

This issue is the consequence of a varcode design discussion (openvax/varcode#372, openvax/varcode#373). varcode is shipping a stopgap varcode.Genome wrapper that adds the same API surface (varcode.Genome(81, fasta="/path/...")) at the consumer layer; once pyensembl natively supports this, the wrapper becomes a one-line delegate or evaporates entirely.

I'm happy to take a stab at the implementation here if there's interest — the slots parallel transcript/protein FASTA handling closely, and I've already worked through the lookup logic on the varcode side. Just say the word.