Skip to content

Tutorial

Jump to bottom
Asaf Zorea edited this page Jan 25, 2026 · 2 revisions

Quick Start & Tutorial

Basic Usage

Opening a FASTA File

from fastaccess import FastaStore

# Open a FASTA file
fa = FastaStore("genome.fa")

On first use, fastaccess builds an index and saves it as genome.fa.fidx. Subsequent opens load from cache (40x faster).

Fetching Sequences

Use 1-based inclusive coordinates (standard in bioinformatics):

# Fetch 1000 bp from chr1, positions 1000-2000 (inclusive)
seq = fa.fetch("chr1", 1000, 2000)
print(len(seq))  # 1001 bases (inclusive range)
print(seq[:50])  # First 50 bases

Note: The range is inclusive on both ends, so fetch("chr1", 1, 100) returns 100 bases.

Common Patterns

Reverse Complement

# Get reverse complement
rc = fa.fetch("chr1", 1000, 2000, reverse_complement=True)

Exploring Available Sequences

# List all sequence names
sequences = fa.list_sequences()
print(sequences)  # ["chr1", "chr2", ..., "chrX", "chrY"]

# Get sequence length
length = fa.get_length("chr1")
print(f"chr1 is {length:,} bp")  # chr1 is 248,956,422 bp

# Get full header description
desc = fa.get_description("chr1")
print(desc)  # May contain additional info from FASTA header

# Get all info at once
info = fa.get_info("chr1")
print(info)
# {'name': 'chr1', 'description': '...', 'length': 248956422}

Batch Fetching

Fetch multiple regions efficiently:

# Define queries as (name, start, stop) tuples
queries = [
    ("chr1", 1000, 2000),
    ("chr1", 5000, 6000),
    ("chr2", 100, 500),
]

# Fetch all at once
sequences = fa.fetch_many(queries)

for seq in sequences:
    print(f"Length: {len(seq)}, First 20bp: {seq[:20]}")

Working with Gzipped Files

Gzipped FASTA files are supported transparently:

# Works exactly the same
fa = FastaStore("genome.fa.gz")
seq = fa.fetch("chr1", 1000, 2000)

Note: Gzip files are fully decompressed to memory on first access (gzip doesn't support random seeking).

Cache Management

Custom Cache Directory

Useful when FASTA files are in read-only locations:

# FASTA in read-only directory, cache in writable location
fa = FastaStore(
    "/readonly/data/genome.fa",
    cache_dir="/tmp/fasta_cache"
)

Disabling Cache

# Don't save or load cache (slower on subsequent opens)
fa = FastaStore("genome.fa", use_cache=False)

Checking Cache Status

# Was this loaded from cache?
if fa.is_cached():
    print("Loaded from cache")
else:
    print("Built new index")

# Does a cache file exist?
if fa.cache_exists():
    print(f"Cache at: {fa.get_cache_path()}")

# Force rebuild
fa.rebuild_index()
print("Index rebuilt and cache updated")

# Delete cache file
if fa.delete_cache():
    print("Cache deleted")

Error Handling

try:
    seq = fa.fetch("chr99", 1, 100)
except KeyError as e:
    print(f"Sequence not found: {e}")

try:
    seq = fa.fetch("chr1", -5, 100)
except ValueError as e:
    print(f"Invalid coordinates: {e}")

try:
    seq = fa.fetch("chr1", 1, 999999999999)
except ValueError as e:
    print(f"Coordinates out of range: {e}")

Advanced Examples

Extracting Gene Regions

from fastaccess import FastaStore

fa = FastaStore("hg38.fa")

# Gene coordinates (example)
gene = {
    "name": "BRCA1",
    "chr": "chr17",
    "start": 43044295,
    "end": 43125483,
    "strand": "-"
}

# Fetch sequence
seq = fa.fetch(
    gene["chr"],
    gene["start"],
    gene["end"],
    reverse_complement=(gene["strand"] == "-")
)

print(f"{gene['name']}: {len(seq)} bp")

Sliding Window Analysis

from fastaccess import FastaStore

fa = FastaStore("genome.fa")

def gc_content(seq):
    """Calculate GC content of a sequence."""
    gc = seq.count('G') + seq.count('C')
    return gc / len(seq) if len(seq) > 0 else 0

# Analyze GC content in 1kb windows
chr_name = "chr1"
window_size = 1000
step = 500  # 50% overlap

chr_len = fa.get_length(chr_name)

for start in range(1, chr_len - window_size + 1, step):
    end = start + window_size - 1
    seq = fa.fetch(chr_name, start, end)
    gc = gc_content(seq)
    print(f"{chr_name}:{start}-{end} GC={gc:.2%}")

Extracting Multiple Genes

from fastaccess import FastaStore

fa = FastaStore("genome.fa")

genes = [
    ("chr1", 1000000, 1005000),
    ("chr1", 2000000, 2008000),
    ("chr2", 500000, 505000),
]

# Batch fetch all genes
sequences = fa.fetch_many(genes)

for (chr, start, end), seq in zip(genes, sequences):
    print(f"{chr}:{start}-{end} = {len(seq)} bp, starts with {seq[:20]}")

Processing Large Genomes Efficiently

from fastaccess import FastaStore

# Large genome (e.g., human)
fa = FastaStore("hg38.fa")  # ~3GB FASTA

# First load: ~2 seconds (builds index)
# Subsequent loads: 0.05 seconds (from cache)

# Check backend
from fastaccess import using_cpp_backend
if using_cpp_backend():
    print("Using fast C++ backend")
else:
    print("Using Python backend (consider building from source)")

# Random access is fast regardless of genome size
seq = fa.fetch("chr1", 100000000, 100001000)  # Instant

Performance Tips

  1. Build from source for C++ backend (13x faster indexing)
  2. Use cache - rebuilding index on every open is slow
  3. Batch fetches with fetch_many() when fetching multiple regions
  4. Custom cache dir for read-only FASTA locations
  5. Keep index files - they're small and speed up loading 40x

Backend Selection

from fastaccess import using_cpp_backend

if using_cpp_backend():
    print("C++ backend active - optimal performance")
else:
    print("Python backend active - works everywhere")
    print("For better performance, install from source with C++ support")

Coordinate System

FASTAccess uses 1-based inclusive coordinates, matching standard bioinformatics tools:

fa = FastaStore("genome.fa")

# First 100 bases of chr1
seq = fa.fetch("chr1", 1, 100)  # Returns 100 bases

# Positions 1000-2000 (inclusive)
seq = fa.fetch("chr1", 1000, 2000)  # Returns 1001 bases

# Single base at position 5000
seq = fa.fetch("chr1", 5000, 5000)  # Returns 1 base

This matches formats like BED (when converted), GFF, and SAMtools.

Next Steps

Clone this wiki locally