2. Read Trimming and Filtering (fastp) - AreliDGS/EPP622-2026-SNPCalling GitHub Wiki

2. Read Trimming and Filtering (fastp)

This section documents the configuration, coding, and execution of the read trimming pipeline. Using fastp, this allows to remove sequencing adapters, filter out low-quality bases, and poly-G tail trimming to clean the simulated datasets before the genome alignment.


Workspace and coding

Follow these step-by-step terminal commands to establish the script execution infrastructure:

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

# 2. Create and open the trimming shell script file using the vi editor
vi scripts/02_read_trimming.sh

# 3. Grant execution permissions to the script to turn it into an executable utility
chmod +x scripts/02_read_trimming.sh

Execution script code

The code scripts/02_read_trimming.sh intent processes all my 4 samples using fastp, and immediately executes a post-trimming FastQC and MultiQC analysis loop on the resulting processed files.

#!/bin/bash

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

# Define paths
OUTPUT_TRIM="results/02_trimmed"
OUTPUT_QC="results/02_fastqc_trimmed"

# Create required subdirectories for clean data organization
mkdir -p "$OUTPUT_TRIM"
mkdir -p "$OUTPUT_QC"

# === Phase 1: Running fastp trimming across all samples ===

# Process each sample with fastp. We have already identified a Poly-G/Adapter contamination, use `--trim_poly_g` for detection. Add `--detect_adapter_for_pe`; adapter detection for paired-end data. 

for sample in A B C D; do
    echo "Processing Sample ${sample}..."
    fastp \
        -i data/raw/Sample_${sample}_R1.fastq.gz \
        -I data/raw/Sample_${sample}_R2.fastq.gz \
        -o "${OUTPUT_TRIM}/Sample_${sample}_R1.trimmed.fastq.gz" \
        -O "${OUTPUT_TRIM}/Sample_${sample}_R2.trimmed.fastq.gz" \
        --html "${OUTPUT_TRIM}/fastp_report_Sample_${sample}.html" \
        --json "${OUTPUT_TRIM}/fastp_report_Sample_${sample}.json" \
        --trim_poly_g \
        --detect_adapter_for_pe \
        --thread 4
done

# === Phase 2: Running Post-Trimming Quality Control ===

# Now we execute FastQC on the cleaned outputs (post-trimming)
fastqc -o "$OUTPUT_QC" -t 4 ${OUTPUT_TRIM}/*.trimmed.fastq.gz

# Compile our report, comparative post-trimming MultiQC evaluation report
multiqc "$OUTPUT_QC" -o results/ -n 02_multiqc_trimmed_report

echo "Trimming and post-QC processing completed successfully."

Usage: 02_read_trimming.sh

# 1. Don't forget to activate the Conda environment containing fastp, fastqc, and multiqc
conda activate snps

# 2. Run the automated trimming script from the root workspace
./scripts/02_read_trimming.sh
# OR: bash scripts/02_read_trimming.sh

# OPTIONAL run in the background using nohup
# nohup: Prevents the script from being killed if the terminal closes
# > 02_fastp.log: Save as readable text log file
# 2>&1: To save the outputs and error logs into the same file
# &: Push the task to background execution instantly
nohup ./scripts/02_read_trimming.sh > results/02_fastp.log 2>&1 &

# 3. To see the real-time progress
tail -f results/02_fastp.log


Downloading post-trimming reports

To compare the quality profile improvement, download the newly compiled interactive HTML report to your local desktop. Open the terminal on your local laptop and execute the following secure copy command:

scp -r [email protected]:/lustre/isaac24/proj/UTK0505/agonza84/EPP622-2026-SNPCalling/results/02_multiqc_trimmed_report.html ~/Desktop/

Results

Command for results extraction

MultiQC gives a rounded number for total reads and does not report total bases per sample. To find out those numbers use this:

# Execute sequence count and base extraction across trimmed datasets
for f in results/02_trimmed/*.trimmed.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/02_multiqc_trimmed_report.html

The post-trimming table below tracks metrics for each sample (R1 + R2 combined). Genome coverage shifts indicate the selective exclusion of low-quality bases and technical artifacts relative to the 12,157,105 bp reference genome.

Sample File Reads Before Reads After Length Before (bp) Length After (bp) Total Bases After
A Sample_A_R1Sample_A_R2 804,650804,650 804,650804,650 150150 150 150 120,685,966120,685,945
B Sample_B_R1Sample_B_R2 160,930160,930 160,930160,930 150 150 150 150 24,138,10124,138,092
C Sample_C_R1Sample_C_R2 804,650804,650 788,629788,629 150 150 125.9125.8 99,262,85899,240,549
D Sample_D_R1Sample_D_R2 804,650804,650 804,650804,650 150 150 150 150 120,667,096120,667,024

Main conclusions from trimming

  • Sample A: This sample was clean. It kept 100% of its reads and bases after running fastp, which proves it did not have any adapter contamination or low-quality issues that needed fixing.
  • Sample B: This sample did not lose a single read or base during filtering. This tells us that the quality of the data is good. The real issue is just a low volume of data, we can see a difference of the starting number of reads (~160k vs ~804k).
  • Sample C: This is the sample that had adapters and poly-G tails, that was successfully trimmed. This result on mean read length dropped from 150 bp to around 125 bp.
  • Sample D: This sample kept all of its reads, but its total bases dropped slightly. This means fastp selectively trimmed off individual low-quality bases and keeping the one of high quality.