RefgetStore Tutorial

This notebook demonstrates how to use the RefgetStore class from the gtars.refget module for managing reference genome sequences using the GA4GH refget standard.

Features Covered

1. Creating and populating a local RefgetStore
2. Retrieving sequences by digest and by collection + name
3. Getting substrings and BED file regions
4. Exporting FASTA from a sequence collection
5. Saving and loading stores
6. Importing many FASTA files into a disk-backed store
7. Opening remote RefgetStores with lazy loading
8. Caching remote sequences locally

Setup

First, import the necessary modules:

import os
import tempfile
from gtars.refget import RefgetStore, digest_fasta

Part 1: Local RefgetStore - Full Workflow

This section demonstrates the complete workflow for creating and using a local RefgetStore.

1.1 Create a Sample FASTA File

We'll create a temporary FASTA file with two chromosomes:

# Create a temporary directory for our example and work inside it,
# so all the paths below are short and relative
original_dir = os.getcwd()
temp_dir = tempfile.mkdtemp()
os.chdir(temp_dir)
print("Working in a temporary directory")

# Create a sample FASTA file
fasta_content = (
    ">chr1\n"
    "ATGCATGCATGCAGTCGTAGC\n"
    ">chr2\n"
    "GGGGAAAA\n"
)
source_fasta_path = "source.fa"
with open(source_fasta_path, "w") as f:
    f.write(fasta_content)

print(f"Created FASTA file: {source_fasta_path}")

Working in a temporary directory
Created FASTA file: source.fa

1.2 Digest the FASTA to Get Collection Information

The digest_fasta function computes GA4GH-compliant digests for each sequence:

# Digest the FASTA to get collection info
collection = digest_fasta(source_fasta_path)
collection_digest = collection.digest

print(f"Collection digest: {collection_digest}")
print(f"Number of sequences: {len(collection)}")
print(f"\nSequences in collection:")
for seq in collection:
    print(f"  - {seq.metadata.name}: {seq.metadata.sha512t24u}")

Collection digest: aep4zTAZpy0wOIfWEQ5o8oykl9M47YVJ
Number of sequences: 2

Sequences in collection:
  - chr1: k4ZQ03Uo81xEjENxmvAkcLnTez9OFlQE
  - chr2: OyXJErGtjgcIVSdobGkHE3sBdQ5faDTf

1.3 Initialize and Populate a RefgetStore

Create a store. You can choose between in-memory or disk-backed:

• RefgetStore.in_memory(): Keeps all sequences in RAM for fast access (shown below)
• RefgetStore.on_disk(path): Writes sequences to disk as they are added, only keeps metadata in RAM (memory-efficient; see section 1.11)

Both start in Encoded storage mode, which packs DNA at 2 bits per base.

# Initialize an in-memory store (Encoded mode by default)
store = RefgetStore.in_memory()
print(f"Initialized store: {store!r}")
print(f"Storage mode: {store.storage_mode}")

# Import the FASTA file into the store.
# This returns the collection metadata and whether it was newly added.
metadata, was_new = store.add_sequence_collection_from_fasta(source_fasta_path)
print(f"\nImported collection {metadata.digest} ({metadata.n_sequences} sequences), new: {was_new}")

# Importing the same file again is a no-op
_, was_new_again = store.add_sequence_collection_from_fasta(source_fasta_path)
print(f"Imported again, new: {was_new_again}")

print(f"\nStore now contains: {store!r}")

Initialized store: RefgetStore(n_sequences=0, n_sequences_in_memory=0, n_collections=0, n_collections_in_memory=0, mode=Encoded, persist=off, quiet=off, memory-only)
Storage mode: StorageMode.Encoded

Imported collection aep4zTAZpy0wOIfWEQ5o8oykl9M47YVJ (2 sequences), new: True
Imported again, new: False

Store now contains: RefgetStore(n_sequences=2, n_sequences_in_memory=2, n_collections=1, n_collections_in_memory=1, mode=Encoded, persist=off, quiet=off, memory-only)
Processing source.fa...
Added aep4zTAZpy0wOIfWEQ5o8oykl9M47YVJ (2 seqs) from source.fa in 0.0s
Imported 1 file(s) in 0.0s (jobs=1)
Processing source.fa...
Skipped aep4zTAZpy0wOIfWEQ5o8oykl9M47YVJ (already exists)
Imported 1 file(s) in 0.0s (jobs=1)

1.4 Retrieve Sequences by Digest or by Name

Get a complete sequence record using its digest with get_sequence(), or look it up by name within a collection with get_sequence_by_name():

