9.4 NGS Quality Metrics, Coverage Depth & Bioinformatics

9.4 NGS Quality Metrics, Coverage Depth & Bioinformatics

Quick Summary: Converting millions of raw optical or electronic sequencing signals into validated clinical diagnostic reports requires a structured bioinformatic pipeline and rigorous quality control checkpoints. The Phred quality score ($Q$-score) quantifies base-calling accuracy, with $Q30$ (99.9% accuracy, 1 in 1,000 error rate) serving as the clinical quality standard. Analytical sensitivity depends on coverage depth (fold coverage) and uniformity, ranging from $30\text{X}–100\text{X}$ for constitutional germline testing to $500\text{X}–2,000\text{X}+$ for somatic oncology panels detecting low-frequency variant allele fractions ($1%–5%$ VAF). The bioinformatic workflow progresses from primary base calling (FASTQ) through secondary alignment (SAM/BAM) and variant calling (VCF), concluding with tertiary clinical annotation against public disease databases (ClinVar, COSMIC, gnomAD) and AMP/ASCO/CAP or ACMG/AMP tier classifications.

---

1. Phred Quality Scores ($Q$-Score) & The FASTQ Architecture

The fundamental metric of base-calling reliability is the Phred Quality Score ($Q$), defined logarithmically as:

$Q = -10 \log_{10}(P_{\text{error}}) \iff P_{\text{error}} = 10^{-\frac{Q}{10}}$

where $P_{\text{error}}$ is the estimated probability that the assigned base call is incorrect.

+---------------------------------------------------------------------------------------------------------+
|                                   PHRED QUALITY SCORE CALIBRATION                                       |
+-------------+-----------------------------+-----------------------+-------------------------------------+
| Phred Score | Probability of Error (P)    | Base Call Accuracy    | Clinical Diagnostic Status          |
+-------------+-----------------------------+-----------------------+-------------------------------------+
| **Q10**     | $1\text{ in } 10$ ($0.10$)  | $90.0\%$              | Unacceptable; high error rate       |
| **Q20**     | $1\text{ in } 100$ ($0.01$) | $99.0\%$              | Minimum acceptable research cutoff  |
| **Q30**     | $1\text{ in } 1,000$ ($0.001$)| $99.9\%$            | **Clinical diagnostic standard**    |
| **Q40**     | $1\text{ in } 10,000$ ($0.0001$)| $99.99\%$          | High-accuracy benchmark             |
+-------------+-----------------------------+-----------------------+-------------------------------------+

Clinical Benchmark: Accredited clinical diagnostic NGS assays generally require that $\ge 80%$ to $85%$ of all sequenced bases in a run achieve $\ge Q30$.

The FASTQ File Format Architecture

Primary analysis software converts raw binary base calls (BCL files on Illumina) into FASTQ text files. Every single sequenced read is represented by a standardized 4-line entry:

@M01234:105:000000000-A8B9C:1:1101:14250:1820 1:N:0:TCGAAGTC
NATTCGATCGATCGATCGATCGATCGATCGATCGATCGATCGATCGA
+
#AAFFFFJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJF

1. Line 1 (Header @): Unique read identifier containing instrument ID (M01234), run number (105), flow cell ID, lane (1), tile (1101), X/Y cluster coordinates (14250:1820), read pair member (1 or 2), filter pass status (N), and the sample index barcode (TCGAAGTC).
2. Line 2 (Sequence): The raw nucleotide base calls (A, C, G, T, or N for ambiguous calls).
3. Line 3 (Separator +): Delimiter line, optionally repeating the header identifier.
4. Line 4 (Quality String): ASCII-encoded Phred quality scores corresponding 1-to-1 with the bases in Line 2. Under the standard Phred+33 offset (Sanger/Illumina 1.8+): $\text{ASCII Character Code} = Q + 33$ (e.g., $Q30 + 33 = 63$, corresponding to the ASCII character ?; $Q40 + 33 = 73$, corresponding to I).

