GPTomics/bio-read-sequences

Version Compatibility

Reference examples tested with: BioPython 1.83+

Before using code patterns, verify installed versions match. If versions differ:

• Python: pip show biopython then help(module.function) to check signatures

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

Read Sequences

Read biological sequence data from files using Biopython's Bio.SeqIO module.

"Read sequences from a file" -> Parse a file into SeqRecord objects exposing id, sequence, and annotations.

• Python: SeqIO.parse() / SeqIO.read() (BioPython)
• R: readDNAStringSet() / readAAStringSet() (Biostrings)

The Governing Principle

Stream by default. SeqIO.parse() yields one record at a time and never holds the whole file in RAM, so it scales to any size. Reach for an in-memory or indexed structure only when the access pattern demands it: load all records (to_dict) only for small files needing random access; build an index (index / index_db) for random access into large files. Never list() a huge file or to_dict() it - that defeats streaming and can exhaust memory.

Which Function to Use

| Method | Returns | Memory model | Random access | Persists | Multi-file | |--------|---------|--------------|---------------|----------|------------| | parse(handle, format) | generator of SeqRecord | one record at a time | no | no | no | | read(handle, format) | one SeqRecord | one record | n/a | no | no | | to_dict(records) | real dict | ALL records in RAM | yes | no | feed combined iterators | | index(filename, format) | dict-like (read-only) | byte offsets only, re-parses on access | yes | no | no | | index_db(idx_file, files, format) | dict-like (read-only) | on-disk SQLite index | yes | yes | yes |

Decision rule: parse for streaming; read for a known single-record file; to_dict when the file is small and random access by ID is needed; index for random access into one large file; index_db for files larger than RAM, many files indexed together, or an index reused across runs.

Behavioral traps these methods hide:

• parse() is a one-pass generator. It is NOT subscriptable (parse(...)[3] raises TypeError), and it EXHAUSTS SILENTLY: a second for loop over the same generator object yields nothing with no error. Re-call parse() for each pass, or list() it once if the file is small.
• read() fails LOUDLY: zero records raise ValueError: No records found in handle; more than one raises ValueError: More than one record found in handle. Use it as an assertion that the file holds exactly one sequence.
• to_dict(), index(), and index_db() all raise ValueError on a DUPLICATE id (Duplicate key '...'). Supply a key_function to derive unique keys when ids collide.
• index() needs a FILENAME, not a handle (it must seek). It stores only byte offsets and re-parses the record from disk on every access, so it returns a fresh object each time and mutations do not persist. It is read-only (__setitem__ raises NotImplementedError).
• index_db() stores the offset index in an on-disk SQLite file. It PERSISTS across sessions (reopen later with just the index filename), and scales beyond RAM and across multiple files (pass a list of filenames). This is the right answer for data larger than memory.

The alphabet= argument still appears in some signatures for back-compatibility but is a no-op since BioPython 1.78; leave it None.

Required Import

from Bio import SeqIO

Reading Records

SeqIO.parse() - Stream Multiple Records

Returns a one-pass iterator of SeqRecord objects. Always pass the format explicitly as the second argument.

for record in SeqIO.parse('sequences.fasta', 'fasta'):
    print(record.id, len(record.seq))

SeqIO.read() - Exactly One Record

Use when the file must contain a single sequence; raises on zero or multiple records.

record = SeqIO.read('single.fasta', 'fasta')

Random Access

SeqIO.to_dict() - Small Files

Loads every record into a dictionary keyed by id. Fast random access, but holds all records in RAM.

records = SeqIO.to_dict(SeqIO.parse('sequences.fasta', 'fasta'))
seq = records['sequence_id'].seq

SeqIO.index() - One Large File

Goal: Random access by id into a large file without loading every record into memory.

Approach: Build an in-memory map of byte offsets keyed by id; each lookup re-parses one record from disk.

Reference (BioPython 1.83+):

records = SeqIO.index('large.fasta', 'fasta')
seq = records['sequence_id'].seq
records.close()