# Get the digest for chr1
seq_digest_chr1 = collection[0].metadata.sha512t24u

# Retrieve the sequence record by digest
record_chr1 = store.get_sequence(seq_digest_chr1)

print(f"Retrieved sequence: {record_chr1.metadata.name}")
print(f"  Length: {record_chr1.metadata.length}")
print(f"  Alphabet: {record_chr1.metadata.alphabet}")
print(f"  Sequence: {record_chr1.decode()}")

# Retrieve a sequence by collection digest + sequence name
record_chr2 = store.get_sequence_by_name(collection_digest, "chr2")
print(f"\nchr2 by name: {record_chr2.decode()}")

Retrieved sequence: chr1
  Length: 21
  Alphabet: dna2bit
  Sequence: ATGCATGCATGCAGTCGTAGC

chr2 by name: GGGGAAAA

1.5 Get Substrings

Extract specific regions from sequences:

# Get a substring from chr1 (positions 5-15)
sub_seq = store.get_substring(seq_digest_chr1, 5, 15)
print(f"Substring chr1[5:15]: {sub_seq}")

# Get a substring from chr2
seq_digest_chr2 = collection[1].metadata.sha512t24u
sub_seq_chr2 = store.get_substring(seq_digest_chr2, 0, 4)
print(f"Substring chr2[0:4]: {sub_seq_chr2}")

Substring chr1[5:15]: TGCATGCAGT
Substring chr2[0:4]: GGGG

1.6 Retrieve Regions from BED Files

Extract multiple regions specified in a BED file:

# Create a BED file with regions of interest
bed_content = (
    "chr1\t0\t10\n"
    "chr2\t2\t6\n"
    "chr_nonexistent\t0\t5\n"  # This will be skipped
)
bed_path = "regions.bed"
with open(bed_path, "w") as f:
    f.write(bed_content)

# Retrieve sequences as a list
retrieved_list = store.substrings_from_regions(collection_digest, bed_path)

print("Retrieved sequences from BED file:")
for rs in retrieved_list:
    print(f"  {rs.chrom_name}[{rs.start}-{rs.end}]: {rs.sequence}")

Retrieved sequences from BED file:
  chr1[0-10]: ATGCATGCAT
  chr2[2-6]: GGAA


Line 4: sequence 'chr_nonexistent' not found in collection 'aep4zTAZpy0wOIfWEQ5o8oykl9M47YVJ': Sequence 'chr_nonexistent' not found in collection

1.7 Write Retrieved Sequences to FASTA

Save BED file regions as a new FASTA file. Unlike substrings_from_regions(), which skips regions on unknown sequences, export_fasta_from_regions() stops with an error, so we give it only valid regions:

# A BED file with only regions on sequences in the collection
valid_bed_path = "valid_regions.bed"
with open(valid_bed_path, "w") as f:
    f.write("chr1\t0\t10\n" "chr2\t2\t6\n")

# Write retrieved sequences to a FASTA file
output_fasta_path = "output_regions.fa"
store.export_fasta_from_regions(collection_digest, valid_bed_path, output_fasta_path)

print(f"Sequences written to: {output_fasta_path}\n")
with open(output_fasta_path, "r") as f:
    print("Content:")
    print(f.read())

Sequences written to: output_regions.fa

Content:
>chr1:0-10
ATGCATGCAT
>chr2:2-6
GGAA

1.8 Export a Complete FASTA from a Sequence Collection

You can reconstruct the entire FASTA file for a sequence collection with export_fasta(). Use sequence_names=None to export all sequences, or pass a list of names to export only some:

# Export the whole collection
extracted_fasta_path = "extracted_complete.fa"
store.export_fasta(collection_digest, extracted_fasta_path, sequence_names=None, line_width=80)

print(f"Extracted complete FASTA to: {extracted_fasta_path}")
with open(extracted_fasta_path, "r") as f:
    print("\nContent:")
    print(f.read())

# Export just chr2
chr2_fasta_path = "chr2_only.fa"
store.export_fasta(collection_digest, chr2_fasta_path, sequence_names=["chr2"], line_width=80)
with open(chr2_fasta_path, "r") as f:
    print("chr2 only:")
    print(f.read())

Extracted complete FASTA to: extracted_complete.fa

Content:
>chr1
ATGCATGCATGCAGTCGTAGC
>chr2
GGGGAAAA

chr2 only:
>chr2
GGGGAAAA

Alternative: You can also extract sequences programmatically:

# Extract all sequences from a collection programmatically
manual_fasta_path = "manual_extraction.fa"

