bio-write-sequences

name	bio-write-sequences
description	Write biological sequences to files (FASTA, FASTQ, GenBank, EMBL) using Biopython Bio.SeqIO. Use when saving sequences, creating new sequence files, or outputting modified records.
origin	openai4s
category	bioskills/sequence-io
metadata	{"tool_type":"python","primary_tool":"Bio.SeqIO","third_party":{"name":"GPTomics/bioSkills","repository":"https://github.com/GPTomics/bioSkills","commit":"d91ed3d563019e649dc854c56ccd62551359488a","license":"MIT"}}

Version Compatibility

Reference examples tested with: BioPython 1.83+

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

• Python: pip show <package> 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.

Write Sequences

"Write sequences to a file" -> Serialize SeqRecord objects into a formatted sequence file.

• Python: SeqIO.write() (BioPython)
• R: writeXStringSet() (Biostrings)

Governing Principle

The write is only as complete as the SeqRecord. Each format reads specific record fields and silently ignores the rest, so what survives a write is decided by which fields are populated before the call, not by the format string. FASTA serializes only id/description+seq; FASTQ additionally requires letter_annotations['phred_quality']; GenBank/EMBL additionally require annotations['molecule_type']. Populate the fields a format needs, or the write either drops data quietly (FASTA) or raises (FASTQ/GenBank).

Required Import

from Bio import SeqIO
from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord

Core Functions

SeqIO.write() - Write Records to File

SeqIO.write(records, 'output.fasta', 'fasta')

• records - Single SeqRecord, list, or iterator of SeqRecords
• handle - Filename (string) or open file handle
• format - Lowercase output format string
• Returns the number of records written (integer)

record.format() - Get Formatted String

formatted = record.format('fasta')

The FASTA Header Trap (id vs description)

FASTA output is built from record.description, NOT record.id. The writer compares the first whitespace token of description to id: if they match it writes description as-is; otherwise it prepends id + a space. When a record is parsed from FASTA, description already leads with the id token, so a round trip is faithful. But when id and description are set independently, a stale id-like token inside the description gets duplicated, and older BioPython releases dropped the id entirely instead of prepending it.

record.id	record.description	Header written
seq1	seq1 kinase domain	>seq1 kinase domain (clean: description leads with id)
seq1	kinase domain	>seq1 kinase domain (id auto-prepended)
seq1	`` (empty)	>seq1 (id used as fallback)
seq1	gene7 kinase domain	>seq1 gene7 kinase domain (stale id duplicated)

To control the header exactly and stay robust across versions, make description begin with id + a space: description=f'{rec_id} kinase domain'. The FASTA writer wraps the sequence at 60 characters per line by default.

Format Field Requirements

Format	String	Record fields read	Hard requirement
FASTA	'fasta'	id/description, seq	none (header trap above)
FASTQ	'fastq'	seq, letter_annotations	phred quality scores
GenBank	'genbank' / 'gb'	seq, annotations, features	molecule_type
EMBL	'embl'	seq, annotations, features	molecule_type
Tab	'tab'	id, seq	none

Creating SeqRecord Objects

Goal: Construct in-memory records that carry the fields the target format requires.

Approach: Build a SeqRecord from a Seq plus id; add letter_annotations['phred_quality'] for FASTQ and annotations['molecule_type'] for GenBank/EMBL.

"Create a sequence record from scratch" -> Wrap a Seq in a SeqRecord with metadata.

• Python: SeqRecord(Seq(...), id=...) (BioPython)

record = SeqRecord(Seq('ATGCGATCGATCG'), id='seq1', description='seq1 example sequence')

Code Patterns

Write Single or Multiple Records

records = [SeqRecord(Seq('ATGC'), id='seq1'), SeqRecord(Seq('GCTA'), id='seq2')]
count = SeqIO.write(records, 'output.fasta', 'fasta')

Write to a File Handle (and Append)

