Step 2 Read Trimming and Filtering (fastp) - alyrcoun/EPP622-2026-SNPCalling GitHub Wiki

WARNING: Because these reads are simulated, fastp fails to identify that they need poly-G tail trimming. Turn that on with a flag.

ADVICE: Its easier to compare the reads pre and post trimming if you run multiqc on the results of fastqc of the trimmed files (and do not use the .json files from fastp).

Analysis Summary

Based on our intial FastQC/MultiQC report: Sample C showed ~35% Illumina adapter contamination and elevated GC content. Sample D's Phred quality scores dropped below Q15 toward 3' ends

FASTP

Run fastp across all samples to ensure read quality is fixed before downstream mapping. Trim adapter sequences, remove low-quality bases, and eliminate poly-G artifacts before aligning reads.

# Navigate to fastp directory.
cd /lustre/isaac24/proj/UTK0505/alyrcoun/epp_622/test2/test2_analyses/fastp_results

# Symbolically link raw FASTQ files into the working directory
ln -s ../../raw_data_test2/*gz .

# Make a directory for trimmed results.
mkdir fastp_trimmed

If you want, look at the '.fastq.gz' files using zcat to verify their formatting without having to manually extract or uncompress it to disk first.

# zcat is a command-line utility used to view uncompressed contents of compressed (.gz) files
zcat Sample_A_R1.fastq.gz | head -n 16

Single Sample Test Run

Before looping over all files, we run a single sample first to ensure parameter settings are working.

# Run fastp on a single sample pair
# --in1 / --in2 = specifies input paired-end FASTQ files
# --out1 / --out2 = specifies output paths for cleaned, trimmed FASTQ files
# --detect_adapter_for_pe = auto-detects adapter sequences for paired-end data
# --trim_poly_g = explicitly forces poly-G tail trimming (necessary for simulated reads)
# --thread 2 = uses 2 threads for faster processing
# --json / --html = generates report files. JSON reports MUST match *.fastp.json for MultiQC to parse them
fastp \
    --in1 Sample_A_R1.fastq.gz \
    --out1 fastp_trimmed/Sample_A_R1.trimmed.fastq.gz \
    --in2 Sample_A_R2.fastq.gz \
    --out2 fastp_trimmed/Sample_A_R2.trimmed.fastq.gz \
    --detect_adapter_for_pe \
    --trim_poly_g \
    --thread 2 \
    --json Sample_A.fastp.json \
    --html Sample_A.fastp.html

JSON report: Sample_A.fastp.json HTML report: Sample_A.fastp.html Sample A did not have Illumina adapter contamination and the results reflect that.

Now let's run all the samples!

Build a run_fastp.sh script. You can build it in steps to verify paths before executing jobs.

Step 1: Testing Basename Extraction and Variable Paths Start simple by iterating over all Read 1 files. Use sed to strip file extensions. This tells us if paths and sample basenames are correct.

## Test 1 - Test loop and check basename extraction. 
# Create a 'run_fastp.sh' script 
nano run1_fastp.sh
# Nano script for Run 1:
for R1 in ../../raw_data_test2/*_R1.fastq.gz
do
    # Use sed to strip '_R1.fastq.gz' suffix from file path
    basename=$(echo "$R1" | sed 's/_R1.fastq.gz//')
    
    # Print variables to screen to confirm correct path targeting
    echo "R1 File: $R1"
    echo "Basename: $basename"
done

# Run it and see if it worked
bash run1_fastp.sh

Step 2: Dry Run Construct the entire fastp command.

## Test 2 - Print out fastp commands as a dry run to check out the generated command string and parameters before execution 

# Write the full `fastp` command inside an `echo` statement with with `xargs basename` to verify variable expansion without running the actual program
nano run2_fastp.sh
# Nano Script for Run 2:
for R1 in ../../raw_data_test2/*_R1.fastq.gz
do
    # Extract filename and strip suffix to isolate clean sample identifier
    basename=$(echo "$R1" | xargs basename | sed 's/_R1.fastq.gz//')
    echo "Targeting sample: $basename"

    # Print out constructed command string to verify variable substitution
    echo "fastp   --in1 ${basename}_R1.fastq.gz \
        --out1 fastp_trimmed/${basename}_R1.trimmed.fastq.gz \
        --in2 ${basename}_R2.fastq.gz \
        --out2 fastp_trimmed/${basename}_R2.trimmed.fastq.gz \
        --detect_adapter_for_pe \
        --trim_poly_g \
        --thread 2 \
        --json ${basename}.fastp.json \
        --html ${basename}.fastp.html"
done

# Run it and see if it worked
bash run2_fastp.sh

Step 3: Full Execution - Final Script Once variable names and paths are confirmed, remove echo statement to run fastp on all samples

# You can scrap the other run_fastp.sh scripts, but I am not since I used run1 and run2 in my example steps.
nano run_fastp.sh
# Nano script for final run
for R1 in ../../raw_data_test2/*_R1.fastq.gz
do
        # Extract sample basename
        basename=$(echo "$R1" | xargs basename | sed 's/_R1.fastq.gz//')
        echo "Processing: $R1"
        echo "Sample Name: $basename"

        # Execute fastp quality filtering and adapter/poly-G trimming
        fastp --in1 ${basename}_R1.fastq.gz \
              --out1 fastp_trimmed/${basename}_R1.trimmed.fastq.gz \
              --in2 ${basename}_R2.fastq.gz \
              --out2 fastp_trimmed/${basename}_R2.trimmed.fastq.gz \
              --detect_adapter_for_pe \
              --trim_poly_g \
              --thread 2 \
              --json ${basename}.fastp.json \
              --html ${basename}.fastp.html
done

# Did it work?
bash run_fastp.sh
# Yes it worked!

Post-trimming validation with FastQC

Compare read quality pre and post trimming by running FastQC directly on the trimmed .fastq.gz files in fastp_trimmed/ After this, we will aggregate them with MultiQC

# Create a directory for the trimmed FastQC results
mkdir -p trimmed_fastqc

# Run FastQC on all trimmed files
for trimmed_file in fastp_trimmed/*.trimmed.fastq.gz; do
    fastqc -t 2 "$trimmed_file" -o trimmed_fastqc/
done

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

# Navigate into your fastp_trimmed directory
cd fastp_trimmed
# Calculate exact reads, total bases, and average read length for FASTQ files
# 'NR%4==2' processes only the 2nd line of every 4-line FASTQ record (sequence string)
# 'r++; b+=length($0)' keeps track of total read count and adds sequence length to running base total
# column -t formats the output into organized columns
for f in *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

This gives the total reads, then the total bases, then the average read length for each file.

Results

MultiQC report path and name: /lustre/isaac24/proj/UTK0505/alyrcoun/epp_622/test2/test2_analyses/02_fastp/trimmed_fastqc/multiqc_report.html

For each sample, report R1 + R2. Use the reference genome size to calculate coverage. Formula: Coverage = total base pairs/genome size base pairs (Saccharomyces cerevisiae) genome size is ~ 12, 100,000 bp and read length is 150bp per read (paired-end)

Sample	Reads before	Reads after	Mean length before	Mean length after
A	1609300-	1609300-	150 bp-	150 bp-
B	321860-	321860-	150 bp-	150 bp-
C	1609300-	1577258-	150 bp-	125.9 bp-
D	1609300-	1609300-	150 bp-	150 bp-