# Bio Applied Ngs Fundamentals

> 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.

- Skill: `pavel-kravchenko/bio-applied-ngs-fundamentals` (Agent Skill)
- Install (CLI): `npx skillmds@latest add pavel-kravchenko/bio-applied-ngs-fundamentals`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/bio-applied-ngs-fundamentals/raw
- Safety review: pending
- Works with: Claude Code, Claude.ai, OpenAI Codex
- Category: Coding & Dev Tools
- Author: pavel-kravchenko (https://skillmd.com/u/pavel-kravchenko)
- Updated: 2026-09-17
- Page: https://skillmd.com/skills/pavel-kravchenko/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

| 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.

```python
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
```

| 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.

```python
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 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 |

```bash
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`.

```python
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.