---

2. Coverage Metrics: Depth, Breadth & Uniformity

Coverage metrics establish the statistical power and analytical limit of detection of an NGS assay.

                                    COVERAGE METRICS DEFINED
                                    
    Target Region:     ================================================ (Gene Exon: 500 bp)
    Read 1:            ------------>
    Read 2:               ------------>
    Read 3:                  ------------>
    Read 4:                        ------------>
    Read 5:                               ------------>
    Read 6:                                      ------------>
    
    - Read Depth (Fold Coverage): Number of overlapping reads covering a single base locus (e.g., 500X)
    - Breadth of Coverage: Percentage of the 500 bp target region covered at or above a threshold (e.g., 100% >= 100X)
    - Uniformity: Evenness of read distribution across all exons (minimal high peaks or dropout valleys)

Mathematical Calculation of Theoretical Fold Coverage ($X$)

$\text{Depth } (X) = \frac{N \times L}{G}$

where $N$ is the total number of reads generated, $L$ is the average read length in base pairs, and $G$ is the total target genome or panel size in base pairs.

Germline vs. Somatic Clinical Depth Standards

• Germline Testing (Hereditary Cancer, Mendelian Genetics):
  • Target Depth: $30\text{X}–100\text{X}$.
  • Rationale: Constitutional variants are present at expected frequencies of either $50%$ (heterozygous) or $100%$ (homozygous). At $50\text{X}$ depth, a heterozygous allele is supported by $\sim 25$ independent reads, providing overwhelming statistical confidence for variant calling.
• Somatic Oncology Testing (Solid Tumors & Hematologic Malignancies):
  • Target Depth: $500\text{X}–2,000\text{X}+$.
  • Rationale: Clinical tumor specimens are heterogeneous mixtures containing non-neoplastic stromal cells, infiltrating lymphocytes, tumor aneuploidy, and subclonal driver mutations. Somatic panels must reliably detect mutations at $1%–5%$ Variant Allele Frequency (VAF). At $1,000\text{X}$ depth, a $2%$ VAF subclonal mutation is supported by 20 independent reads, meeting statistical thresholds for confident clinical reporting.
• Liquid Biopsy / Cell-Free Tumor DNA (ctDNA) & Minimal Residual Disease (MRD):
  • Target Depth: $5,000\text{X}–20,000\text{X}+$ combined with Unique Molecular Identifiers (UMIs) to detect mutant fragments at $0.01%–0.1%$ VAF.

---

3. The Secondary Bioinformatics Pipeline: FASTQ to BAM

Secondary analysis transforms raw, unaligned reads into aligned, coordinate-sorted, deduplicated alignments.

                     THE SECONDARY BIOINFORMATICS PIPELINE WORKFLOW
                     
   +-------------------------------------------------------------------------+
   | 1. Quality Control & Trimming: FastQC -> Adapter Trimming (Trimmomatic) |
   +------------------------------------+------------------------------------+
                                        |
                                        v
   +-------------------------------------------------------------------------+
   | 2. Read Alignment: BWA-MEM maps reads to Reference Genome (hg19 / hg38) |
   |    Generates Sequence Alignment/Map (SAM) -> Binary Alignment Map (BAM) |
   +------------------------------------+------------------------------------+
                                        |
                                        v
   +-------------------------------------------------------------------------+
   | 3. Post-Alignment Processing: Coordinate Sorting & Indexing (BAI)        |
   |    PCR Duplicate Marking/Removal (Picard) + Base Recalibration (BQSR)   |
   +------------------------------------+------------------------------------+
                                        |
                                        v
   +-------------------------------------------------------------------------+
   | 4. Variant Calling: GATK HaplotypeCaller (Germline) / Mutect2 (Somatic) |
   |    Generates Variant Call Format (VCF) File                             |
   +-------------------------------------------------------------------------+

Key Pipeline Stages and Algorithms

