Read Sequences

---
name: bio-read-sequences
description: Read biological sequence files (FASTA, FASTQ, GenBank, EMBL, ABI, SFF) with Biopython Bio.SeqIO, choosing between streaming, in-memory, and on-disk-indexed access. Use when parsing sequence files, iterating multi-record files, randomly accessing records by ID in large files, or maximizing parse throughput.
tool_type: python
primary_tool: Bio.SeqIO
---

## 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

```python
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.

```python
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.

```python
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.

```python
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+):**
```python
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):
```python
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+):**
```python
# 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+):**
```python
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+):**
```python
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:

```python
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
```python
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
```python
for record in SeqIO.parse('reads.sff', 'sff'):
    print(record.id, len(record.seq))
```

### Reading PDB Sequences
```python
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
```python
count = sum(1 for _ in SeqIO.parse('sequences.fasta', 'fasta'))
```

### Read GenBank with Features
```python
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
```python
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
```python
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