sequence-io/format-conversion/SKILL.md
Convert between sequence file formats (FASTA, FASTQ, GenBank, EMBL, Stockholm) and re-encode FASTQ quality offsets using Biopython Bio.SeqIO. Use when changing a file format for a downstream tool, fixing FASTQ quality encoding (Phred+33 vs Phred+64 vs Solexa), or when a conversion risks silently dropping annotations or quality scores.
npx skillsauth add GPTomics/bioSkills bio-format-conversionInstall this skill globally with one command. Works with Claude Code, Cursor, and Windsurf.
3 of 9 scanners reported clean
Some scanners were skipped, did not run, or reported a non-clean status. Review each row below.
Reference examples tested with: BioPython 1.83+
Before using code patterns, verify installed versions match. If versions differ:
pip show <package> then help(module.function) to check signaturesIf code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.
"Convert this file to a different format" -> Read records in one format, optionally add or drop annotations, and write in the target format.
SeqIO.convert() for a streaming one-shot conversion, or SeqIO.parse() + SeqIO.write() when records need modification (BioPython)seqkit seq (SeqKit) for FASTA/FASTQ; samtools view for SAM/BAM/CRAMA conversion is lossy whenever the target format cannot represent the source's information. The conversion still succeeds with no error and no warning. FASTA stores only id + description + sequence, so it is the most lossy common target: converting GenBank, EMBL, or FASTQ to FASTA silently discards everything the richer format carried. Before converting, decide whether the destination can hold what the source contains. If it cannot, treat the conversion as a deliberate downgrade, not a neutral reformat.
SeqIO.convert('in.gb', 'genbank', 'out.fasta', 'fasta') discards all features, annotations, qualifiers, and dbxrefs. The genes, CDS coordinates, /product and /gene qualifiers, organism, taxonomy, references, and molecule_type are all gone. There is no error, no warning, and the record count is unchanged, so the loss is invisible unless the output is inspected. FASTA encodes only record.id, record.description, and record.seq; everything in record.features, record.annotations, and record.dbxrefs has nowhere to go.
If the features matter, do not convert to FASTA. Extract feature sequences first (see sequence-manipulation/sequence-slicing) or keep the GenBank file as the source of record and use the FASTA only as a sequence-only derivative for tools that demand FASTA.
| From | To | What is lost (silently) |
|------|-----|-------------------------|
| GenBank / EMBL | FASTA | All features, qualifiers, annotations, dbxrefs; keeps id + description + seq |
| GenBank | EMBL (or reverse) | Usually lossless; both hold features and annotations |
| FASTQ | FASTA | Per-base quality scores (phred_quality) |
| FASTQ Phred+64 | FASTQ Phred+33 | Nothing if offsets handled correctly; corruption if the wrong parser is used |
| FASTQ Phred | FASTQ Solexa | Precision at low quality (round-trip lossy below ~Q10); warns when max Solexa exceeded |
| Stockholm | FASTA | Alignment columns (gaps), consensus, per-column annotation; keeps ungapped seqs |
| Any rich format | FASTA | Everything except id + description + seq |
The general pattern: rich -> flat loses the richness. The conversion succeeds regardless.
For a plain conversion with no record modification, use SeqIO.convert(). It streams one record at a time from input to output (memory-efficient, never loads the whole file) and is preferred over parse() + write(), which is only needed when records must be changed en route.
from Bio import SeqIO
count = SeqIO.convert('input.gb', 'genbank', 'output.fasta', 'fasta')
print(f'Converted {count} records')
Parameters: in_file, in_format, out_file, out_format (filenames or handles; format strings are lowercase). Returns the number of records written. Reach for parse() + write() only when injecting annotations, transforming sequences, or filtering during the conversion.
FASTQ quality is one ASCII character per base, but the offset and score type differ across instrument generations. Re-encoding between them is a conversion, not a copy: the bytes in the quality line change.
| Format string | For | Offset | Score type |
|---------------|-----|--------|-----------|
| fastq (alias of fastq-sanger) | Sanger and modern Illumina 1.8+ | 33 | Phred 0-93 |
| fastq-sanger | same as above | 33 | Phred 0-93 |
| fastq-illumina | Illumina 1.3-1.7 | 64 | Phred 0-62 |
| fastq-solexa | pre-1.3 Solexa | 64 | Solexa odds -5..62 |
Re-encode old Illumina 1.3+ (Phred+64) to modern Sanger (Phred+33) by naming both variants. SeqIO.convert() reads with the input offset and writes with the output offset:
from Bio import SeqIO
SeqIO.convert('illumina13.fastq', 'fastq-illumina', 'sanger.fastq', 'fastq-sanger')
Never re-encode without a verified source encoding. Quality encoding cannot be auto-detected in general: ASCII >= 64 is legal in every variant, so a high-quality Sanger file and a low-quality Illumina-1.3 file can be byte-identical in their quality lines. Two failure modes follow from guessing wrong:
ValueError noting the quality string is not in the correct range for the chosen QualityIO parser.Confirm the encoding from the sequencing pipeline (or FastQC's inferred encoding) before re-encoding; do not let the agent guess the offset.
Solexa is doubly lossy. Solexa uses an odds score, Q = -10 log10(P/(1-P)), not Phred's Q = -10 log10(P), which is why Solexa scores go negative. Phred <-> Solexa conversions round a float to one ASCII char per base, so the round trip is many-to-one and lossy below ~Q10 (for example Solexa 9 and 10 both map to Phred 10). Writing fastq-solexa from a Phred-only record forces an on-the-fly lossy conversion and emits a BiopythonWarning when max(qualities) >= 62.5. There is no clean Phred -> Solexa path that avoids the loss; only re-encode toward Solexa when a legacy tool truly requires it.
FASTA has no molecule_type and no quality, so converting FASTA up to a richer format means supplying what FASTA lacked. Stream records through a generator that injects the missing field.
Goal: Convert FASTA to GenBank, which the writer refuses to produce without molecule_type.
Approach: Stream records through a generator that sets record.annotations['molecule_type'], then write as GenBank.
Reference (BioPython 1.83+):
from Bio import SeqIO
def add_molecule_type(records, mol_type='DNA'):
for record in records:
record.annotations['molecule_type'] = mol_type
yield record
records = SeqIO.parse('input.fasta', 'fasta')
SeqIO.write(add_molecule_type(records), 'output.gb', 'genbank')
Goal: Convert FASTA to FASTQ by assigning placeholder per-base quality.
Approach: Stream records through a generator that adds a phred_quality list of the right length, then write as FASTQ.
Reference (BioPython 1.83+):
from Bio import SeqIO
def add_quality(records, quality=40):
for record in records:
record.letter_annotations['phred_quality'] = [quality] * len(record.seq)
yield record
records = SeqIO.parse('input.fasta', 'fasta')
SeqIO.write(add_quality(records), 'output.fastq', 'fastq')
Placeholder quality is fabricated, not measured: downstream QC and variant callers will treat it as real. Use it only to satisfy a tool's format requirement, never to imply the bases were measured at that quality. letter_annotations is length-locked to the sequence, so the list length must equal len(record.seq).
Goal: Convert every file of one format in a directory to another format.
Approach: Glob the input files, apply SeqIO.convert() to each, and report per-file counts.
Reference (BioPython 1.83+):
from pathlib import Path
from Bio import SeqIO
for gb_file in Path('.').glob('*.gb'):
fasta_file = gb_file.with_suffix('.fasta')
count = SeqIO.convert(str(gb_file), 'genbank', str(fasta_file), 'fasta')
print(f'{gb_file.name}: {count} records')
When the conversion must also transform sequences, parse and write explicitly rather than using convert().
from Bio import SeqIO
from Bio.SeqRecord import SeqRecord
def uppercase_record(rec):
return SeqRecord(rec.seq.upper(), id=rec.id, description=rec.description)
records = SeqIO.parse('input.fasta', 'fasta')
SeqIO.write((uppercase_record(rec) for rec in records), 'output.fasta', 'fasta')
Seq is case-preserving, so lowercase soft-masking survives a plain conversion; call .upper() explicitly only when the destination tool requires uppercase.
Sequence formats drop gaps and alignment columns. To convert between alignment formats (Stockholm, PHYLIP, Clustal, FASTA-alignment) keeping the columns, use AlignIO, not SeqIO.
from Bio import AlignIO
AlignIO.convert('alignment.sto', 'stockholm', 'alignment.phy', 'phylip')
| Symptom | Cause | Fix |
|---------|-------|-----|
| GenBank features missing after conversion | Target was FASTA, which cannot hold features | Expected and silent; keep the GenBank as source, or extract features before converting |
| ValueError about missing molecule_type | Writing GenBank/EMBL from records that lack it (e.g. from FASTA) | Set record.annotations['molecule_type'] before writing |
| ValueError about quality scores | Writing FASTQ from records with no phred_quality | Add phred_quality to letter_annotations (length must equal the sequence) |
| ValueError mentioning the QualityIO parser | A quality char is outside the named parser's range (wrong FASTQ variant) | Use the correct variant: fastq-sanger, fastq-illumina, or fastq-solexa |
| FASTQ scores all off by 31 with no error | Read Phred+33 as fastq-illumina or Phred+64 as fastq-sanger (overlap region) | Confirm the true encoding from the pipeline; re-read with the right variant |
| BiopythonWarning "Data loss - max Solexa quality" | Writing fastq-solexa from Phred scores above ~62 | Expected lossy conversion; only write Solexa when a legacy tool requires it |
| Alignment columns/gaps lost | Used SeqIO on an alignment | Use AlignIO.convert() to preserve columns |
tools
End-to-end CLIP-seq pipeline from FASTQ to ENCODE-compliant binding sites, single-nucleotide crosslink maps, annotation, motifs, and (optionally) differential binding. Use when running the full Yeo lab eCLIP / iCLIP / iCLIP2 / iCLIP3 / irCLIP / PAR-CLIP analysis with SMInput control, protocol-specific UMI extraction, ENCODE STAR parameters, CLIPper or Skipper peak calling with stringent log2 FC and -log10 p thresholds, IDR rescue and self-consistency QC, and downstream motif registration with mCross or PEKA.
development
Detect, date, and contextualize whole-genome duplication (WGD / paleopolyploidy) events using wgd v2 (Chen et al 2024), KsRates (Sensalari 2022 substitution-rate-corrected Ks dating), DupGen_finder (Qiao 2019), MAPS (Li 2018 phylogenomic), POInT (Conant 2008 ordered-block), SLEDGe (2024 ML-based), Whale.jl (Bayesian DL+WGD), and synteny-anchored paranome construction. Use when identifying ancient polyploidy from Ks distributions and synteny block analysis, positioning WGD events relative to speciation, distinguishing tandem from segmental from WGD duplications, dating the 2R/3R vertebrate / fish / salmonid WGDs, building paranome and Ks-age mixture models, applying KsRates substitution-rate correction across lineages, or testing alternative biased-fractionation / dosage-balance models post-WGD.
tools
Build whole-genome alignments using Progressive Cactus (Armstrong 2020 reference-free clade-level WGA), Minigraph-Cactus (Hickey 2024 pangenome-aware), LASTZ chain/net (UCSC pipeline), MUMmer4 (Marçais 2018 pairwise), minimap2 -x asm5/10/20 (Li 2018 fast pairwise), AnchorWave (Song 2022 WGD-aware), and Mauve / progressiveMauve (bacterial). Operates the HAL toolkit (Hickey 2013) for downstream extraction including halSynteny, halLiftover, halBranchMutations, and hal2maf. Use when constructing multi-species alignments for comparative-annotation projection (TOGA), synteny detection, conservation analyses (phyloP / PhastCons), or pangenome graph construction; selecting between reference-free (Cactus) and reference-anchored (LASTZ chains/nets) approaches; tuning sensitivity for closely vs distantly related genomes; or producing HAL files for genome-wide downstream tools.
development
Detect syntenic blocks and structural rearrangements between genomes using MCScanX (Wang 2012), JCVI/MCScan (Tang 2008 Python), GENESPACE (Lovell 2022) for orthology-anchored riparian visualization, SyRI for structural variation, AnchorWave for sequence-level synteny, i-ADHoRe 3.0 for highly diverged species, SynNet for synteny networks, and ntSynt for multi-genome macrosynteny. Use when identifying collinear gene blocks across species, distinguishing macrosynteny from microsynteny, detecting inversions/translocations/duplications, anchoring orthology in WGD lineages, producing publication riparian plots, computing synteny block age via Ks (cross-references whole-genome-duplication), or running synteny-aware ortholog inference in polyploids.