bio-batch-processing

name	bio-batch-processing
description	Process many sequence files in batch (count, merge, split, convert, summarize) with memory-safe streaming and on-disk indexing using Biopython, pysam, or pyfastx. Use when iterating over a directory of FASTA/FASTQ files, merging or splitting datasets, building random access across many or huge files, or automating per-file operations without exhausting RAM.
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+ (alternatives: pysam 0.22+, pyfastx 2.0+)

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.

Batch Processing

"Process all my sequence files in a directory" -> Iterate, merge, split, convert, and summarize across multiple sequence files without loading everything into RAM.

• Python: SeqIO.parse() + Path.glob() (BioPython, pathlib) for streaming
• Python: SeqIO.index_db() (BioPython) for persistent random access across many files
• Python: pysam.FastxFile (pysam) or pyfastx for fast iteration over huge FASTQ

The Governing Principle

list(SeqIO.parse(...)) materializes every SeqRecord in RAM at once. On a directory of large files this causes OOM. SeqIO.parse() itself returns a generator that holds one record at a time, so streaming is the default for batch work: iterate, never list(), unless the file is known-small and needs multiple passes.

For random access across many or huge files, do not load them. SeqIO.index_db() builds one on-disk SQLite index over a list of files that persists across sessions. That, not to_dict(), is the batch random-access tool.

For tens of millions of reads, SeqIO is slow by design: it constructs a full SeqRecord (a Seq, id/name/description, and a letter_annotations dict of per-base qualities) for every read. When the job is plain linear iteration, a thinner reader wins.

Choosing a Reader

Reader	Per-record object	Random access	Best for
Bio.SeqIO.parse	full SeqRecord (rich API)	no (one-pass generator)	small/medium data needing the Biopython record API
Bio.SeqIO.index_db	reparsed SeqRecord on access	yes, on-disk SQLite, multi-file, persists	batch random access across many/huge files
pysam.FastxFile	thin entry (.name/.sequence/.comment/.quality)	no (linear, gzip sequential)	fast linear iteration over huge FASTQ
pyfastx	tuple/object via SQLite index	yes, into plain or gzipped FASTA/Q	random access + indexed reuse of gzipped files

pysam.FastxFile exposes .name, .sequence, .comment, .quality, and .get_quality_array() (offset-removed int Phred, but it always subtracts 33, so it is correct only for Phred+33 data - for legacy Phred+64/Solexa stay on SeqIO with the explicit variant string). pyfastx builds a persistent .fxi/.fqi SQLite index and reads random records out of plain or gzipped files without re-bgzipping.

Required Imports

from pathlib import Path
from Bio import SeqIO

Iterate and Count Across Files

Count by iterating, never by building a list. len(list(SeqIO.parse(f))) loads the whole file; sum(1 for _ in ...) holds one record at a time.

for fasta_file in Path('data/').glob('*.fasta'):
    count = sum(1 for _ in SeqIO.parse(fasta_file, 'fasta'))
    print(f'{fasta_file.name}: {count} sequences')

Recursive search uses rglob:

for gb_file in Path('data/').rglob('*.gb'):
    print(f'Found: {gb_file}')

For huge FASTQ where only sequence content matters, skip SeqRecord construction entirely:

import pysam

with pysam.FastxFile('reads.fastq.gz') as fh:
    count = sum(1 for _ in fh)

Random Access Across Many Files

Goal: Look up records by id across a whole directory of files, repeatedly, without holding them in RAM.

Approach: Build one persistent on-disk SQLite index over the file list with index_db. Reopen later with just the index path; lookups reparse single records from disk on demand.

Reference (BioPython 1.83+):

from pathlib import Path
from Bio import SeqIO

files = [str(p) for p in Path('data/').glob('*.fasta')]
records = SeqIO.index_db('combined.idx', files, 'fasta')

print(len(records))            # total across all files
record = records['seq_00042']  # random access by id
records.close()

The index file persists. A later session calls SeqIO.index_db('combined.idx') with no file list and reopens instantly. Ids must be unique across the merged set: a collision raises ValueError: Duplicate key. index_db also indexes BGZF-compressed files; plain gzip is not seekable and cannot be indexed.

Merge Files

Goal: Concatenate sequences from many files into one output without loading them all.

Approach: Chain per-file generators with yield from and stream straight into SeqIO.write, which consumes the generator one record at a time.

Reference (BioPython 1.83+):

def all_records(directory, pattern, format):
    for filepath in Path(directory).glob(pattern):
        yield from SeqIO.parse(filepath, format)

count = SeqIO.write(all_records('data/', '*.fasta', 'fasta'), 'merged.fasta', 'fasta')
print(f'Merged {count} records')

Merge with Source Tracking

Goal: Combine sequences from multiple files, tagging each record with its source filename.

Approach: Stream records through a generator that appends source metadata to the description before writing.

Reference (BioPython 1.83+):

def records_with_source(directory, pattern, format):
    for filepath in Path(directory).glob(pattern):
        for record in SeqIO.parse(filepath, format):
            record.description = f'{record.description} [source={filepath.name}]'
            yield record

SeqIO.write(records_with_source('data/', '*.fasta', 'fasta'), 'merged_tracked.fasta', 'fasta')

