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

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

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

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)

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

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