Bio Fastq Quality

Version Compatibility

Reference examples tested with: BioPython 1.83+

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.

FASTQ Quality Scores

"Filter my FASTQ reads by quality score" -> Access, analyze, and filter per-base quality scores, trim low-quality bases, and generate per-position quality profiles.

• Python: SeqIO.parse() with record.letter_annotations['phred_quality'] (BioPython)
• CLI alternative: pysam.FastxFile (.get_quality_array() returns offset-removed Phred ints, but ALWAYS subtracts 33 - it cannot read Phred+64/Solexa correctly, so use it only on confirmed Phred+33 data; for legacy encodings stay on SeqIO with the explicit variant string)

The Governing Principle: Never Guess the Offset

A FASTQ file does not record which quality encoding it uses. The same ASCII byte means different Phred scores under different encodings, and the offsets (33 vs 64) differ by exactly 31. Choosing the wrong format string has two failure modes:

• LOUD (safe): a quality character lies outside the chosen parser's legal range -> ValueError naming the wrong QualityIO parser. The run stops.
• SILENT (dangerous): every character lies in the ASCII overlap region legal for both encodings -> no error, and every score is off by exactly 31. Reading Phred+33 data as fastq-illumina makes all scores 31 too LOW; reading Phred+64 data as fastq-sanger makes them 31 too HIGH. QC, filtering, and trimming silently operate on garbage scores.

Auto-detection is provably ambiguous: ASCII >= 64 is legal in every variant, so a high-quality Sanger file (all Q >= 31) and a low-quality Illumina-1.3 file can be byte-identical in their quality lines. There is no reliable way to detect the encoding from content alone (Biopython docs: this "cannot be detected reliably automatically"). Determine the encoding from the sequencing instrument and run metadata, not by guessing. Scanning for the minimum ASCII byte can only RULE OUT encodings (see "Ruling Out Encodings" below), never confirm one.

The Four FASTQ Encodings

Almost all data produced since 2011 is Phred+33 ('fastq'). Phred+64 and Solexa appear only in legacy datasets, but the cost of misreading them is silent corruption, so the encoding must be confirmed before parsing legacy files.

'fastq' is an alias for 'fastq-sanger'; both are Phred+33. The wrong choice produces a LOUD ValueError only when an out-of-range character appears, and SILENT 31-shifted scores otherwise.

Phred vs Solexa: Two Different Score Definitions

The Solexa encoding is not just a different offset; it uses a different score formula, which is why it needs a separate parser.

• Phred: Q = -10 * log10(P_error). Always >= 0.
• Solexa: Q = -10 * log10(P/(1 - P)) - an ODDS score. It goes NEGATIVE when P > 0.5 (floor -5), which is why Solexa quality strings include ASCII 59-63.

The two scales are asymptotically equal at high quality (rounded scores above ~Q10-13 are interchangeable) but diverge for poor-quality bases. The round trip Phred -> Solexa -> Phred is LOSSY in that low-quality region: Cock et al. (2010) note that Solexa scores 9 and 10 both map to Phred 10. Do not convert legacy Solexa data to Phred and back if the low-Q values matter.

Accessing Quality Scores

Quality scores live in record.letter_annotations['phred_quality'] as a list of ints. The attribute is letter_annotations (NOT per_letter_annotations, which does not exist). Solexa data parsed with 'fastq-solexa' stores record.letter_annotations['solexa_quality'] instead, and those values can be negative.

from Bio import SeqIO

for record in SeqIO.parse('reads.fastq', 'fastq'):
    quals = record.letter_annotations['phred_quality']
    print(record.id, quals[:10])

letter_annotations is length-locked to len(record.seq): assigning a list of the wrong length raises. To edit sequence and quality together, slice the record (slicing keeps qualities in sync) or build a fresh record.

Code Patterns

Calculate Average Quality per Read

for record in SeqIO.parse('reads.fastq', 'fastq'):
    quals = record.letter_annotations['phred_quality']
    print(f'{record.id}: {sum(quals) / len(quals):.1f}')

Filter Reads by Mean Quality

def high_quality_reads(records, min_avg_qual=20):
    for record in records:
        quals = record.letter_annotations['phred_quality']
        if sum(quals) / len(quals) >= min_avg_qual:
            yield record

records = SeqIO.parse('reads.fastq', 'fastq')
SeqIO.write(high_quality_reads(records, 25), 'filtered.fastq', 'fastq')

Filter by Minimum Quality at Any Position

def all_bases_above(records, min_qual=20):
    for record in records:
        if min(record.letter_annotations['phred_quality']) >= min_qual:
            yield record

Trim Low-Quality 3' End

Goal: Drop trailing bases below a quality cutoff while keeping qualities aligned to the trimmed sequence.

Approach: Walk inward from the 3' end to the first base that meets the cutoff, then slice the record; slicing a SeqRecord trims letter_annotations in step with the sequence.

Reference (BioPython 1.83+):

def trim_low_quality(record, min_qual=20):
    quals = record.letter_annotations['phred_quality']
    trim_pos = len(quals)
    for i in range(len(quals) - 1, -1, -1):
        if quals[i] >= min_qual:
            trim_pos = i + 1
            break
    return record[:trim_pos]

records = SeqIO.parse('reads.fastq', 'fastq')
SeqIO.write((trim_low_quality(r) for r in records), 'trimmed.fastq', 'fastq')

Sliding Window Quality Trim

Goal: Truncate a read at the first position where average quality in a sliding window drops below a threshold (the Trimmomatic SLIDINGWINDOW model).

Approach: Slide a fixed-size window across the quality list; when the window mean falls below the cutoff, slice the record at that position.

Reference (BioPython 1.83+):

def sliding_window_trim(record, window_size=5, min_avg_qual=20):
    quals = record.letter_annotations['phred_quality']
    for i in range(len(quals) - window_size + 1):
        if sum(quals[i:i + window_size]) / window_size < min_avg_qual:
            return record[:i] if i > 0 else None
    return record

Per-Position Quality Profile

Goal: Compute mean quality at each read position to spot systematic drops (typically 3' degradation).

Approach: Accumulate scores by position across reads, then average each position. NovaSeq binning (see below) makes per-position values cluster at a few discrete levels - expected, not a defect.

Reference (BioPython 1.83+):

from collections import defaultdict

position_quals = defaultdict(list)
for record in SeqIO.parse('reads.fastq', 'fastq'):
    for i, q in enumerate(record.letter_annotations['phred_quality']):
        position_quals[i].append(q)

for pos in sorted(position_quals)[:20]:
    quals = position_quals[pos]
    print(f'Position {pos}: mean={sum(quals) / len(quals):.1f}')

Count Reads by Quality Threshold

thresholds = [20, 25, 30, 35]
counts = {t: 0 for t in thresholds}
for record in SeqIO.parse('reads.fastq', 'fastq'):
    avg = sum(record.letter_annotations['phred_quality']) / len(record.seq)
    for t in thresholds:
        if avg >= t:
            counts[t] += 1