Module 5: Long-Read Quality Control: NanoPlot, Read Length Distributions & N50 Metrics

1. Interpreting Long-Read Quality Diagnostics

Unlike Illumina short reads where all reads have identical lengths ($150\text{ bp}$) and quality decays predictable across the cycle number, long-read sequencing exhibits highly skewed length distributions and independent single-molecule quality profiles.

Key metrics that define a successful long-read run include:

• Read Length N50: The sequence length such that $50\%$ of the total sequenced bases reside in reads equal to or longer than this value. For high-molecular-weight (HMW) DNA extractions, an $N50 \ge 25 - 40\text{ kb}$ is expected.
• Mean Quality Score ($\text{Q}$): The arithmetic mean of the Phred quality scores across all sequenced bases. A successful R10.4.2 run with SUP basecalling yields a mean $\text{Q} \ge 18 - 20$.
• Yield & Active Pore Count: The total gigabases ($\text{Gb}$) generated over time, showing flow cell pore decay kinetics and reagent degradation.

2. Generating Interactive Diagnostic Reports with NanoPlot

NanoPlot is the standard open-source quality control engine designed specifically for Oxford Nanopore and PacBio data:

# Run NanoPlot on raw basecalled FASTQ files
NanoPlot \
    --fastq basecalled_reads.fastq.gz \
    --threads 8 \
    --maxlength 100000 \
    --plots hex dot \
    --title "HG002 Long-Read Quality Audit" \
    --outdir ./qc_reports/

# Generate an instantaneous tabular summary with NanoStat
NanoStat --fastq basecalled_reads.fastq.gz > ./qc_reports/nanostat_summary.txt

3. Mathematical Derivation of Read Length N50

Given an array of read lengths $L = [l_1, l_2, \dots, l_n]$ sorted in descending order ($l_1 \ge l_2 \ge \dots \ge l_n$), let $G = \sum_{i=1}^n l_i$ represent the total sequenced throughput. The $N50$ is the length $l_k$ such that:

$$\sum_{i=1}^{k-1} l_i < 0.5 \times G \quad \text{and} \quad \sum_{i=1}^k l_i \ge 0.5 \times G$$

# Pure Python implementation of N50 calculation
def calculate_n50(lengths):
    lengths.sort(reverse=True)
    total_bases = sum(lengths)
    cumulative = 0
    for length in lengths:
        cumulative += length
        if cumulative >= total_bases * 0.5:
            return length
    return 0