Batch Sequence File Processing

Hard rules

1. No fabricated citations. Every cited work must resolve to a verifiable
2. No claim without provenance. Every quantitative or factual claim
3. No silent failure. Every script invocation, API call, or tool use must declare its exit status and what to do on non-zero. A skill that silently swallows errors is a violation.

When to use

• You have a directory tree of FASTA/FASTQ/GenBank files and need to merge them, split them, convert formats, or compute aggregate stats.
• You need to organize records by some attribute (sample ID, GC content, taxonomy prefix) into per-bin output files.
• You're prototyping before porting to nf-core / Snakemake and want a clear Python baseline.

When NOT to use

• Real production pipelines at petabyte scale — use nf-core/rnaseq, nf-core/sarek, or Snakemake with pysam-backed streaming.
• Variant-level work (use the variant-calling skills in this category instead).
• Tasks that are fundamentally one-file — see the individual read-sequences/write-sequences skills.

Prerequisites

• Python 3.10+
• biopython>=1.83 (pip install biopython)
• For Parquet/Polars summaries: polars>=1.0

Core workflow

1. Inventory first. Use Path.glob() (recursive when needed) to enumerate inputs. Always sort the result — glob order is filesystem-dependent.
2. Stream, don't slurp. Wrap file iteration in a generator that yields SeqRecord objects; never build a list of millions of records in memory.
3. Compute on the way through. Aggregate counters and statistics in a single pass.
4. Write once. SeqIO.write accepts an iterator, so the read pipeline can feed the write pipeline directly.
5. Persist summaries as Parquet/CSV for downstream analysis (Polars, pandas, DuckDB).

Code patterns

Stream and count across a directory

from pathlib import Path
from Bio import SeqIO

def count_records(directory: Path, pattern: str, fmt: str):
    for fp in sorted(directory.glob(pattern)):
        n = sum(1 for _ in SeqIO.parse(fp, fmt))"
        yield {"file": fp.name, "records": n, "format": fmt}

for row in count_records(Path("data/"), "*.fasta.gz", "fasta"):
    print(row)

Merge with provenance in the description

When the source file matters later (downstream filters, audits), tag each record:

from pathlib import Path
from Bio import SeqIO

def with_source(directory: Path, pattern: str, fmt: str):
    for fp in sorted(directory.glob(pattern)):
        for rec in SeqIO.parse(fp, fmt):
            rec.description = f"{rec.description} [source={fp.name}]"
            yield rec

n = SeqIO.write(with_source(Path("data/"), "*.fa", "fasta"),
                "merged.fasta", "fasta")
print(f"Wrote {n} records")

Split a huge FASTA into N-record chunks

from itertools import islice
from Bio import SeqIO

def split_by_count(in_path: str, fmt: str, n_per_file: int, prefix: str):
    it = SeqIO.parse(in_path, fmt)
    i = 0
    while True:
        batch = list(islice(it, n_per_file))
        if not batch:
            return
        i += 1
        out = f"{prefix}_{i:04d}.{fmt}"
        SeqIO.write(batch, out, fmt)
        print(f"{out}: {len(batch)} records")

split_by_count("huge.fasta", "fasta", 10_000, "chunk")

Split by sample/chromosome prefix

from collections import defaultdict
from Bio import SeqIO

groups = defaultdict(list)
for rec in SeqIO.parse("input.fasta", "fasta"):
    groups[rec.id.split("|")[0]].append(rec)

for prefix, recs in groups.items():
    SeqIO.write(recs, f"{prefix}.fasta", "fasta")

Batch format convert

from pathlib import Path
for gb in sorted(Path("genbank").glob("*.gb")):
    out = Path("fasta") / gb.with_suffix(".fasta").name
    n = SeqIO.convert(str(gb), "genbank", str(out), "fasta")
    print(f"{gb.name} -> {out.name}: {n}")

Aggregate per-file stats to Parquet

from pathlib import Path
import polars as pl
from Bio import SeqIO
from Bio.SeqUtils import gc_fraction