with open(manual_fasta_path, "w") as f:
    for seq in collection:
        # Get the sequence data
        seq_data = store.get_substring(
            seq.metadata.sha512t24u,
            0,
            seq.metadata.length
        )

        # Write FASTA header and sequence
        f.write(f">{seq.metadata.name}\n")
        f.write(f"{seq_data}\n")

print(f"Manually extracted FASTA to: {manual_fasta_path}")
with open(manual_fasta_path, "r") as f:
    print("\nContent:")
    print(f.read())

Manually extracted FASTA to: manual_extraction.fa

Content:
>chr1
ATGCATGCATGCAGTCGTAGC
>chr2
GGGGAAAA

1.9 Save Store to Disk

Persist the store for later use:

# Save the store to a directory
saved_store_path = "my_refget_store"
store.write_store_to_dir(saved_store_path, "sequences/%s2/%s.seq")

print(f"Store saved to: {saved_store_path}")
print(f"\nStore structure:")
for root, dirs, files in os.walk(saved_store_path):
    level = root.replace(saved_store_path, '').count(os.sep)
    indent = ' ' * 2 * level
    print(f"{indent}{os.path.basename(root)}/")
    subindent = ' ' * 2 * (level + 1)
    for file in files:
        print(f"{subindent}{file}")

Store saved to: my_refget_store

Store structure:
my_refget_store/
  collections.rgci
  sequences.rgsi
  rgstore.json
  sequences/
    k4/
      k4ZQ03Uo81xEjENxmvAkcLnTez9OFlQE.seq
    Oy/
      OyXJErGtjgcIVSdobGkHE3sBdQ5faDTf.seq
  collections/
    aep4zTAZpy0wOIfWEQ5o8oykl9M47YVJ.rgsi


Writing store to directory: my_refget_store; Using seqdata path template: sequences/%s2/%s.seq

1.10 Open a Store from Disk

Open a previously saved store with open_local(). Only the lightweight index files are read up front; sequence data is loaded on demand:

# Open the store (with lazy loading)
loaded_store = RefgetStore.open_local(saved_store_path)

print(f"Store opened from: {saved_store_path}")
print(f"Loaded store: {loaded_store!r}")

# Verify we can retrieve sequences (data read on demand)
seq = loaded_store.get_substring(seq_digest_chr1, 0, 10)
print(f"\nRetrieved sequence from loaded store: {seq}")

Store opened from: my_refget_store
Loaded store: RefgetStore(n_sequences=2, n_sequences_in_memory=0, n_collections=1, n_collections_in_memory=0, mode=Encoded, persist=on, quiet=off, cache='my_refget_store')

Retrieved sequence from loaded store: ATGCATGCAT

1.11 Import Many FASTA Files into a Disk-Backed Store

For larger jobs, use a disk-backed store and import many FASTA files at once with add_sequence_collections_from_fastas(). It accepts a path, a glob, a directory, or a list of these, and returns an ImportReport describing what the import actually did. Sequences are content-addressed, so a sequence shared between files is stored only once:

# A second FASTA that shares chr2 with the first one
with open("other.fa", "w") as f:
    f.write(">chrA\nTTTTCCCCGGGG\n>chr2\nGGGGAAAA\n")

disk_store = RefgetStore.on_disk("disk_store")
report = disk_store.add_sequence_collections_from_fastas(["source.fa", "other.fa"])

print(f"Collections added: {report.n_collections_new}")
print(f"Sequences written: {report.n_sequences_written}")
print(f"Sequences deduplicated: {report.n_sequences_deduped}")
for meta, was_new in report.collections:
    print(f"  {meta.digest}: {meta.n_sequences} sequences, new: {was_new}")

Processing source.fa...
Processing other.fa...
Collections added: 2
Sequences written: 3
Sequences deduplicated: 1
  aep4zTAZpy0wOIfWEQ5o8oykl9M47YVJ: 2 sequences, new: True
  TEZrSPigoLfsjqAC6wszPlHDFq3nZvLt: 2 sequences, new: True


Added TEZrSPigoLfsjqAC6wszPlHDFq3nZvLt (2 seqs) from other.fa in 0.0s
Added aep4zTAZpy0wOIfWEQ5o8oykl9M47YVJ (2 seqs) from source.fa in 0.0s
Imported 2 file(s) in 0.1s (jobs=2)

stats() is different: it is a snapshot of what the store holds right now, including how much is resident in memory. A disk-backed store keeps sequence bytes on disk, so n_sequences_in_memory is 0:

disk_store.stats()

{'n_sequences': '3',
 'logical_sequence_bytes': '11',
 'storage_mode': 'Encoded',
 'n_sequences_in_memory': '0',
 'n_collections_in_memory': '2',
 'n_collections': '2'}