When merging files that may share ids, decide upfront: write-then-merge tolerates duplicates (FASTA allows repeated ids), but any later index_db/to_dict over the merged file raises on the duplicate.

Split Files

Split by Number of Records

Goal: Divide a large file into chunks of N records each.

Approach: Consume the parse generator in fixed-size batches with islice, writing each batch to a numbered file. islice pulls only N records into memory per chunk, so an arbitrarily large input streams safely.

Reference (BioPython 1.83+):

from itertools import islice

def split_file(input_file, format, records_per_file, output_prefix):
    records = SeqIO.parse(input_file, format)
    file_num = 1
    while True:
        batch = list(islice(records, records_per_file))
        if not batch:
            break
        output_file = f'{output_prefix}_{file_num}.{format}'
        SeqIO.write(batch, output_file, format)
        print(f'Wrote {len(batch)} records to {output_file}')
        file_num += 1

split_file('large.fasta', 'fasta', 1000, 'split')

On Python 3.12+, itertools.batched(records, records_per_file) yields the same fixed-size tuples without the manual while/islice loop.

Split by Sequence ID Prefix

Goal: Group sequences into separate files by a shared id prefix (sample or chromosome).

Approach: Route each record to a per-prefix open output handle while streaming, so no group is fully held in RAM.

Reference (BioPython 1.83+):

handles = {}
for record in SeqIO.parse('input.fasta', 'fasta'):
    prefix = record.id.split('_')[0]
    if prefix not in handles:
        handles[prefix] = open(f'{prefix}.fasta', 'w')
    SeqIO.write(record, handles[prefix], 'fasta')

for handle in handles.values():
    handle.close()

Batch Convert

for gb_file in Path('genbank/').glob('*.gb'):
    fasta_file = Path('fasta/') / gb_file.with_suffix('.fasta').name
    count = SeqIO.convert(str(gb_file), 'genbank', str(fasta_file), 'fasta')
    print(f'{gb_file.name} -> {fasta_file.name}: {count} records')

SeqIO.convert streams internally and never loads the whole file. GenBank-to-FASTA silently drops features, annotations, and qualifiers (FASTA stores only id, description, and sequence); see sequence-io/format-conversion before converting away annotated formats.

Parallel Processing

For CPU-bound per-file work, distribute whole files across processes. Each worker streams its own file, so peak memory is one file's records per process, not the whole directory.

from multiprocessing import Pool

def process_file(filepath):
    total = 0
    bp = 0
    for record in SeqIO.parse(filepath, 'fasta'):
        total += 1
        bp += len(record.seq)
    return {'file': filepath.name, 'count': total, 'total_bp': bp}

files = list(Path('data/').glob('*.fasta'))
with Pool(4) as pool:
    results = pool.map(process_file, files)

Use concurrent.futures.ThreadPoolExecutor instead for I/O-bound work (gzip decode, network filesystems); the GIL makes threads pointless for CPU-bound parsing.

Summary Statistics

Goal: Build a per-file CSV of counts and length stats for a directory.

Approach: Stream each file once, accumulating count, total, min, and max as integers rather than collecting a length list per file.

Reference (BioPython 1.83+):

import csv

summaries = []
for fasta_file in Path('data/').glob('*.fasta'):
    count = total = 0
    min_len = None
    max_len = 0
    for record in SeqIO.parse(fasta_file, 'fasta'):
        n = len(record.seq)
        count += 1
        total += n
        max_len = max(max_len, n)
        min_len = n if min_len is None else min(min_len, n)
    summaries.append({'file': fasta_file.name, 'sequences': count, 'total_bp': total,
                      'min_len': min_len or 0, 'max_len': max_len,
                      'avg_len': total / count if count else 0})

with open('summary.csv', 'w', newline='') as f:
    writer = csv.DictWriter(f, fieldnames=summaries[0].keys())
    writer.writeheader()
    writer.writerows(summaries)

Common Errors

Symptom	Cause	Fix
MemoryError / process killed on a directory	list(SeqIO.parse(...)) materializes every record at once	Stream the generator; iterate or sum(1 for _ in ...); never list() a large file
Counting/merge job runs for minutes on tens of millions of reads	SeqIO builds a full SeqRecord per read	Use pysam.FastxFile for linear iteration, or pyfastx for indexed access
ValueError: Duplicate key from index_db/to_dict	Same id appears in more than one merged file	Make ids unique (prefix by filename) or supply a key_function
Second loop over SeqIO.parse(...) yields nothing	The generator is one-pass and exhausts silently	Re-create the generator per pass, or use index_db for repeated access
index_db fails on a .gz file	Plain gzip is not seekable	Re-compress with bgzip; only BGZF is indexable (sequence-io/compressed-files)
Annotations missing after batch convert	GenBank-to-FASTA drops all features silently	Keep an annotated format, or extract needed qualifiers first

Related Skills

• read-sequences - parse, index, and index_db semantics for each file
• filter-sequences - apply per-record filters while streaming a batch
• sequence-statistics - N50 and length distributions across files
• format-conversion - batch format conversion and its data-loss traps
• compressed-files - BGZF vs plain gzip for indexable batch random access
• paired-end-fastq - keep R1/R2 synchronized when batch-filtering mates
• database-access/entrez-fetch - batch download sequences from NCBI