1. Initial Quality Control (FastQC MultiQC) - AreliDGS/EPP622-2026-SNPCalling GitHub Wiki

1. Initial Quality Control (FastQC + MultiQC)

This section documents the process of preparing, coding, and executing the baseline quality control assessment across all raw simulated Illumina paired-end datasets (Samples A, B, C, and D) using FastQC and MultiQC.

Workspace and coding

To start our analysis, we will set up a workspace scripts/ directory to store our pipeline code. Instead of repeatedly typing long command sequences into the terminal, we compile all instructions into a single file, making our pipeline faster and more reliable. Follow these steps:

# 1. Make sure you are in the project root directory
cd /lustre/isaac24/proj/UTK0505/agonza84/EPP622-2026-SNPCalling

# 2. Create the directory `scripts/`
mkdir -p scripts

# 3. Create and open the shell script file using the vi editor
# Inside vi, press 'i' to enter Insert mode, paste the script code, 
# then press 'ESC' and type ':wq' to save and exit.
vi scripts/01_fastqc_multiqc.sh

# 4. Give execute permissions to the script
chmod +x scripts/01_fastqc_multiqc.sh

Execution script code

#!/bin/bash

# Exit immediately if any command, undefined variable, or pipe fails
set -euo pipefail

# ========= Step 1 =========
# Create the subdirectory for QC output files
mkdir -p results/01_fastqc_raw

# ========= Step 2 =========
# Execute FastQC across all raw datasets, by calling `fastqc`
# -o: Specifies the directory for output files
# -t 4: Allocates 4 parallel processing threads to speed up runtime
# data/raw/*.fastq.gz: Indicates the location and format of the input files
fastqc -o results/01_fastqc_raw -t 4 data/raw/*.fastq.gz

# ========= Step 3 =========
# Call `multiqc` to compile all diagnostic summaries into a single report
# results/01_fastqc_raw/: Targets the folder containing the FastQC logs
# -o results/: Sets the output target directory 
# -n 01_multiqc_raw_report: Defines the name of the generated HTML report
multiqc results/01_fastqc_raw/ -o results/ -n 01_multiqc_raw_report

Usage: 01_fastqc_multiqc.sh

# 1. Activate the Conda environment containing the bioinformatic software tools
conda activate snps

# 2. Run the script
./scripts/01_fastqc_multiqc.sh
# OR: bash scripts/01_fastqc_multiqc.sh

Results

  • MultiQC report name: 01_multiqc_raw_report.html

The interactive MultiQC summary is generated as an HTML file (results/01_multiqc_raw_report.html). It must be transferred from the remote cluster to your local computer for visualization in a web browser.

Open a new terminal session on your local laptop and execute the following command:

# Copy the entire results folder from the remote server to your local desktop
scp -r [email protected]:/lustre/isaac24/proj/UTK0505/agonza84/EPP622-2026-SNPCalling/results/ ~/Desktop/

Command for results extraction

MultiQC gives a rounded number for total reads and does not report total bases per sample. To extract exact numbers from compressed FASTQ structures, the following bash loop was executed relative to the raw data subdirectory:

# Execute sequence count and base extraction across raw datasets
for f in data/raw/*.fastq.gz; do
  echo -ne "$f\t"
  zcat "$f" | awk 'NR%4==2{r++; b+=length($0)} END{printf "%d\t%d\t%.1f\n", r, b, b/r}'
done | column -t
  • MultiQC report path and name: results/multiqc_raw_report.html

The statistics below summarize the metrics for each sample, merging R1 + R2 allocations. Genome coverage was calculated based on the Saccharomyces cerevisiae reference genome size of 12,157,105 bp.

Sample Total Reads Total Bases Duplication (%) GC content (%) Quality Issues?
A 1,609,300 241,395,000 1.8-1.9 38 No problems
B 321,860 48,279,000 0.7 38 Low coverage
C 1,609,300 241,395,000 3.7-3.8 42 Poly-G & Adapter contamination
D 1,609,300 241,395,000 1.1-1.3 39 Low quality (Low Phred scores)

Visualizations

The plot below illustrates the distribution of average Phred quality scores per read across all raw datasets, revealing discrepancies, such as:

  • Baseline quality: Samples A, B, and C exhibit ideal, normal distribution curves peaking sharply at a Phred score of Q38. This signifies a low sequencing error probability.
  • Low-coverage diagnostic: Sample B is compressed on the Y-axis. This represents a low coverage defect, related to a library with less raw molecular data despite maintaining high sequence accuracy.
  • Low-quality: Sample D shows a prominent downward shift toward lower average Phred scores. This pattern confirms a low-quality library.
  • Adapter modification signal: Sample C's curve remains within high-quality boundaries but shows a distinct widening at the peak baseline. This is correlated with high duplicates (3.7-3.8%) and synthetic adapter/poly-G expansions.

Figure 1: Per sequence quality scores distribution

Moreover the plot below illustrates the cumulative percentage of sequencing reads contaminated with artificial adapter sequences as a function of the nucleotide position along the 150 bp read length. Some points to highligth are:

  • Flawless: The curves for samples A, B, and D remain completely flat across the entire 150 bp sequencing read. This baseline indicates optimal insert fragment lengths during library preparation.
  • Adapter Contamination: Conversely, Sample C displays an anomaly. The cumulative adapter contamination baseline takes off vertically around position 53 bp, peaking at 35% of sequences.

Figure 2: Adapter content cumulative distribution