Part 2: Remote RefgetStore with Lazy Loading

The open_remote() method allows you to work with RefgetStores hosted on remote servers without downloading everything upfront. Sequence data is fetched on demand, and can be cached locally.

2.1 Understanding Remote Stores

Key Concepts:

• Metadata loaded first: Small index files (rgstore.json, collections.rgci, and sequences.rgsi when first needed) are fetched from the remote
• Sequence data loaded on demand: get_substring() uses HTTP range requests to read only the bytes covering the region you ask for
• Local caching: load_sequence() downloads a whole sequence and caches it on disk, for repeated access
• User-controlled cache: You specify where cached data is stored

2.2 Simulating a Remote Store

For this example, we'll serve the store we created earlier over HTTP on localhost. In practice, you'd use an https:// URL pointing to an actual remote server (or an s3:// bucket).

Because get_substring() relies on HTTP range requests, the server must support them. Python's built-in http.server does not, so we start a tiny static server that does:

import socket
import subprocess
import sys
import time

# A minimal static file server that honors `Range: bytes=start-end` requests
server_code = '''
import os, re, sys
from http.server import SimpleHTTPRequestHandler, ThreadingHTTPServer

class RangeHandler(SimpleHTTPRequestHandler):
    def send_head(self):
        m = re.fullmatch(r"bytes=(\\d+)-(\\d*)", self.headers.get("Range", ""))
        path = self.translate_path(self.path)
        if not m or not os.path.isfile(path):
            return super().send_head()
        size = os.path.getsize(path)
        start = int(m.group(1))
        end = min(int(m.group(2) or size - 1), size - 1)
        f = open(path, "rb")
        f.seek(start)
        self.send_response(206)
        self.send_header("Content-Range", f"bytes {start}-{end}/{size}")
        self.send_header("Content-Length", str(end - start + 1))
        self.end_headers()
        self.remaining = end - start + 1
        return f

    def copyfile(self, source, outputfile):
        if hasattr(self, "remaining"):
            outputfile.write(source.read(self.remaining))
        else:
            super().copyfile(source, outputfile)

os.chdir(sys.argv[2])
ThreadingHTTPServer(("127.0.0.1", int(sys.argv[1])), RangeHandler).serve_forever()
'''
with open("range_server.py", "w") as f:
    f.write(server_code)

# Pick a free port and serve the saved store from it
with socket.socket() as s:
    s.bind(("127.0.0.1", 0))
    port = s.getsockname()[1]

server = subprocess.Popen(
    [sys.executable, "range_server.py", str(port), saved_store_path],
    stdout=subprocess.DEVNULL,
    stderr=subprocess.DEVNULL,
)
time.sleep(1)

# Simulated remote URL (in practice, this would be https://...)
remote_url = f"http://127.0.0.1:{port}"
print(f"Serving {saved_store_path} over HTTP")

# Create a cache directory for remote data
cache_dir = "refget_cache"
print(f"Cache directory: {cache_dir}")

Serving my_refget_store over HTTP
Cache directory: refget_cache

2.3 Open Remote Store

Open a store from a remote URL with local caching:

def show_cache():
    for root, dirs, files in os.walk(cache_dir):
        for file in sorted(files):
            full_path = os.path.join(root, file)
            rel_path = os.path.relpath(full_path, cache_dir)
            print(f"  - {rel_path} ({os.path.getsize(full_path)} bytes)")

# Open remote store
# In a real scenario, remote_url would be something like:
# "https://s3.amazonaws.com/genomes/hg38" or "https://refget.example.com/store"
remote_store = RefgetStore.open_remote(cache_dir, remote_url)

print(f"Remote store opened: {remote_store!r}")
print(f"\nCache directory contents (metadata only at this point):")
show_cache()

Remote store opened: RefgetStore(n_sequences=0, n_sequences_in_memory=0, n_collections=1, n_collections_in_memory=0, mode=Encoded, persist=on, quiet=off, cache='refget_cache', remote='http://127.0.0.1:36637')

Cache directory contents (metadata only at this point):
  - collections.rgci (379 bytes)
  - rgstore.json (420 bytes)

2.4 Access Sequences (Triggers Lazy Loading)

When you request a substring, only the bytes covering that region are fetched. Nothing is written to the sequence cache:

print("Substring access (reads just the needed bytes from the remote):")
seq1 = remote_store.get_substring(seq_digest_chr1, 0, 10)
print(f"  Sequence: {seq1}")

print(f"\nCache directory contents (after substring access):")
show_cache()

