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: ```python 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: ```python # 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: ```python # 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. ```python # 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()`: ```python # 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: ```python # 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: ```python # 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: ```python # 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: ```python # 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: ```python # 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: ```python # 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: ```python # 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: ```python # 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: ```python 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: ```python 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: ```python 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: ```python 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: ```python # 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: ```python 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: ```python # 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: ```python 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](/gtars/python/refget-store.md) and the [gtars repository](https://github.com/databio/gtars).