SeqRepo

Core Principle: SeqRepo is a disk-backed, immutable sequence repository providing high-throughput random access to reference genomes and transcript sequences by accession.

Why it exists

Validation and normalization need the actual bases. To know whether NM_000546.6:c.215C>G is valid, you must know that base 215 really is a C. To 3'-shift a deletion, you must be able to read the surrounding sequence.

Fetching that from NCBI per variant is slow, rate-limited, and: crucially: not reproducible, because remote records change. SeqRepo solves all three by keeping a versioned local snapshot.

Installing

pip install seqrepo
seqrepo --root-directory /usr/local/share/seqrepo pull
seqrepo --root-directory /usr/local/share/seqrepo list-local-instances

export HGVS_SEQREPO_DIR=/usr/local/share/seqrepo/2024-02-20

Roughly 15 GB for a full human snapshot. Sizeable, but you download it once and it never changes: immutability is the design.

Without HGVS_SEQREPO_DIR, hgvs falls back to fetching sequences remotely (via SeqRepo REST or NCBI), which works and is slow.

Direct use

from biocommons.seqrepo import SeqRepo

sr = SeqRepo("/usr/local/share/seqrepo/2024-02-20")

sr["NM_000546.6"]                        # the full sequence
sr["NM_000546.6"][214:215]               # 0-based slice → the base at c. position 215
sr["NC_000017.11"][7676153:7676154]      # a genomic base
len(sr["NC_000017.11"])

sr.translate_alias("NM_000546.6")        # every known alias for this sequence
sr.translate_identifier("refseq:NM_000546.6")

Slicing is 0-based half-open (Python convention) while HGVS positions are 1-based. Hence [214:215] for HGVS position 215. This is the Coordinate Systems conversion, and it is exactly the kind of thing you should let the library do rather than doing by hand.

The namespace system

SeqRepo keys sequences by digest, not by accession, and treats accessions as aliases:

sr.translate_alias("NM_000546.6")
# [{'namespace': 'refseq',    'alias': 'NM_000546.6'},
#  {'namespace': 'ensembl',   'alias': 'ENST00000269305.9'},
#  {'namespace': 'ga4gh',     'alias': 'SQ.xxxxx...'},
#  {'namespace': 'MD5',       'alias': '...'},
#  ...]

Two consequences worth understanding:

1. Deduplication. RefSeq and Ensembl often distribute byte-identical sequences under different accessions. SeqRepo stores one copy.
2. Content-addressed identity. The ga4gh:SQ. digest identifies the sequence by its content, so two accessions with the same digest are provably the same sequence. This is the foundation of GA4GH VRS: the successor scheme to string-based variant identity.

Immutability

Once a snapshot is created it never changes. New releases are new directories:

/usr/local/share/seqrepo/
    2021-01-29/
    2023-05-15/
    2024-02-20/

This is the property that makes reproducibility possible. A pipeline that records HGVS_SEQREPO_DIR=.../2024-02-20 can be re-run in five years and produce identical output. One that fetches from NCBI cannot.

Use with hgvs

Set the environment variable and hgvs uses it automatically:

export HGVS_SEQREPO_DIR=/usr/local/share/seqrepo/2024-02-20

Or explicitly:

'C'

The fetcher resolution order is: HGVS_SEQREPO_DIR → SeqRepo REST service → NCBI E-utilities. The first is fast and reproducible; the last is neither.

Verifying reference bases

The single most useful thing SeqRepo enables:

A mismatch means one of: the wrong accession version, the wrong assembly (GRCh37 vs GRCh38), a strand error, or a genuinely bad record. All four are common, and all four are silent without this check. hgvs Validator does this for you.

Adding your own sequences

seqrepo --root-directory /path/to/seqrepo init -i my_instance
seqrepo --root-directory /path/to/seqrepo load -i my_instance -n mynamespace custom.fasta

Useful for custom references, patient-specific assemblies, or transcripts absent from RefSeq.

SeqRepo vs the alternatives

	SeqRepo	pysam.FastaFile	Bio.SeqIO.index
Random access	yes	yes	yes
Cross-namespace aliases	yes	no	no
Content-addressed	yes	no	no
Immutable/versioned	yes	no	no
Transcripts + genome together	yes	one file at a time	one file at a time
Ecosystem	biocommons/VRS	samtools	Biopython

For variant work, SeqRepo's advantages are decisive: one lookup interface for transcripts and chromosomes, alias resolution across RefSeq/Ensembl, and immutability. For a plain reference FASTA in a read-alignment pipeline, pysam.FastaFile is simpler.

Disk and setup cost

The 15 GB and the setup are the main objections to SeqRepo, and they are real. The honest cost–benefit:

• Exploring a few variants: skip it. Use remote fetching via hgvs easy.
• Any pipeline, any clinical work, anything you will publish: install it. The reproducibility argument alone justifies the disk, and the speed difference on batch work is large.

Common mistakes

• Not setting HGVS_SEQREPO_DIR, silently falling back to slow remote fetching.
• 1-based/0-based confusion when slicing directly.
• Not recording the snapshot date in your outputs.
• Assuming an accession is present. Older or suppressed versions may not be in a given snapshot.
• Mixing snapshots between development and production.
• Fetching sequences remotely in a loop.
• Pointing HGVS_SEQREPO_DIR at the root rather than at a specific dated instance directory.

See also

UTA and Data Providers · hgvs Validator · hgvs Normalizer · hgvs easy · Coordinate Systems · Applied - Reproducible Environment

Exercise

SeqRepo is not installed in this sandbox, so you will build the store yourself: sequences keyed by digest, accessions held as aliases pointing at those digests, and lookups going through the alias table. Writing it out is the shortest route to understanding why SeqRepo can hold one copy of a sequence that RefSeq and Ensembl both distribute, and where the 1-based to 0-based conversion has to happen.

Write check_refs(text), taking the whole of stdin and returning a list of output lines. Above a line -- stdin gives seq <digest> <sequence> lines and alias <namespace>:<accession> <digest> lines. Below it come queries of the form <namespace>:<accession> <start> <ref>, where start is a 1-based position and ref is the reference bases the caller claims sit there. Emit <query> no_data when the alias is unknown, <query> out_of_range when the claimed interval is not wholly inside the sequence, <query> mismatch:<actual> when the bases differ from the sequence, and <query> ok when they agree.

Finish with queries=<n> and then aliases=<a> digests=<d>, where a counts the alias entries and d counts the distinct stored sequences: the pair that shows the deduplication. Both lines are printed even when no queries follow the --.