rows = []
for fa in sorted(Path("data/").glob("*.fasta")):
    lengths = [len(r.seq) for r in SeqIO.parse(fa, "fasta")]
    if not lengths:
        continue
    rows.append({
        "file": fa.name,
        "n_records": len(lengths),
        "total_bp": sum(lengths),
        "min_len": min(lengths),
        "median_len": sorted(lengths)[len(lengths) // 2],
        "max_len": max(lengths),
    })

pl.DataFrame(rows).write_parquet("summary.parquet")

Organize by GC content bin

from pathlib import Path
from Bio import SeqIO
from Bio.SeqUtils import gc_fraction

high, low = Path("high_gc"), Path("low_gc")
high.mkdir(exist_ok=True)
low.mkdir(exist_ok=True)

for fa in Path("input").glob("*.fasta"):
    recs = list(SeqIO.parse(fa, "fasta"))
    avg_gc = sum(gc_fraction(r.seq) for r in recs) / len(recs)
    dest = high if avg_gc >= 0.55 else low
    SeqIO.write(recs, dest / fa.name, "fasta")

Parallel count with concurrent.futures

ThreadPoolExecutor is usually fine for I/O-bound sequence parsing; switch to ProcessPoolExecutor if you do CPU-heavy work per record (translation, alignment, etc.).

from concurrent.futures import ThreadPoolExecutor
from pathlib import Path
from Bio import SeqIO

def count(fp: Path) -> tuple[str, int]:
    return fp.name, sum(1 for _ in SeqIO.parse(fp, "fastq"))

with ThreadPoolExecutor(max_workers=8) as ex:
    for name, n in ex.map(count, sorted(Path(".").glob("*.fastq.gz"))):
        print(f"{name}\t{n}")

Common pitfalls

• Memory blow-up on list(SeqIO.parse(...)). If the file is >1 GB, materialize only the records you need; otherwise stream.
• Non-deterministic order. Path.glob() results are not guaranteed sorted. Always sorted().
• Hidden '.fasta.gz' formats. SeqIO reads .gz natively, but you must pass "fasta" (not "fasta-gz") and accept the slower text path.
• Empty summary.csv writes. Don't write a header if summaries is empty — handle the zero-record case.
• Mixing GenBank and FASTA with the same parser. Use SeqIO.convert(src, src_fmt, dst, dst_fmt) for format-boundary code; it returns record count for free.

Validation

• Run a row-count check: total records out equals total records in (modulo your filter).
• For round-trip conversions (GenBank → FASTA → GenBank), verify the FASTA round-trips back to a valid GenBank.
• Spot-check: head -n 4 merged.fasta | grep '^>' should show the [source=...] provenance tags.

Open alternatives

• nf-core modules for production: nf-core/modules has a samtools_faidx and seqkit_stats you can call directly.
• SeqKit (seqkit stats *.fasta) is faster than Python for stats-only tasks and is installed in most BioContainers images.
• SeqKit + GNU parallel if you want shell-only batch processing: seqkit stats -j 8 *.fasta | tee summary.tsv.

References

• Biopython SeqIO tutorial: https://biopython.org/wiki/SeqIO
• Polars IO docs: https://pola-rs.github.io/polars/py-polars/html/reference/io.html
• Companion: ors-bioinformatics-sequence-compressed-sequence-files for .gz/.bgz handling.

Changelog

• 1.0.0 (2026-06-10): Initial adaptation by Pradyumna Jayaram from bio-batch-processing (bioSkills-main/sequence-io/batch-processing).

Cross-references

Other skills in this category:

• bowtie2-alignment
• bwa-alignment
• bwa-mem2-alignment
• codon-usage
• compressed-sequence-files
• fastq-quality-scores
• filter-sequences
• format-conversion
• hisat2-alignment
• motif-search
• paired-end-fastq
• pysam-genomics
• read-write-sequences
• reverse-complement
• sam-bam-basics
• samtools-bam-processing
• seq-objects
• sequence-properties
• sequence-slicing
• sequence-statistics
• star-alignment
• transcription-translation