with open('output.fasta', 'w') as handle:
    SeqIO.write(records, handle, 'fasta')

with open('output.fasta', 'a') as handle:
    SeqIO.write(new_records, handle, 'fasta')

Write Modified Records via Generator

Goal: Transform sequences in memory and write the modified versions to a new file.

Approach: Parse input, map a transform over a generator, write the generator. Streaming avoids loading every record into RAM.

"Modify sequences and save" -> Parse records, transform each, write with SeqIO.write().

def uppercase_record(rec):
    return SeqRecord(rec.seq.upper(), id=rec.id, description=rec.description)

records = SeqIO.parse('input.fasta', 'fasta')
modified = (uppercase_record(rec) for rec in records)
SeqIO.write(modified, 'output.fasta', 'fasta')

Write FASTQ with Quality Scores

FASTQ requires letter_annotations['phred_quality'] as a list of ints. letter_annotations is length-locked to len(seq): assigning a list whose length differs from the sequence raises. Set the sequence first, then the quality list of matching length.

record = SeqRecord(Seq('ATGCGATCG'), id='read1')
record.letter_annotations['phred_quality'] = [40] * len(record.seq)
SeqIO.write(record, 'output.fastq', 'fastq')

Quality Encoding on Write (Phred vs Solexa)

When both phred_quality and solexa_quality keys are present, the writer uses Phred. Writing 'fastq-solexa' from a Phred-only record forces an on-the-fly lossy conversion (the scales diverge in the low-quality region) and emits a BiopythonWarning once any score reaches the high end (max quality >= ~62). For modern data, write plain 'fastq' (Sanger/Phred+33); only use 'fastq-solexa'/'fastq-illumina' when a tool explicitly demands that legacy encoding.

Write GenBank Format

GenBank and EMBL writing requires annotations['molecule_type'] (the alphabet that once carried this was removed in BioPython 1.78). Missing it raises on write.

record = SeqRecord(Seq('ATGCGATCGATCG'), id='SEQ001', name='example')
record.annotations['molecule_type'] = 'DNA'
record.annotations['topology'] = 'linear'
record.annotations['organism'] = 'Example organism'
SeqIO.write(record, 'output.gb', 'genbank')

Common Errors

Symptom	Cause	Fix
Header has a duplicated or mangled id	FASTA builds the header from description; it does not lead with id + space	Set description=f'{rec.id} ...' or leave description empty to fall back to id
ValueError: No suitable quality scores found in letter_annotations of SeqRecord (id=...) on FASTQ write	Record has no letter_annotations['phred_quality']	Assign record.letter_annotations['phred_quality'] = [q]*len(seq)
TypeError: Any per-letter annotation should be a Python sequence ... of the same length	Quality list length != len(seq) (annotations are length-locked)	Set seq first, then a quality list of matching length
ValueError: missing molecule_type ... on GenBank/EMBL write	No annotations['molecule_type'] since the 1.78 alphabet removal	Add record.annotations['molecule_type'] = 'DNA' (or 'RNA'/'protein')
BiopythonWarning: Data loss - max Solexa quality ...	Writing 'fastq-solexa' from a high Phred-only record forces lossy conversion	Write plain 'fastq' unless a tool requires the Solexa encoding
TypeError passing a raw str/Seq to write	SeqIO.write expects SeqRecord(s)	Wrap the sequence in a SeqRecord first
ValueError: Sequences must all be the same length	PHYLIP/alignment format with unequal lengths	Align, pad, or trim to equal length first

Related Skills

• read-sequences - Read sequences before modifying and writing
• format-conversion - Direct format conversion without intermediate processing
• filter-sequences - Filter sequences before writing a subset
• fastq-quality - Phred/Solexa encodings and quality-score handling
• sequence-manipulation/seq-objects - Create SeqRecord objects to write
• alignment-files/sam-bam-basics - For SAM/BAM output, use samtools/pysam