A key_function maps the id STRING to a custom key (note: to_dict's key_function receives the whole record instead):

def get_accession(identifier):
    return identifier.split('.')[0]  # drop the version suffix

records = SeqIO.index('sequences.fasta', 'fasta', key_function=get_accession)

SeqIO.index_db() - Huge / Multiple Files

Goal: Random access into data larger than RAM, or across many files, with the index reusable across runs.

Approach: Persist the offset index in an on-disk SQLite database; reopen it later without re-parsing.

Reference (BioPython 1.83+):

# First call parses the file(s) and builds the SQLite index
records = SeqIO.index_db('index.sqlite', 'large.fasta', 'fasta')
seq = records['sequence_id'].seq
records.close()

# Later sessions reopen instantly with just the index filename
records = SeqIO.index_db('index.sqlite')

# Index multiple files as one database
records = SeqIO.index_db('combined.sqlite', ['file1.fasta', 'file2.fasta'], 'fasta')

High-Performance Parsing

For maximum throughput on large files, low-level parsers (SimpleFastaParser, FastqGeneralIterator) yield raw tuples and skip SeqRecord construction, so they run substantially faster than SeqIO.parse.

SimpleFastaParser

Goal: Parse large FASTA files at maximum speed without SeqRecord overhead.

Approach: Iterate (title, sequence) string tuples directly from the handle.

Reference (BioPython 1.83+):

from Bio.SeqIO.FastaIO import SimpleFastaParser

with open('large.fasta') as handle:
    for title, sequence in SimpleFastaParser(handle):
        if len(sequence) > 1000:
            seq_id = title.split()[0]  # first whitespace token is the id

FastqGeneralIterator

Goal: Parse large FASTQ files at maximum speed.

Approach: Iterate (title, sequence, quality_string) string tuples; decode quality manually if needed.

Reference (BioPython 1.83+):

from Bio.SeqIO.QualityIO import FastqGeneralIterator

with open('reads.fastq') as handle:
    for title, sequence, quality in FastqGeneralIterator(handle):
        avg_qual = sum(ord(c) - 33 for c in quality) / len(quality)  # Phred+33

SeqRecord Attributes

After parsing, each record exposes:

record.id          # first whitespace token of the header (string)
record.name        # same first token (for FASTA, name == id)
record.description # the ENTIRE header after '>', including the id token
record.seq         # sequence data (Seq object; case-preserving)
record.features    # list of SeqFeature objects (GenBank/EMBL)
record.annotations # dict of annotations (organism, molecule_type, ...)
record.letter_annotations  # per-letter dict (e.g. 'phred_quality' list)
record.dbxrefs     # database cross-references

id vs name vs description - the first-space split

A FASTA header >FIRST rest of the line parses to: id = FIRST (the first whitespace token), name = FIRST (same token), description = FIRST rest of the line (the WHOLE header after >, including the id). So >seq1 some desc gives id seq1, name seq1, description seq1 some desc. The id is therefore the leading word of the description, not a separate field - relevant when writing records back out.

Common Formats

| Format | String | Typical Extension | Notes | |--------|--------|-------------------|-------| | FASTA | 'fasta' | .fasta, .fa, .fna, .faa | Most common | | FASTA 2-line | 'fasta-2line' | .fasta | One line per sequence (no wrapping) | | FASTQ | 'fastq' | .fastq, .fq | Alias of fastq-sanger (Phred+33) | | FASTQ Solexa | 'fastq-solexa' | .fastq | Old Solexa (Solexa+64, scores -5..62) | | FASTQ Illumina | 'fastq-illumina' | .fastq | Illumina 1.3-1.7 (Phred+64) | | GenBank | 'genbank' or 'gb' | .gb, .gbk | With features/annotations | | EMBL | 'embl' | .embl | European format with features | | Swiss-Prot | 'swiss' | .dat | UniProt format |

FASTQ quality encoding cannot be auto-detected reliably: the same quality line can be valid Phred+33 and Phred+64. Picking the wrong string can silently shift every score by 31. Confirm the encoding before parsing; see fastq-quality for the full encoding decision.

Specialized Formats

| Format | String | Use Case | |--------|--------|----------| | ABI | 'abi' | Sanger sequencing trace files (.ab1) | | ABI Trimmed | 'abi-trim' | ABI with low-quality ends trimmed | | SFF | 'sff' | 454/Ion Torrent flowgram data | | SFF Trimmed | 'sff-trim' | SFF with adapter/quality trimming | | QUAL | 'qual' | Quality scores file (pairs with FASTA) | | PDB SEQRES | 'pdb-seqres' | Protein sequences from PDB SEQRES records | | PDB ATOM | 'pdb-atom' | Sequences from ATOM records in PDB | | SnapGene | 'snapgene' | SnapGene .dna files |

Reading ABI Trace Files

record = SeqIO.read('sample.ab1', 'abi')
qualities = record.letter_annotations['phred_quality']
record_trimmed = SeqIO.read('sample.ab1', 'abi-trim')  # low-quality ends removed

Reading 454/Ion Torrent SFF

for record in SeqIO.parse('reads.sff', 'sff'):
    print(record.id, len(record.seq))

Reading PDB Sequences

for record in SeqIO.parse('structure.pdb', 'pdb-seqres'):
    print(record.id, record.seq)

Alignment Formats (Read-Only)

| Format | String | Notes | |--------|--------|-------| | PHYLIP | 'phylip' | Interleaved; 'phylip-relaxed' allows longer names | | Clustal | 'clustal' | ClustalW output | | Stockholm | 'stockholm' | Rfam/Pfam alignments | | NEXUS | 'nexus' | PAUP/MrBayes format | | MAF | 'maf' | Multiple Alignment Format |

Code Patterns

Count Records Without Loading All

count = sum(1 for _ in SeqIO.parse('sequences.fasta', 'fasta'))

Read GenBank with Features

for record in SeqIO.parse('sequence.gb', 'genbank'):
    for feature in record.features:
        if feature.type == 'CDS':
            product = feature.qualifiers.get('product', ['Unknown'])[0]
            cds_seq = feature.extract(record.seq)  # spliced feature sequence

Access FASTQ Quality Scores

for record in SeqIO.parse('reads.fastq', 'fastq'):
    qualities = record.letter_annotations['phred_quality']
    avg_quality = sum(qualities) / len(qualities)

Read From a File Handle

with open('sequences.fasta') as handle:
    for record in SeqIO.parse(handle, 'fasta'):
        print(record.id)

Common Errors

| Symptom | Cause | Fix | |---------|-------|-----| | Second loop over a parser yields nothing, no error | parse() generator exhausted after the first pass | Re-call parse() per pass, or list() once for small files | | TypeError: 'generator' object is not subscriptable | Indexed/sliced a parse() result | Wrap in list(), or use to_dict/index for keyed access | | ValueError: More than one record found in handle | read() on a multi-record file | Use parse() | | ValueError: No records found in handle | read() on an empty/zero-record file | Check the file and format string; use parse() if multi-record | | ValueError: Duplicate key '...' | to_dict/index/index_db hit a repeated id | Pass a key_function that derives unique keys | | Random access by id silently slow / re-reads disk | index() re-parses each access; mutations don't persist | Expected; cache needed records, or use to_dict for small files | | MemoryError / process killed on a huge file | list() or to_dict() loaded everything into RAM | Stream with parse(); use index_db() for random access | | ValueError: unknown format | Misspelled format string | Use a lowercase string from the format tables | | ValueError/AssertionError naming the LOCUS line | GenBank parser reads fixed LOCUS columns (molecule type ~44-54, topology ~55-63); ICE/SnapGene/Ensembl/assembler LOCUS lines violate the spec | Biologically valid content can still fail the strict column parse; fix the LOCUS columns or re-export from a spec-compliant writer | | FASTQ scores all off by ~31 with no error | Wrong FASTQ variant string (Phred+33 vs +64 overlap) | Confirm encoding; see fastq-quality | | AttributeError referencing .alphabet | Code assumes pre-1.78 alphabet API | Drop alphabet usage; molecule type lives in annotations['molecule_type'] |

Related Skills

• write-sequences - Write parsed sequences to new files
• filter-sequences - Filter sequences by criteria after reading
• format-conversion - Convert between formats (GenBank->FASTA silently drops annotations)
• compressed-files - Read gzip/bzip2/BGZF compressed files; only BGZF supports indexed random access
• fastq-quality - FASTQ encoding (Phred vs Solexa) and offset selection
• sequence-manipulation/seq-objects - Work with parsed SeqRecord and Seq objects
• database-access/entrez-fetch - Fetch sequences from NCBI instead of local files
• alignment-files/sam-bam-basics - For SAM/BAM/CRAM alignment files, use samtools/pysam