NGS Data Quality Control: FastQC & Trimmomatic Workflow Guide

1. Introduction

Next-Generation Sequencing (NGS) instruments generate billions of base pairs per run, but raw sequencing reads are inherently subject to technical artifacts: optical sensor noise, optical duplicates, phasing/pre-phasing decay, sequence-specific hexamer bias, and residual adapter read-through.

Quality control (QC) is the foundational gatekeeper of any computational genomics pipeline. Processing unfiltered reads directly through alignment or assembly leads to spurious mismatch inflation, false-positive variant calls, and skewed expression quantification.

This production guide provides an automated end-to-end paired-end bash pipeline combining FastQC, Trimmomatic, and MultiQC, accompanied by an in-depth diagnostic guide for interpreting FastQC metrics.

2. Production Bash Pipeline: FastQC -> Trimmomatic -> MultiQC

The following production bash script processes paired-end Illumina FASTQ files through pre-QC, adapter/quality trimming, post-QC, and unified MultiQC report generation.

#!/usr/bin/env bash
# ==============================================================================
# NGS Paired-End Quality Control Pipeline (FastQC -> Trimmomatic -> MultiQC)
# ==============================================================================
set -euo pipefail

# Configuration and Resource Allocation
THREADS=8
ADAPTER_FILE="/usr/share/trimmomatic/TruSeq3-PE.fa" # Update to local adapter path
RAW_DIR="data/raw_fastq"
QC_PRE_DIR="results/qc_pre"
TRIMMED_DIR="results/trimmed_fastq"
QC_POST_DIR="results/qc_post"
MULTIQC_DIR="results/multiqc_report"

# Create output directories
mkdir -p "${QC_PRE_DIR}" "${TRIMMED_DIR}" "${QC_POST_DIR}" "${MULTIQC_DIR}"

