bio-applied-ngs-fundamentals

name	bio-applied-ngs-fundamentals
description	Decode Phred+33 FASTQ quality scores, compute FastQC-style per-position QC stats, and sliding-window trim reads in Python. Use when parsing FASTQ, decoding quality ASCII, or choosing Illumina/PacBio/Nanopore.
tool_type	python
primary_tool	FastQC

NGS Fundamentals

When to Use

• Choosing a sequencing platform (Illumina vs PacBio HiFi vs Oxford Nanopore) for a project.
• Parsing raw FASTQ files or decoding Phred+33 ASCII quality strings by hand.
• Reproducing or explaining FastQC modules (per-base quality, per-sequence quality, GC content, adapter content).
• Writing a custom sliding-window quality trimmer or adapter trimmer (Trimmomatic/fastp-style logic).
• Debugging "read count mismatch" or "quality score looks wrong" issues in an upstream pipeline.

Version Compatibility

FastQC ≥0.12, fastp ≥0.23, Trimmomatic ≥0.39, Python ≥3.10, NumPy ≥1.24. Phred+33 is universal on all instruments shipped since 2011 (Illumina CASAVA ≥1.8, all PacBio, all Nanopore).

Prerequisites

• pip install numpy
• Concepts: FASTA/FASTQ format, ASCII encoding, basic statistics (mean/percentile).
• Related skill: bio-read-qc-fastp-workflow for running trimming as a CLI step instead of by hand.

Platform Comparison

Feature	Illumina	PacBio HiFi	Oxford Nanopore
Read length	50–300 bp	10–25 kb	10 kb – 1 Mb+
Accuracy	~99.9%	~99.9%	~99% (R10.4.1 simplex)
Error type	Substitutions	Random (CCS averages indels)	Homopolymer indels (R9), substitutions (R10)
Throughput	Up to 6 Tb/run	~30 Gb/cell	50–200 Gb/cell
Best for	WGS, RNA-seq, ChIP-seq	De novo assembly, SVs	Structural variants, field diagnostics

Selection guide: High-coverage WGS/RNA-seq → Illumina. De novo assembly/complex SVs → PacBio HiFi. Rapid diagnostics/ultra-long reads → Nanopore. Best assemblies: hybrid Illumina + long-read.

Phred Quality Scores and FASTQ Parsing

Goal: convert between Phred quality scores, error probabilities, and the ASCII characters stored in FASTQ files, then parse a FASTQ file into (header, sequence, quality) records.

Approach: all modern platforms use Phred+33 encoding (ord(char) - 33). Each read is 4 lines: @header, sequence, +, quality string of identical length. Older Illumina CASAVA <1.8 used Phred+64 — check the FastQC encoding warning if scores look implausibly high or negative.

import math

def phred_to_error_prob(phred):
    """Convert a Phred quality score to a base-call error probability."""
    return 10 ** (-phred / 10)

def error_prob_to_phred(prob):
    """Convert an error probability back to a Phred quality score."""
    if prob <= 0:
        return 40  # cap at Q40
    return -10 * math.log10(prob)

def ascii_to_phred(char, offset=33):
    """Decode one FASTQ quality character to its Phred score (Phred+33 default)."""
    return ord(char) - offset

def phred_to_ascii(phred, offset=33):
    """Encode a Phred score back to its FASTQ ASCII quality character."""
    return chr(phred + offset)

def parse_fastq(filepath, max_reads=None):
    """
    Yield (header, sequence, quality_string) tuples from a FASTQ file.
    For gzipped input, pass a file handle opened with gzip.open(path, 'rt').
    """
    count = 0
    with open(filepath, 'r') as f:
        while True:
            header = f.readline().strip()
            if  header:
                
            sequence = f.readline().strip()
            f.readline()  
            quality = f.readline().strip()
             header[:], sequence, quality  
            count += 
             max_reads  count >= max_reads:
                

Phred	Error prob	Accuracy	ASCII char
10	10%	90%	+
20	1%	99%	5
30	0.1%	99.9%	?
40	0.01%	99.99%	I

QC Summary and Per-Position Quality (FastQC-style)

Goal: reproduce FastQC's headline numbers (read count, GC%, Q20/Q30 fraction) and the per-base quality profile, entirely in Python, for a stream of (header, seq, qual) reads.