1. Read Alignment (Mapping):
   • Reads are mapped to the human reference genome—GRCh37 (hg19) or GRCh38 (hg38).
   • BWA-MEM (Burrows-Wheeler Aligner): The clinical standard for short-read DNA alignment, utilizing the Burrows-Wheeler transform and Ferragina-Manzini (FM) index to rapidly align millions of reads with high mapping quality ($MQ$).
   • Splice-Aware Aligners (STAR, HISAT2): Used for RNA-sequencing to allow reads to span large non-coding introns across exon-exon junctions.
2. File Formats (SAM, BAM, BAI):
   • SAM (Sequence Alignment/Map): Human-readable, tab-delimited text file containing alignment coordinates, mapping quality scores, and CIGAR strings (e.g., 100M = 100 bp match; 50M2I48M = 2 bp insertion).
   • BAM (Binary Alignment Map): Highly compressed, BGZF binary representation of SAM. Coordinates must be sorted numerically by chromosome and position.
   • BAI (BAM Index): Accompanying index file enabling visualization software like the Integrated Genomics Viewer (IGV) to instantly retrieve and render reads at any genomic locus.
3. PCR Duplicate Removal (Deduplication):
   • Multiple reads sharing identical 5' mapping coordinates arise from in vitro PCR amplification of the same original DNA fragment during library prep.
   • Algorithms (Picard MarkDuplicates, Sambamba) flag duplicate reads, retaining only a single read with the highest quality score for variant calling. This prevents PCR amplification artifacts from artificially inflating coverage depth or skewing variant allele frequencies.

---

4. Variant Call Format (VCF) Architecture & Format Fields

Variant calling algorithms (e.g., GATK HaplotypeCaller for germline, Mutect2 / Strelka2 for somatic) analyze aligned BAM files and output candidate variants in a standardized VCFv4.2 file.

                                   ANATOMY OF A CLINICAL VCF ENTRY
                                   
  #CHROM  POS       ID          REF  ALT  QUAL  FILTER  INFO              FORMAT        PATIENT_TUMOR
  chr7    140453136 rs121913342 T    A    1500  PASS    DP=1000;AF=0.25   GT:AD:DP:GQ   0/1:750,250:1000:99

Dissecting the VCF Fields

• #CHROM & POS: Chromosome and 1-based coordinate (e.g., chr7:140453136, the locus of the BRAF V600E mutation).
• REF & ALT: Reference nucleotide (T) and alternate mutant nucleotide (A).
• FILTER: PASS indicates the variant met all clinical quality thresholds; failures list specific flags (e.g., LowQual, StrandBias, LowDepth).
• FORMAT and Sample Genotype Sub-Fields:
  • GT (Genotype): 0/0 = Homozygous Reference; 0/1 = Heterozygous Alternate; 1/1 = Homozygous Alternate.
  • AD (Allele Depth): Raw count of reads supporting each allele in the order Ref_Count,Alt_Count (e.g., 750,250).
  • DP (Depth): Total filtered read depth at the locus (1000).
  • GQ (Genotype Quality): Phred-scaled confidence in the genotype call (99 = maximum confidence).

Calculating Variant Allele Frequency (VAF)

Variant Allele Frequency represents the proportion of sequencing reads supporting the mutant allele:

$\text{VAF} = \frac{\text{Alt Read Depth}}{\text{Ref Read Depth} + \text{Alt Read Depth}} = \frac{\text{AD}[1]}{\text{AD}[0] + \text{AD}[1]} = \frac{250}{750 + 250} = \mathbf{0.25 \quad (25.0\%)}$

---

5. Tertiary Analysis, Annotation & Clinical Variant Interpretation

Tertiary analysis interprets called variants within clinical and pathological contexts.

