Applied - FASTA and FASTQ Workflows

In one line: streaming, filtering, trimming and summarising sequence files without loading them into memory.

Uses: SeqIO · SeqRecord · Seq · Bio.SeqUtils · DataFrame

The streaming principle

A FASTQ file can be 50 GB. SeqIO.parse is a lazy iterator and SeqIO.write accepts a generator, so a filter is constant-memory regardless of file size:

import gzip
from Bio import SeqIO

def quality_filter(records, min_mean_q=25, min_len=50, max_n=2):
    for rec in records:
        q = rec.letter_annotations["phred_quality"]
        if len(rec) < min_len:
            continue
        if sum(q) / len(q) < min_mean_q:
            continue
        if str(rec.seq).upper().count("N") > max_n:
            continue
        yield rec

with gzip.open("in.fastq.gz", "rt") as fin, gzip.open("out.fastq.gz", "wt") as fout:
    n = SeqIO.write(quality_filter(SeqIO.parse(fin, "fastq")), fout, "fastq")
print(f"kept {n:,} reads")

Generator in, generator out. "rt"/"wt" for text mode — binary mode gives you a confusing parse error.

Trimming

def trim_quality(rec, min_q=20, window=4):
    """Trailing-window trim, like Trimmomatic's SLIDINGWINDOW."""
    q = rec.letter_annotations["phred_quality"]
    end = len(q)
    while end >= window and sum(q[end-window:end]) / window < min_q:
        end -= 1
    return rec[:end]

def trim_adapter(rec, adapter="AGATCGGAAGAG", min_overlap=8):
    s = str(rec.seq)
    i = s.find(adapter[:min_overlap])
    return rec[:i] if i != -1 else rec

Slicing a SeqRecord slices the quality string too — that automatic co-slicing is the whole reason to work with records rather than parallel lists. See SeqRecord.

For production read trimming, use fastp, cutadapt or Trimmomatic. They are faster, handle paired-end synchronisation, and are better tested. Write your own only when you need something unusual.

Paired-end reads

def paired(r1_path, r2_path):
    with gzip.open(r1_path, "rt") as f1, gzip.open(r2_path, "rt") as f2:
        for a, b in zip(SeqIO.parse(f1, "fastq"), SeqIO.parse(f2, "fastq")):
            assert a.id.split()[0] == b.id.split()[0], f"desync: {a.id} vs {b.id}"
            yield a, b

def filter_pairs(r1, r2, out1, out2, **kw):
    kept = 0
    with gzip.open(out1, "wt") as o1, gzip.open(out2, "wt") as o2:
        for a, b in paired(r1, r2):
            if passes(a, **kw) and passes(b, **kw):     # BOTH must pass
                SeqIO.write(a, o1, "fastq"); SeqIO.write(b, o2, "fastq")
                kept += 1
    return kept

The assertion is not optional. Filtering R1 and R2 independently desynchronises them, and every downstream aligner assumes the nth record of R1 pairs with the nth of R2. A desynced pair of files aligns "successfully" and produces complete nonsense.

Random access

idx = SeqIO.index("genome.fasta", "fasta")
chr17 = idx["chr17"]
region = chr17.seq[7_668_401:7_687_550]      # 0-based! see [[Coordinate Systems]]
idx.close()

idx = SeqIO.index_db("genome.idx", ["chr1.fa", "chr2.fa"], "fasta")   # persistent

SeqIO.index stores file offsets only — a 3 GB genome costs a few MB of RAM. Contrast SeqIO.to_dict, which loads everything.

For heavy reference access, pysam.FastaFile with a samtools .fai index is faster and is what the rest of the ecosystem uses.

Summarising a FASTA

import pandas as pd
from Bio.SeqUtils import gc_fraction

def summarise(path):
    return pd.DataFrame([
        {"id": r.id,
         "length": len(r),
         "gc": gc_fraction(r.seq),
         "n_frac": str(r.seq).upper().count("N") / len(r),
         "description": r.description}
        for r in SeqIO.parse(path, "fasta")
    ])

df = summarise("contigs.fasta")

def assembly_stats(lengths):
    L = sorted(lengths, reverse=True)
    total = sum(L)
    cum = 0
    for i, x in enumerate(L, 1):
        cum += x
        if cum >= total / 2:
            return {"n_contigs": len(L), "total_bp": total, "n50": x, "l50": i,
                    "longest": L[0], "mean": total / len(L)}