Approach: accumulate per-position Phred scores in a dict[int, list[int]], then take the median and IQR at each position — this is exactly what the FastQC "per base sequence quality" plot shows.

from collections import Counter, defaultdict
import numpy as np

def compute_qc_summary(reads):
    """Generate a FastQC-like summary report from parsed FASTQ reads."""
    total_reads, total_bases, gc_count = 0, 0, 0
    lengths, mean_quals = [], []
    base_counts = Counter()

    for _, seq, qual in reads:
        total_reads += 1
        total_bases += len(seq)
        lengths.append(len(seq))
        gc_count += seq.count('G') + seq.count('C')
        base_counts.update(seq)
        mean_quals.append(np.mean([ascii_to_phred(c) for c in qual]))

    q20 = sum(1 for q in mean_quals if q >= 20)
    q30 = sum(1 for q in mean_quals if q >= 30)
    return {
        'total_reads': total_reads,
        'gc_content_pct': 100 * gc_count / total_bases,
        'mean_quality': np.mean(mean_quals),
        'pct_q20': 100 * q20 / total_reads,
        'pct_q30': 100 * q30 / total_reads,
        'base_counts': dict(base_counts),
    }

def per_position_quality(reads):
    
    pos_quals = defaultdict()
     _, _, qual  reads:
         pos, qchar  (qual):
            pos_quals[pos].append(ascii_to_phred(qchar))
    positions = (pos_quals)
     {
        : positions,
        : [np.median(pos_quals[p])  p  positions],
        : [np.percentile(pos_quals[p], )  p  positions],
        : [np.percentile(pos_quals[p], )  p  positions],
    }

FastQC module	Pass criteria	Common failure cause
Per-base quality	Q28+ across all positions	Quality drop at 3' end (normal)
Per-sequence quality	Peak at Q30+	Bimodal = subset of failed reads
GC content	Normal distribution	Shifted = contamination
Duplication level	<20%	High in targeted / PCR-heavy libs
Adapter content	<5% at ends	>10% → trim with Trimmomatic/fastp

fastqc sample_R1.fastq.gz sample_R2.fastq.gz -o qc_output/ -t 4

Sliding-Window Quality Trimming

Goal: trim low-quality 3' ends the way Trimmomatic's SLIDINGWINDOW does, without shelling out.

Approach: slide a fixed-size window from 5'→3'; the first window whose mean Phred drops below the threshold marks the trim point. Drop the read entirely if what remains is shorter than min_length.

def sliding_window_trim(sequence, quality_str, window_size=4, min_quality=20, min_length=36):
    """Trim a read's 3' end once a window's mean quality drops below min_quality."""
    quals = [ascii_to_phred(c) for c in quality_str]
    trim_pos = len(quals)
    for i in range(len(quals) - window_size + 1):
        if np.mean(quals[i:i + window_size]) < min_quality:
            trim_pos = i
            break
    trimmed_seq, trimmed_qual = sequence[:trim_pos], quality_str[:trim_pos]
    if len(trimmed_seq) < min_length:
        return None, None  # read too short after trimming — drop it
    return trimmed_seq, trimmed_qual

Pitfalls

• Coordinate systems: BED = 0-based half-open; VCF/GFF = 1-based inclusive — mixing causes off-by-one errors downstream.
• Phred+33 vs Phred+64: Old Illumina CASAVA <1.8 used +64. Check the FastQC encoding warning; ord('@') - 64 == 0 under +64 but ord('!') - 33 == 0 under +33.
• Quality drop at 3' end: Normal for sequencing-by-synthesis chemistry. Trim with fastp --cut_tail or Trimmomatic SLIDINGWINDOW:4:20.
• Paired-end read order: R1 and R2 must stay in sync — dropping a read from R1 requires dropping the same-index read from R2.
• error_prob_to_phred on prob=0: capped at Q40 above; don't feed it directly into math.log10 without the guard, or it raises ValueError.
• Multiple testing on QC metrics across many samples: apply FDR correction (Benjamini-Hochberg) before flagging outlier samples.

See Also

• bio-read-qc-fastp-workflow — running fastp as a production trimming/QC pipeline step.
• bio-read-qc-quality-reports — aggregating FastQC/MultiQC reports across samples.
• bio-sequence-io-fastq-quality — deeper FASTQ I/O and quality-filtering patterns.
• bio-read-alignment-bwa-alignment — next pipeline step after trimming/QC.