echo "=== [Step 1/4] Running Initial FastQC on Raw Data ==="
fastqc \
  --threads "${THREADS}" \
  --outdir "${QC_PRE_DIR}" \
  "${RAW_DIR}"/*_R1.fastq.gz "${RAW_DIR}"/*_R2.fastq.gz

echo "=== [Step 2/4] Running Trimmomatic on Paired-End Samples ==="
for R1 in "${RAW_DIR}"/*_R1.fastq.gz; do
  # Extract base sample identifier (e.g., sample01_R1.fastq.gz -> sample01)
  SAMPLE=$(basename "${R1}" _R1.fastq.gz)
  R2="${RAW_DIR}/${SAMPLE}_R2.fastq.gz"

  echo "Trimming Sample: ${SAMPLE}"

  trimmomatic PE \
    -threads "${THREADS}" \
    -phred33 \
    "${R1}" "${R2}" \
    "${TRIMMED_DIR}/${SAMPLE}_R1_paired.fastq.gz" \
    "${TRIMMED_DIR}/${SAMPLE}_R1_unpaired.fastq.gz" \
    "${TRIMMED_DIR}/${SAMPLE}_R2_paired.fastq.gz" \
    "${TRIMMED_DIR}/${SAMPLE}_R2_unpaired.fastq.gz" \
    ILLUMINACLIP:"${ADAPTER_FILE}":2:30:10:2:keepBothReads \
    LEADING:3 \
    TRAILING:3 \
    SLIDINGWINDOW:4:15 \
    MINLEN:36
done

echo "=== [Step 3/4] Running Post-Trimming FastQC ==="
fastqc \
  --threads "${THREADS}" \
  --outdir "${QC_POST_DIR}" \
  "${TRIMMED_DIR}"/*_paired.fastq.gz

echo "=== [Step 4/4] Aggregating Reports with MultiQC ==="
multiqc \
  --outdir "${MULTIQC_DIR}" \
  --title "NGS QC Aggregate Report" \
  --filename "ngs_qc_summary.html" \
  "${QC_PRE_DIR}" "${TRIMMED_DIR}" "${QC_POST_DIR}"

echo "=== QC Pipeline Completed Successfully! ==="
echo "Inspect report at: ${MULTIQC_DIR}/ngs_qc_summary.html"

3. FastQC Diagnostic Metrics Guide

FastQC runs a battery of modular quality assessments. Understanding which warnings are fatal and which are expected biological phenomena prevents over-trimming and data loss.

3.1 Per Base Sequence Quality & Phred Scores

Illumina quality scores ($Q$) are logarithmic representations of the probability of an erroneous base call ($P$):

$$Q = -10 \log_{10}(P)$$

Quality (Phred Q)
   40 +---------------------------------------------------------+
      |  [GREEN ZONE: High Quality > Q28]                       |
   30 |===================================\                     |
      |                                    \=========           |
   20 + - - - - - - - - - - - - - - - - - - \ - - - - - - - - - + [YELLOW ZONE: Q20 - Q28]
      |                                      \                  |
   10 + - - - - - - - - - - - - - - - - - - - \ - - - - - - - - + [RED ZONE: < Q20]
      |                                        \========        |
    0 +---------------------------------------------------------+
      1    10   20   30   40   50   60   70   80   90  100  150  (Read Cycle / Position)

Diagnostic Rule: A modest quality drop towards the 3’ end of the read is normal due to fluorophore consumption and laser degradation over dozens of sequencing cycles. When the median Phred drops below Q20, sliding window trimming is mandatory.

3.2 Per Base Sequence Content: The RNA-seq False Alarm

In an unbiased random genome library, the base proportions (A, T, G, C) should run completely parallel across all cycles. FastQC triggers an automatic Warning if difference $> 10%$ and a Failure if difference $> 20%$.

% Base Composition
  50% +---------------------------------------------------------+
      |  \  /                                                   |
  40% |   \/                                                    |
  30% |---/\----------------------------------------------------| A / T (~30%)
  20% |--/--\---------------------------------------------------| G / C (~20%)
  10% | /    \                                                  |
   0% +---------------------------------------------------------+
      1    5   10   15   20   25   30   35   40   45   50 (Cycles)
      [Hexamer Priming Bias]        [Stable Random Representation]

The RNA-Seq Phenomenon: In bulk RNA-seq, cDNA first-strand synthesis requires random hexamer priming. Hexamer binding is thermodynamically non-random, causing an inevitable A/T vs G/C spike across the first 10 to 12 bases.

• Action: If analyzing RNA-seq, DO NOT trim these bases unless specific adapter sequences are detected. Trimming the first 12 bases shifts alignment boundaries and discards valid biological reads without improving accuracy.

3.3 Sequence Duplication Levels: DNA-Seq vs. RNA-Seq

FastQC calculates the percentage of reads found multiple times in the library.

• In DNA-seq / WGS / ChIP-seq: High duplication indicates PCR over-amplification or low starting library complexity. Reads must be marked and removed using Picard MarkDuplicates or Sambamba.
• In RNA-seq: Highly expressed genes (e.g., GAPDH, ACTB, ribosomal RNAs) naturally generate millions of identical transcripts. High duplication in RNA-seq is expected biology—do not filter duplicate reads unless Unique Molecular Identifiers (UMIs) were incorporated during library prep.

4. Trimmomatic Parameter Ordering & Mechanics

The execution order of arguments in Trimmomatic is critical. Trimmomatic evaluates each read through the sequence of filters strictly from left to right:

[Raw Paired Reads]
       |
       v
1. ILLUMINACLIP  ---> Clips adapter read-through and palindromic fragments FIRST
       |
       v
2. LEADING        ---> Strips low quality (Q < 3) bases from 5' end
       |
       v
3. TRAILING       ---> Strips low quality (Q < 3) bases from 3' end
       |
       v
4. SLIDINGWINDOW  ---> Scans 4-base window, cuts when average Q < 15
       |
       v
5. MINLEN         ---> Discards any read trimmed shorter than 36 bp
       |
       v
[High-Quality Clean Reads: Paired & Unpaired Output]

Why ILLUMINACLIP Must Precede SLIDINGWINDOW

If you run SLIDINGWINDOW first, a low-quality region inside an adapter fragment will be cleaved, truncating the adapter sequence. The truncated adapter fragment will then fail to meet the minimum alignment seed score required by ILLUMINACLIP, leaving partial adapter sequences attached to the read that contaminate downstream alignments.

Parameter Breakdown

• ILLUMINACLIP:TruSeq3-PE.fa:2:30:10:2:keepBothReads:
  • 2: Seed mismatches allowed (max 2 mismatches).
  • 30: Palindrome clip threshold (score required to clip adapter in paired reads).
  • 10: Simple clip threshold (score required for single reads).
  • 2: Min adapter length to keep.
  • keepBothReads: Retains read information for downstream reverse alignment when adapter trimming palindromes.
• LEADING:3 and TRAILING:3: Removes extremely poor bases ($Q < 3$) from both termini.
• SLIDINGWINDOW:4:15: Scans in 4-base windows; truncates the read once average window quality falls below Phred 15 (~97% accuracy).
• MINLEN:36: Completely discards reads shorter than 36 bases, preventing ambiguous mapping of micro-fragments against the reference genome.

5. Experimental Use Cases & QC Profiles

Quality control thresholds vary significantly depending on the underlying genomic assay:

5.1 Whole Genome & Whole Exome Sequencing (WGS / WES)

• Primary Concerns: Low PCR duplication, uniform coverage depth across GC-extreme exons, and minimal adapter read-through.
• Action: Always mark and remove duplicate reads using Picard MarkDuplicates or Sambamba before calling germline or somatic single-nucleotide variants (SNVs).

5.2 RNA-Seq Transcriptomics

• Primary Concerns: Ribosomal RNA (rRNA) contamination, 3’ sequence degradation, and library strandedness.
• Action: Do not remove duplicate reads without UMIs. Expect and ignore FastQC failure in Per Base Sequence Content across cycles 1–12 resulting from random hexamer priming.

5.3 ChIP-Seq & Epigenomics

• Primary Concerns: Library complexity, cross-correlation metrics (NSC and RSC), and read enrichment at transcription factor binding sites.
• Action: High duplication severely distorts peak callers like MACS3. Filter duplicates and verify fragment length distribution before peak calling.

6. Frequently Asked Questions (FAQ)

Why does FastQC show a failure in Per Base Sequence Content for RNA-seq reads?

In RNA-seq, cDNA first-strand synthesis utilizes random hexamer priming. Hexamer binding is thermodynamically biased toward specific GC-rich motifs, producing a non-random nucleotide composition across the first 10–12 bases. This triggers an automated warning or failure in FastQC. This is a known technical feature of the assay and does not indicate adapter contamination or poor library quality.

Should duplicate reads always be removed in NGS quality control?

No. In DNA sequencing (WGS/WES) and ChIP-seq, duplicate reads usually stem from PCR amplification artifacts and should be marked and removed to avoid false-positive variant calls. In contrast, in RNA-seq, highly expressed genes naturally yield identical sequence fragments from genuine biological transcription. Removing duplicates in RNA-seq without UMIs severely distorts dynamic range and statistical power.

What is the optimal order for Trimmomatic parameters in paired-end mode?

Trimmomatic processes directives sequentially from left to right. You must always specify ILLUMINACLIP first, followed by LEADING, TRAILING, and SLIDINGWINDOW, concluding with MINLEN. Running quality-based trimming prior to adapter clipping truncates adapter fragments, preventing the seed algorithm from recognizing and stripping the full adapter sequence.