Sequencing Reads and Quality Control

You already know how to read a FASTQ file and move around a shell, so this lesson is about the decisions that turn raw Illumina reads into clean input for downstream analysis. We will inspect paired-end output, quantify quality with FastQC and MultiQC, remove primers and adapters, screen out host DNA, and connect every trimming choice to what denoising expects.

Illumina Paired-End Output

Demultiplexing with bcl2fastq or BCL Convert splits raw base calls into per-sample FASTQ files, one for the forward read (R1) and one for the reverse read (R2), named like S01_S1_L001_R1_001.fastq.gz. Each read spans four lines: a header, the sequence, a plus separator, and one quality character per base. Quality is a Phred score Q = -10 log10(P) written as ASCII with a +33 offset, so Q30 means a 1-in-1000 error probability (99.9 percent accuracy) and Q20 means 1-in-100.

@M00123:45:000000000-A1B2C:1:1101:15589:1332 1:N:0:ACGTACGT
GTGCCAGCAGCCGCGGTAACGTAGGGTGCAAGCGTTAATCGGAATTACTGGGCGTAAAGC
+
CCCCCGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFFFFEEEDDCCBA@?<

FastQC and MultiQC

mkdir -p qc_raw
fastqc -t 4 -o qc_raw/ reads/*_R1_001.fastq.gz reads/*_R2_001.fastq.gz
multiqc qc_raw/ -n multiqc_raw -o qc_raw/

Started analysis of S01_R1_001.fastq.gz
Approx 25% complete for S01_R1_001.fastq.gz
Approx 100% complete for S01_R1_001.fastq.gz
Analysis complete for S01_R2_001.fastq.gz

[INFO   ]         multiqc : This is MultiQC v1.21
[INFO   ]          fastqc : Found 24 reports
[INFO   ]         multiqc : Report      : qc_raw/multiqc_raw.html

In the MultiQC report, the per-base sequence quality panel shows where mean Q drops along the read, and Illumina reverse reads (R2) almost always decay earlier and lower than R1. The adapter content plot flags read-through into Illumina adapters near the 3-prime end, while the sequence length distribution confirms whether reads are still uniform length after trimming. Read all three before choosing any cutoff, because the numbers you pick here propagate downstream.

Trimming and Host Removal

# Amplicon: strip 16S primers (515F / 806R) from R1 and R2
cutadapt -g GTGYCAGCMGCCGCGGTAA -G GGACTACNVGGGTWTCTAAT \
  --discard-untrimmed \
  -o trimmed/S01_R1.fastq.gz -p trimmed/S01_R2.fastq.gz \
  reads/S01_R1_001.fastq.gz reads/S01_R2_001.fastq.gz

# Quality- and length-trim (fastp; Trimmomatic is an equivalent choice)
fastp -i trimmed/S01_R1.fastq.gz -I trimmed/S01_R2.fastq.gz \
  -o clean/S01_R1.fastq.gz -O clean/S01_R2.fastq.gz \
  --cut_tail --cut_tail_mean_quality 20 --length_required 100 \
  --detect_adapter_for_pe

# Shotgun only: drop host reads by mapping to a reference (KneadData wraps this)
bowtie2 -x /db/GRCh38 -1 clean/S01_R1.fastq.gz -2 clean/S01_R2.fastq.gz \
  --very-sensitive -p 8 --un-conc-gz hostfree/S01_R%.fastq.gz -S /dev/null

Your QC read-outs set the truncation lengths for denoising: in DADA2 you pass truncLen (QIIME 2 uses --p-trunc-len-f and --p-trunc-len-r) to cut each read just before its quality curve collapses. Truncate too aggressively and paired reads no longer share the minimum 12 bp of overlap DADA2 needs to merge (leave a ~20 bp margin to be safe); truncate too little and error-rich tails inflate spurious ASVs. DADA2 also assumes primer-free reads of uniform length, which is exactly why cutadapt and fixed-length trimming must run first.