Step 3 Alignment (bwa‐mem2 samtools) - alyrcoun/EPP622-2026-SNPCalling GitHub Wiki

Now that our reads are trimmed and verified we will align them to the Saccharomyces cerevisiae reference genome

Analysis Overview

Use bwa-mem2 to align our trimmed FASTQ reads to the reference genome. After alignment, use samtools to convert the alignments into sorted BAM files and index them,

Data linking

# Navigate into the bwamem2_results directory
cd bwamem2_results

# Link to reference genome FASTA and index files
ln -s /lustre/isaac24/proj/UTK0505/test2/ref/sacCer3.fa* .

# Link to trimmed FASTQ files from the fastp_results directory
 ln -s ../fastp_results/fastp_trimmed/*trimmed.fastq.gz .

# Check that the reference genome and trimmed FASTQ files are in the directory
ls

Shell Script

Write a shell script with built-in safety checks to map all four trimmed read sets automatically

nano run_bwa.sh
# Nano script:

# Loop over every trimmed R1 file.
for R1 in ../fastp_results/fastp_trimmed/*_R1.trimmed.fastq.gz
do
    # Get R2 name by switching R1 with R2, retaining the full path
    R2="${R1/_R1/_R2}"

    # Find the short sample name - no path
    # first, extract the filename from a full file path.
    # The ## tells Bash to find the longest possible match from
    # the beginning of the string and delete it.
    # The */ matches any text ending with a forward slash.
    SAMPLE="${R1##*/}"
    # second, remove _R1.trimmed.fastq.gz from the file name
    SAMPLE="${SAMPLE%_R1.trimmed.fastq.gz}"

    echo "$R1"
    echo "$R2"
    echo "$SAMPLE"

    # Map reads using bwa-mem2, add Read Group tags, pipe to samtools sort
    # sacCer3.fa is the reference genome file to align reads against
    # -t 2 allocates 2 multi-threading threads to speed up alignment
    # -R stamps sample name into BAM header so bcftools can separate samples
    # samtools sort -@ 2 -o sorts alignments by coordinate and sorted BAM file 
    bwa-mem2 mem \
        sacCer3.fa \
        -t 2 \
        -R "@RG\tID:${SAMPLE}\tSM:${SAMPLE}\tPL:ILLUMINA" \
        "$R1" \
        "$R2" \
      | samtools sort \
        -@ 2 \
        -o "${SAMPLE}.sorted.bam" \
        -

    ## Build an index (.bai) so  IGV/JBrowse can jump to regions
    samtools index "${SAMPLE}.sorted.bam"

done

echo "All done mapping!"

Execute the script inside a screen session to ensure it continues if the SSH session drops

# Start a detached screen session names bwa
screen -S bwa
# Run the script
bash run_bwa.sh
# Detach from screen session: Press Ctrl + a then d
# Re-attach to check progress later
screen -r bwa

Did our mapping work? Alignment Statistics and Quality Check

We need to know if our mapping worked. We are going to use samtools stats and samtools flagstat reports for each BAM file. Then aggregrate them using multiqc. Put this in a script and run it:

nano stats.sh
# Nano script
for bamfile in *.sorted.bam
do
        # The ${VAR%pattern} syntax removes "pattern" from the END of the value.
        SAMPLE=${bamfile%.sorted.bam}
        echo $SAMPLE

        samtools stats -@ 2 $bamfile > $SAMPLE.stats.txt
        samtools flagstat -@ 2 $bamfile > $SAMPLE.flagstat.txt
        samtools coverage $bamfile > $SAMPLE.coverage.txt

done

# Aggregate alignment results with MultiQC
multiqc .

Results

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

Sample	% mapped	% properly paired	Mean depth
A	100-	100-	20x-
B	100-	100-	4x-
C	100-	100-	16.4x-
D	100-	100-	19.1x-