Bio Applied Ngs Fundamentals

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

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 not header:
                break
            sequence = f.readline().strip()
            f.readline()  # '+' separator line
            quality = f.readline().strip()
            yield header[1:], sequence, quality  # strip leading '@'
            count += 1
            if max_reads and count >= max_reads:
                break

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):
    """Compute per-position quality percentiles (FastQC 'per base quality')."""
    pos_quals = defaultdict(list)
    for _, _, qual in reads:
        for pos, qchar in enumerate(qual):
            pos_quals[pos].append(ascii_to_phred(qchar))
    positions = sorted(pos_quals)
    return {
        'positions': positions,
        'median': [np.median(pos_quals[p]) for p in positions],
        'q25': [np.percentile(pos_quals[p], 25) for p in positions],
        'q75': [np.percentile(pos_quals[p], 75) for p in positions],
    }

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.