Module 2.2 fastp, bwa mem, samtools - mestato/EPP622_2026 GitHub Wiki
Starting on sphinx from this directory: /lustre/isaac24/proj/UTK0505/yourusesrid with conda using the snps environment.
Fastp
Get set up. make a new analysis folder for this and symbolically link the input files
cd analyses
mkdir 02_fastp
cd 02_fastp
ln -s ../../raw_data/*gz .
mkdir fastp_trimmed
Lets run one sample first. We can generate two types of output - html and json. If we want to use MultiQC later, we need to name the json output files matching the pattern *.fastp.json
fastp --in1 21_FP_OH_FSB_F_10_R1.fastq.gz \
--out1 21_FP_OH_FSB_F_10_R1.trimmed.fastq.gz \
--in2 21_FP_OH_FSB_F_10_R2.fastq.gz \
--out2 21_FP_OH_FSB_F_10_R2.trimmed.fastq.gz \
--detect_adapter_for_pe \
--thread 2 \
--json 21_FP_OH_FSB_F_10_R1.fastp.json \
--html 21_FP_OH_FSB_F_10_R1.fastp.html
It autodetected our adapter which is convenient. Secure copy the html to your laptop to check out the results.
Now lets run all the samples. Create a run_fastp.sh file. this took about a minute per file for me.
Here is a script that I wrote. I start simple:
## Run all files through fastp
for R1 in ../../raw_data/*_R1.fastq.gz
do
basename=$(echo "$R1" | sed 's/_R1.fastq.gz//')
echo $R1
echo $basename
done
Now start add the fastp part and check out those variable names
## Run all files through fastp
for R1 in ../../raw_data/*_R1.fastq.gz
do
basename=$(echo "$R1" | sed 's/_R1.fastq.gz//')
echo $R1
echo $basename
echo "fastp --in1 ${basename}_R1.fastq.gz \
--out1 ${basename}_R1.trimmed.fastq.gz \
--in2 ${basename}_R2.fastq.gz \
--out2 ${basename}_R2.trimmed.fastq.gz \
--detect_adapter_for_pe \
--thread 2 \
--json ${basename}_R1.fastp.json \
--html ${basename}_R1.fastp.html
"
done
Now remove the echo and see if it works
## Run all files through fastp
for R1 in ../../raw_data/*_R1.fastq.gz
do
basename=$(echo "$R1" | sed 's/_R1.fastq.gz//')
echo $R1
echo $basename
fastp --in1 ${basename}_R1.fastq.gz \
--out1 ${basename}_R1.trimmed.fastq.gz \
--in2 ${basename}_R2.fastq.gz \
--out2 ${basename}_R2.trimmed.fastq.gz \
--detect_adapter_for_pe \
--thread 2 \
--json ${basename}_R1.fastp.json \
--html ${basename}_R1.fastp.html
done
But if we let AI loose we get a lot more complexity...
#!/bin/bash
# run_fastp.sh
# Run fastp on every pair of paired-end FASTQ files in this folder.
# Usage: bash run_fastp.sh
# Stop the script immediately if any command fails.
set -e
# Make a folder to hold the fastp reports (the -p means "don't complain if it already exists").
mkdir -p fastp_reports
mkdir -p trimmed
# Loop over every R1 file in the current directory.
# The * is a wildcard: it matches ANY text, so this matches all 10 R1 files.
for R1 in *_R1.fastq.gz
do
# Build the matching R2 filename by swapping "_R1" for "_R2" in the name.
# The ${VAR/old/new} syntax replaces the first match of "old" with "new".
R2="${R1/_R1/_R2}"
# Build a short sample name by chopping "_R1.fastq.gz" off the end.
# The ${VAR%pattern} syntax removes "pattern" from the END of the value.
# Example: 21_FP_OH_FSB_F_10_R1.fastq.gz -> 21_FP_OH_FSB_F_10
SAMPLE="${R1%_R1.fastq.gz}"
# Print a message so we can watch the progress in the terminal.
echo "Now processing sample: $SAMPLE"
# Safety check: make sure the R2 file actually exists before running fastp.
if [ ! -f "$R2" ]; then
echo " WARNING: no R2 file found for $R1 -- skipping."
continue # "continue" jumps to the next sample in the loop
fi
# This is the fastp command, with the filenames replaced by variables.
fastp \
--in1 "$R1" \
--out1 trimmed/"$SAMPLE"_R1.trimmed.fastq.gz \
--in2 "$R2" \
--out2 trimmed/"$SAMPLE"_R2.trimmed.fastq.gz \
--detect_adapter_for_pe \
--thread 2 \
--json fastp_reports/"$SAMPLE".fastp.json \
--html fastp_reports/"$SAMPLE".fastp.html
done
echo "All done! Reports are in the fastp_reports/ folder."
What are the positives and negatives of using AI? How could you make sure you understand what it is doing and if it is doing what you expect?
Screen
What if you want this to run and you want to walk away and return later? We can use screen for that. GNU Screen is a terminal multiplexer that lets you run a bash shell or other commands and keep them active in the background even if you disconnect.
Basic Screen Commands
- Start a new session: Type
screen. - Start a named session: Type
screen -S myname - List active sessions: Type
screen -ls - Reattach to a session: Type
screen -rorscreen -r myname - Detach from inside a session: Press Ctrl + A then D.
Example for our bash script above.
screen -S fastp
bash run_fastq.sh
Ctrl + A D
Now you can exit. Come back later to see if the job finished or had any errors.
screen -r fastp
For Homework:
- Run multiqc on the output json files from fastp and see what it says.
- Rerun fastqc and multiqc on the trimmed files. Did fastp fix the reports? Doing this is best practice to make sure fastq fixed anything that needed fixing.
BWA-MEM2
Next we will run bwa-mem2 to map the trimmed files to the reference genome. From analyses dir
mkdir 03_bwamem2
cd 03_bwamem2
If you were starting from scratch, you would index the reference genome. I have already done this for you, so you don't have to do it and wait on it (it took about 5 minutes). Here is what i did:
YOU DO NOT HAVE TO DO THIS
ln -s /lustre/isaac24/proj/UTK0505/module2_data/reference_genome/Fraxinus_pennsylvanica_PE-57.fasta .
bwa-mem2 index Fraxinus_pennsylvanica_PE-57.fasta
INSTEAD YOU CAN LINK TO THE FILES HERE
ln -s /lustre/isaac24/proj/UTK0505/module2_data/reference_genome/Fraxinus_pennsylvanica_PE-57.fasta* .
We also need our trimmed files.
ln -s ../02_fastq/trimmed/*trimmed.fastq.gz .
Now lets map one read set to see how it works.
bwa-mem2 mem \
Fraxinus_pennsylvanica_PE-57.fasta \
../02_fastp/trimmed/21_FP_OH_FSB_F_10_R1.trimmed.fastq.gz \
../02_fastp/trimmed/21_FP_OH_FSB_F_10_R2.trimmed.fastq.gz \
> 21_FP_OH_FSB_F_10.sam
We are doing this so we can look at the sam (do a head and a tail) but it is actually a very bad idea.
Now, lets see how to go straight to a sorted bam.
bwa-mem2 mem \
-t 2 \
-R "@RG\tID:21_FP_OH_FSB_F_10\tSM:21_FP_OH_FSB_F_10\tPL:ILLUMINA"
Fraxinus_pennsylvanica_PE-57.fasta \
../02_fastp/trimmed/21_FP_OH_FSB_F_10_R1.trimmed.fastq.gz \
../02_fastp/trimmed/21_FP_OH_FSB_F_10_R2.trimmed.fastq.gz \
| samtools sort \
-@ 2 \
-o 21_FP_OH_FSB_F_10_R1.sorted.bam \
-
# Build an index (.bai) so genome browsers like IGV can jump to a region
# without reading the whole file.
samtools index 21_FP_OH_FSB_F_10_R1.sorted.bam
How much smaller is the bam than the sam?
Now we can get a very similar script going as we did with fastp, to utilize a for loop. You'll also want to run this in screen.
Here is a script I wrote (other than the sample name extraction, which I stole from AI. I normally use sed, this seems easier).
# Loop over every trimmed R1 file.
for R1 in ../02_fastp/trimmed/*_R1.trimmed.fastq.gz
do
# Find the short sample name.
SAMPLE="${R1##*/}"
SAMPLE="${SAMPLE%_R1.trimmed.fastq.gz}"
echo "$SAMPLE"
bwa-mem2 mem \
Fraxinus_pennsylvanica_PE-57.fasta -t 2 \
"$R1" \
${SAMPLE}_R2.trimmed.fastq.gz \
| samtools sort \
-@ 2 \
-o ${SAMPLE}.sorted.bam" \
-
# Build an index
samtools index "${SAMPLE}.sorted.bam"
done
Here is the much longer AI version. What is being added?
#!/bin/bash
# run_bwa.sh
# Map every pair of trimmed FASTQ files to the green ash reference genome.
# Usage: bash run_bwa.sh
# ----------------------------------------------------------------------
# SETTINGS - these are the only lines you should need to change.
# Putting them at the top (instead of burying them in the loop) means
# you edit ONE place instead of hunting through the whole script.
# ----------------------------------------------------------------------
# The reference genome. bwa-mem2 looks for the index files that sit NEXT to
# this file (they end in .0123, .amb, .ann, .bwt.2bit.64, .pac).
# You give it the name of the FASTA; it finds the rest on its own.
REF="Fraxinus_pennsylvanica_PE-57.fasta"
# Where the trimmed reads from the fastp step live.
INDIR="../02_fastp/trimmed"
# Where the finished BAM files should go.
OUTDIR="mapped"
# How many CPU threads to use. Ask for no more than your computer/node has.
THREADS=2
# ----------------------------------------------------------------------
# SAFETY CHECKS - fail early with a clear message instead of failing
# confusingly halfway through 10 samples.
# ----------------------------------------------------------------------
# Does the reference FASTA exist?
if [ ! -f "$REF" ]; then
echo "ERROR: cannot find the reference file: $REF"
exit 1
fi
# Does the index exist? We test for one of the index files. The .bwt.2bit.64
# file is the big one, so if it is there the indexing step finished.
if [ ! -f "${REF}.bwt.2bit.64" ]; then
echo "ERROR: cannot find the bwa-mem2 index for $REF"
echo " Expected a file named ${REF}.bwt.2bit.64"
echo " Do NOT re-index - ask your instructor for the correct path."
exit 1
fi
# Make the output folder (-p means "don't complain if it already exists").
mkdir -p "$OUTDIR"
# ----------------------------------------------------------------------
# THE LOOP
# ----------------------------------------------------------------------
# Loop over every trimmed R1 file. The * is a wildcard: it matches ANY text,
# so this picks up all 10 samples without you typing any filenames.
for R1 in "$INDIR"/*_R1.trimmed.fastq.gz
do
# Build the matching R2 filename by swapping "_R1" for "_R2".
# ${VAR/old/new} replaces the FIRST match of "old" with "new".
R2="${R1/_R1/_R2}"
# Build a short sample name.
# Step 1: ${VAR##*/} removes everything up to and including the last "/",
# turning "trimmed/21_FP_OH_FSB_F_10_R1.trimmed.fastq.gz"
# into "21_FP_OH_FSB_F_10_R1.trimmed.fastq.gz"
# Step 2: ${VAR%pattern} removes "pattern" from the END,
# leaving "21_FP_OH_FSB_F_10"
SAMPLE="${R1##*/}"
SAMPLE="${SAMPLE%_R1.trimmed.fastq.gz}"
echo "Now mapping sample: $SAMPLE"
# Make sure the R2 file actually exists before we start.
if [ ! -f "$R2" ]; then
echo " WARNING: no R2 file found for $R1 -- skipping."
continue # "continue" jumps to the next sample in the loop
fi
# The read group (-R) stamps the sample name INSIDE the BAM file.
# Without it, later tools cannot tell which reads came from which sample.
# \t means "tab" - the pieces have to be separated by tabs.
# ID = a unique name for this batch of reads
# SM = the sample name (this is the one downstream tools care about)
# PL = the sequencing platform
RG="@RG\tID:${SAMPLE}\tSM:${SAMPLE}\tPL:ILLUMINA"
# Now the actual work. The "|" character is a PIPE: it hands the output
# of one command straight to the next one without writing a file to disk.
# We do this because the raw SAM output is enormous; piping it keeps us
# from ever storing it.
#
# bwa-mem2 mem ... aligns the reads, prints SAM to the screen (stdout)
# samtools sort takes that SAM and writes a coordinate-sorted BAM
#
# The "-" in the samtools command means "read from the pipe, not a file".
# The backslash at the end of a line means "this command continues below".
bwa-mem2 mem \
-t "$THREADS" \
-R "$RG" \
"$REF" \
"$R1" \
"$R2" \
| samtools sort \
-@ $THREADS \
-o "$OUTDIR/${SAMPLE}.sorted.bam" \
-
# Build an index (.bai) so genome browsers like IGV can jump to a region
# without reading the whole file.
samtools index "$OUTDIR/${SAMPLE}.sorted.bam"
echo " Finished $SAMPLE"
done
echo ""
echo "All done! Sorted BAM files are in the $OUTDIR/ folder."
echo ""
echo "Quick check - how well did the reads map?"
echo "Run this on any one sample:"
echo " samtools flagstat $OUTDIR/21_FP_OH_FSB_F_10.sorted.bam"