Applied - FASTA and FASTQ Workflows

Core Principle: High-throughput sequence processing pipelines stream FASTA/FASTQ records, execute base-quality filtering, and calculate read composition metrics efficiently without loading entire datasets 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

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

Exercise

Biopython is not installed in this sandbox, so parse the FASTQ yourself: four lines per record, the header, the sequence, a + line and the quality string in Phred+33. Implementing the trimming loop by hand is the point -- it is what trim_quality in the lesson does under the wrapper.

Write process(text), returning a list of output lines. The first input line holds four whitespace-separated values: the minimum mean quality, the window size, the minimum retained length, and the adapter sequence, which may be the literal none. Read the records that follow four lines at a time. A record whose sequence and quality strings differ in length is malformed: count it and skip it. Otherwise trim the adapter first -- search the sequence for the adapter's first eight characters, or the whole adapter when it is shorter, and cut from the first occurrence -- then trim quality from the right, shrinking the read while at least window bases remain and the mean quality of the last window bases is below the threshold. Discard a read whose surviving length is below the minimum.

Print reads=<n> malformed=<n> adapter_trimmed=<n> quality_trimmed=<n> kept=<n>, counting a read as adapter trimmed when the adapter was found and as quality trimmed when the quality loop removed at least one base. Then print mean_len=<v>, the mean surviving length over the kept reads rounded to two decimals. Then one line per kept read in input order, holding the read name -- the first whitespace-separated token of the header with the leading @ removed -- the surviving length, and the mean quality of the surviving bases rounded to two decimals, tab separated. If no read is kept, print no reads kept in place of the mean length line and print no read lines.