Substring access (reads just the needed bytes from the remote):
  Sequence: ATGCATGCAT

Cache directory contents (after substring access):
  - collections.rgci (379 bytes)
  - rgstore.json (420 bytes)
  - sequences.rgsi (214 bytes)


Downloading sequence index sequences.rgsi...

If you will read from the same sequence many times, download it once with load_sequence(). It is cached on disk, and later reads are served locally:

# Download the whole sequence and cache it
remote_store.load_sequence(seq_digest_chr1)

print(f"Cache directory contents (after load_sequence):")
show_cache()

# Later access uses the cached data (no network request)
print("\nSecond access (uses cache, no network request):")
seq2 = remote_store.get_substring(seq_digest_chr1, 5, 15)
print(f"  Sequence: {seq2}")

Cache directory contents (after load_sequence):
  - collections.rgci (379 bytes)
  - rgstore.json (420 bytes)
  - sequences.rgsi (214 bytes)
  - sequences/k4/k4ZQ03Uo81xEjENxmvAkcLnTez9OFlQE.seq (6 bytes)

Second access (uses cache, no network request):
  Sequence: TGCATGCAGT


Downloading sequence k4ZQ03Uo81xEjENxmvAkcLnTez9OFlQE...

Collections are lazy too. Before looking up a sequence by name in a remote store, load its collection. The lookup returns a record with metadata only; use its digest to fetch the bases:

remote_store.load_collection(collection_digest)
record = remote_store.get_sequence_by_name(collection_digest, "chr2")
print(f"Found {record.metadata.name}: {record.metadata.length} bp, loaded: {record.is_loaded}")

# The record is a lightweight stub; fetch the bases by digest
chr2_seq = remote_store.get_substring(record.metadata.sha512t24u, 0, record.metadata.length)
print(f"chr2 from remote store: {chr2_seq}")

Found chr2: 8 bp, loaded: False
chr2 from remote store: GGGGAAAA


Downloading collection metadata aep4zTAZpy0wOIfWEQ5o8oykl9M47YVJ...

2.5 Real-World Example Template

Here's how you would use this with an actual remote RefgetStore:

# Example: Opening a remote genome (pseudocode - URL would need to be real)
'''
# Open human genome hg38 from a remote server
cache_path = "/data/refget_cache/hg38"
remote_url = "https://refget-server.example.com/hg38"

hg38_store = RefgetStore.open_remote(cache_path, remote_url)

# Get a region of chr1 (only the bytes for this region are fetched)
chr1_digest = "..."
sequence = hg38_store.get_substring(chr1_digest, 1000000, 1001000)

# For many reads on chr1, cache the whole sequence first
hg38_store.load_sequence(chr1_digest)
another_region = hg38_store.get_substring(chr1_digest, 2000000, 2001000)
'''
print("See code cell for real-world usage example")

See code cell for real-world usage example

Part 3: Comparison - Local vs Remote Loading

When to Use Each Method

Method	Use Case	Pros	Cons
open_local(path)	Store is already on disk	Fast, no network	Requires full download
open_remote(cache, url)	Store hosted remotely	Only download what you need	Slower first access

Performance Characteristics

Local Store:

• Initial load: Fast (metadata only)
• Sequence access: Fast (reads only the needed bytes from local disk)

Remote Store:

• Initial load: A few small HTTP requests for the index files
• Substring access: 1 HTTP range request for just the bytes in the region
• load_sequence(): 1 HTTP request + cache write; later access is fast (cached on disk)

Benefits of Lazy Loading

1. Memory efficient: Only loaded sequences consume memory
2. Bandwidth efficient: Only download what you use
3. Fast startup: No waiting for full genome download
4. Selective caching: Control which sequences to cache locally

Cleanup

Stop the HTTP server and remove temporary files:

import shutil

# Stop the local HTTP server
server.terminate()
server.wait()

# Clean up temporary directory
os.chdir(original_dir)
shutil.rmtree(temp_dir)
print("Cleaned up temporary directory")

Cleaned up temporary directory

Summary

This notebook demonstrated:

1. ✅ Creating and populating a RefgetStore from FASTA files
2. ✅ Retrieving sequences by digest and by collection + name
3. ✅ Extracting substrings and BED file regions
4. ✅ Exporting complete FASTA files from sequence collections
5. ✅ Saving stores to disk and opening them again
6. ✅ Importing many FASTA files into a disk-backed store
7. ✅ Opening remote RefgetStores with lazy loading
8. ✅ On-demand sequence fetching with local caching

For more information, see the RefgetStore Python reference and the gtars repository.