print(assembly_stats(df["length"]))

N50 is the length such that contigs of that length or longer contain half the assembly. L50 is how many contigs that takes. Reporting N50 without total assembly size is meaningless — you can raise N50 by discarding short contigs.

The contamination check

import seaborn as sns
sns.set_theme(style="ticks")

fig, ax = plt.subplots(figsize=(6, 4))
sns.scatterplot(data=df, x="length", y="gc", size="n_frac",
                alpha=0.5, ax=ax, sizes=(5, 100))
ax.set_xscale("log")
ax.set(xlabel="Contig length (bp)", ylabel="GC fraction")

A length-vs-GC scatter of an assembly is the standard first contamination screen. A distinct cloud at a different GC fraction is usually a different organism — bacterial contamination in a eukaryotic assembly shows up immediately.

Common transformations

# format conversion
SeqIO.convert("in.fastq", "fastq", "out.fasta", "fasta")

# subset by an ID list
wanted = set(pd.read_csv("ids.txt", header=None)[0])
SeqIO.write((r for r in SeqIO.parse("all.fasta", "fasta") if r.id in wanted),
            "subset.fasta", "fasta")

# rename, dropping the description so headers stay clean
def rename(records, mapping):
    for r in records:
        r.id = mapping.get(r.id, r.id)
        r.description = ""            # ← or the old id appears twice
        yield r

# split into chunks
def chunks(records, size=1000):
    batch = []
    for r in records:
        batch.append(r)
        if len(batch) == size:
            yield batch; batch = []
    if batch:
        yield batch

for i, batch in enumerate(chunks(SeqIO.parse("all.fasta", "fasta"))):
    SeqIO.write(batch, f"chunk_{i:04d}.fasta", "fasta")

# deduplicate by sequence
def dedup(records):
    seen = set()
    for r in records:
        h = hash(str(r.seq))
        if h not in seen:
            seen.add(h); yield r

# six-frame translation
def six_frames(rec, min_aa=50):
    for strand, nuc in [(+1, rec.seq), (-1, rec.seq.reverse_complement())]:
        for frame in range(3):
            trimmed = nuc[frame: len(nuc) - (len(nuc) - frame) % 3]
            for i, prot in enumerate(trimmed.translate().split("*")):
                if len(prot) >= min_aa:
                    yield SeqRecord(prot, id=f"{rec.id}_s{strand}_f{frame}_{i}",
                                    description="")

Per-position quality, for a QC plot

import numpy as np

def quality_matrix(path, n_reads=100_000, max_len=300):
    q = np.full((n_reads, max_len), np.nan, dtype=np.float32)
    with gzip.open(path, "rt") as fh:
        for i, rec in enumerate(SeqIO.parse(fh, "fastq")):
            if i >= n_reads:
                break
            scores = rec.letter_annotations["phred_quality"][:max_len]
            q[i, :len(scores)] = scores
    return q

q = quality_matrix("reads.fastq.gz")
mean_q = np.nanmean(q, axis=0)          # per POSITION → axis=0

axis=0 because reads are rows and positions are columns, and you want one value per position. Getting this backwards gives you per-read quality, which is a different (also useful) plot. See Axes and Reductions and Applied - Quality Control Plots.

Performance

Biopython's parsers are correct and general, not fastest. For tens of millions of reads:

Tool	Advantage
pyfastx	much faster, built-in random access
pysam.FastxFile	fast, already a dependency in most pipelines
seqkit	CLI; usually the right answer for pure filter/subset/stats
fastp	CLI; QC + trimming + reporting in one pass

Use Biopython when you need its object model or format breadth. For "filter this 40 GB FASTQ", seqkit in a subprocess call will beat anything you write.

Common mistakes

• list(SeqIO.parse(...)) on a large file.
• Iterating a parse result twice.
• gzip.open in binary mode.
• Desynchronising paired reads.
• Duplicated FASTA headers from an uncleared description.
• SeqIO.to_dict on a genome.
• Wrong FASTQ quality offset.
• Ignoring soft-masked lowercase.
• N50 without assembly size.
• Not closing a SeqIO.index.

See also

SeqIO · SeqRecord · Applied - Quality Control Plots · Applied - Sequence Composition Analysis · Applied - Genomic File Formats · Coordinate Systems