+---------------------------------------------------------------------------------------------------------+
|                               CLINICAL VARIANT CLASSIFICATION FRAMEWORKS                                |
+---------------------+-----------------------------------+-----------------------------------------------+
| Framework           | Tier / Classification Levels      | Diagnostic & Therapeutic Actionability        |
+---------------------+-----------------------------------+-----------------------------------------------+
| **AMP / ASCO / CAP**| **Tier I:** Strong Clinical Sig.  | Biomarkers predicting response to FDA-approved|
| **Somatic Oncology**| (Therapeutic, Diagnostic, Prog.)  | therapies in standard-of-care guidelines.     |
| **System (2017)**   +-----------------------------------+-----------------------------------------------+
|                     | **Tier II:** Potential Clin. Sig. | Investigational therapies in clinical trials; |
|                     |                                   | off-label therapeutic indications.            |
|                     +-----------------------------------+-----------------------------------------------+
|                     | **Tier III:** Unknown Sig. (VUS)  | Variants of Uncertain Significance.           |
|                     +-----------------------------------+-----------------------------------------------+
|                     | **Tier IV:** Benign / Likely Ben. | Common population polymorphisms (gnomAD >1%). |
+---------------------+-----------------------------------+-----------------------------------------------+
| **ACMG / AMP**      | **Class 5:** Pathogenic           | >99% certainty of causing disease.            |
| **Germline System** +-----------------------------------+-----------------------------------------------+
| **(2015)**          | **Class 4:** Likely Pathogenic    | >90% certainty of causing disease.            |
|                     +-----------------------------------+-----------------------------------------------+
|                     | **Class 3:** VUS                  | Uncertain Significance (insufficient evidence)|
|                     +-----------------------------------+-----------------------------------------------+
|                     | **Class 2:** Likely Benign        | >90% probability of being benign.             |
|                     +-----------------------------------+-----------------------------------------------+
|                     | **Class 1:** Benign               | Documented non-pathogenic polymorphism.       |
+---------------------+-----------------------------------+-----------------------------------------------+

Essential Annotation Databases

• Population Polymorphism Databases: gnomAD (Genome Aggregation Database) — variants with population allele frequencies $>1%$ in gnomAD are generally classified as benign polymorphisms rather than causative germline disease alleles.
• Somatic Cancer Knowledgebases: COSMIC (Catalogue of Somatic Mutations in Cancer), OncoKB, CIViC.
• Germline Mutation Archives: ClinVar, HGMD (Human Gene Mutation Database).

---

6. Clinical NGS Quality Assurance & Bioinformatic Failure Modes

1. Sample Cross-Contamination: Evaluated using bioinformatic tools such as VerifyBamID, which detect excess heterozygous base calls at invariant homozygous single nucleotide loci. Sample contamination exceeding $1%–2%$ requires specimen repeat.
2. GC Coverage Bias: Regions with extreme GC content ($>70%$ GC, such as the CEBPA or TERT promoter; or $<25%$ GC) amplify poorly during PCR, resulting in severe coverage valleys or complete allele dropouts.
3. Pseudogene Interference: Highly homologous pseudogenes (e.g., PMS2 vs. PMS2CL, SMN1 vs. SMN2, CYP2D6 vs. CYP2D7) share $>98%$ sequence identity with functional genes. Short $150\text{ bp}$ NGS reads frequently misalign to the pseudogene, creating false-positive mutation calls or masking true pathogenic variants. Resolution requires long-read sequencing (PacBio HiFi), Multiplex Ligation-dependent Probe Amplification (MLPA), or specialized pseudogene-masking bioinformatic pipelines.
4. FFPE Cytosine Deamination Artifacts: Formalin fixation induces hydrolytic deamination of unmethylated cytosines into uracils ($C \rightarrow U$). During library PCR amplification, DNA polymerase reads uracil as thymine, creating artificial $C>T$ / $G>A$ transition artifacts at low variant allele frequencies ($1%–5%$). Treating FFPE DNA with Uracil-DNA Glycosylase (UDG) prior to library construction excises uracil bases, eliminating false-positive deamination calls.