3. Alignment to Reference Genome (bwa-mem2)

This section documents the execution loop for index compilation, sequence mapping against the Saccharomyces cerevisiae reference genome using bwa-mem2, and conversion into position-sorted binary alignments using samtools.

---

Workspace and coding

Follow the steps to process large-scale multi-sample alignments efficiently, with a executable script.

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

# 2. Create and open the alignment script using the vi editor
vi scripts/03_alignment.sh

# 3. Give execution permissions
chmod +x scripts/03_alignment.sh

---

Execution script code

The code below scripts/03_alignment.sh, systematically compiles indices, passes paired-end files through bwa-mem2 mem, and pipes output streams into compressed, sorted, and indexed .bam files.

#!/bin/bash

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

# Create subdirectory for alignment output
mkdir -p results/03_alignment


#Index the reference genome for BWA-MEM2
echo "Indexing reference genome for BWA-MEM2"
bwa-mem2 index data/reference/s_cerevisiae_ref.fasta
echo "Reference genome indexing completed."
echo "If you skip this step, you will get a ERROR!"


# === SAMPLE A ALIGNMENT + DIAGNOSTICS ===

echo "Processing alignment for sample A"
bwa-mem2 mem -t 4 data/reference/s_cerevisiae_ref.fasta \
    results/02_trimmed/Sample_A_R1.trimmed.fastq.gz \
    results/02_trimmed/Sample_A_R2.trimmed.fastq.gz | \
    samtools sort -@ 4 -o results/03_alignment/Sample_A.sorted.bam -

samtools index results/03_alignment/Sample_A.sorted.bam

samtools stats -@ 4 results/03_alignment/Sample_A.sorted.bam > results/03_alignment/Sample_A.stats.txt
samtools flagstat -@ 4 results/03_alignment/Sample_A.sorted.bam > results/03_alignment/Sample_A.flagstat.txt
samtools coverage results/03_alignment/Sample_A.sorted.bam > results/03_alignment/Sample_A.coverage.txt

# === SAMPLE B ALIGNMENT + DIAGNOSTICS ===
echo "Processing alignment for sample B"
bwa-mem2 mem -t 4 data/reference/s_cerevisiae_ref.fasta \
    results/02_trimmed/Sample_B_R1.trimmed.fastq.gz \
    results/02_trimmed/Sample_B_R2.trimmed.fastq.gz | \
    samtools sort -@ 4 -o results/03_alignment/Sample_B.sorted.bam -

samtools index results/03_alignment/Sample_B.sorted.bam

samtools stats -@ 4 results/03_alignment/Sample_B.sorted.bam > results/03_alignment/Sample_B.stats.txt
samtools flagstat -@ 4 results/03_alignment/Sample_B.sorted.bam > results/03_alignment/Sample_B.flagstat.txt
samtools coverage results/03_alignment/Sample_B.sorted.bam > results/03_alignment/Sample_B.coverage.txt

# === SAMPLE C ALIGNMENT + DIAGNOSTICS ===
echo "Processing alignment for sample C"
bwa-mem2 mem -t 4 data/reference/s_cerevisiae_ref.fasta \
    results/02_trimmed/Sample_C_R1.trimmed.fastq.gz \
    results/02_trimmed/Sample_C_R2.trimmed.fastq.gz | \
    samtools sort -@ 4 -o results/03_alignment/Sample_C.sorted.bam -

samtools index results/03_alignment/Sample_C.sorted.bam

samtools stats -@ 4 results/03_alignment/Sample_C.sorted.bam > results/03_alignment/Sample_C.stats.txt
samtools flagstat -@ 4 results/03_alignment/Sample_C.sorted.bam > results/03_alignment/Sample_C.flagstat.txt
samtools coverage results/03_alignment/Sample_C.sorted.bam > results/03_alignment/Sample_C.coverage.txt

# === SAMPLE D ALIGNMENT + DIAGNOSTICS ===
echo "Processing alignment for sample D"
bwa-mem2 mem -t 4 data/reference/s_cerevisiae_ref.fasta \
    results/02_trimmed/Sample_D_R1.trimmed.fastq.gz \
    results/02_trimmed/Sample_D_R2.trimmed.fastq.gz | \
    samtools sort -@ 4 -o results/03_alignment/Sample_D.sorted.bam -

samtools index results/03_alignment/Sample_D.sorted.bam

samtools stats -@ 4 results/03_alignment/Sample_D.sorted.bam > results/03_alignment/Sample_D.stats.txt
samtools flagstat -@ 4 results/03_alignment/Sample_D.sorted.bam > results/03_alignment/Sample_D.flagstat.txt
samtools coverage results/03_alignment/Sample_D.sorted.bam > results/03_alignment/Sample_D.coverage.txt

# === Get report with MULTIQC ===
echo "Compiling final alignment metrics report with MultiQC"
multiqc results/03_alignment/ -o results/ -n 03_multiqc_alignment_report

echo "Genomic alignment pipeline completed successfully."

---

Usage: 03_alignment.sh

# 1. Activate the genomic Conda environment
conda activate snps

# 2. Execute the alignment pipeline via background deployment
nohup ./scripts/03_alignment.sh > results/03_alignment.log 2>&1 &

# 3. See thereal-time progress through log file
tail -f results/03_alignment.log

---

Results

• MultiQC report path and name: results/03_multiqc_alignment_report.html

Because the Saccharomyces cerevisiae genome is structurally divided into 16 separate chromosomes, the outputs within the *.coverage.txt are written row-by-row for each chromosome. Next step is sum all aligned segments and average the molecular read depth.

# Get the properly paired alignment
grep -E "mapped \(|properly paired \(" results/03_alignment/Sample_A.flagstat.txt
grep -E "mapped \(|properly paired \(" results/03_alignment/Sample_B.flagstat.txt
grep -E "mapped \(|properly paired \(" results/03_alignment/Sample_C.flagstat.txt
grep -E "mapped \(|properly paired \(" results/03_alignment/Sample_D.flagstat.txt

# Average mean genomic depth:
awk 'NR>1 {sum+=$7} END {print "Sample A avg depth: ", sum/16}' results/03_alignment/Sample_A.coverage.txt
awk 'NR>1 {sum+=$7} END {print "Sample B avg depth: ", sum/16}' results/03_alignment/Sample_B.coverage.txt
awk 'NR>1 {sum+=$7} END {print "Sample C avg depth: ", sum/16}' results/03_alignment/Sample_C.coverage.txt
awk 'NR>1 {sum+=$7} END {print "Sample D avg depth: ", sum/16}' results/03_alignment/Sample_D.coverage.txt

What's doing the code?

1. grep -E "mapped \(|properly paired \(" ...Sample_A.flagstat.txt This code use grep to scan the alignment diagnostic file to extract the results. Add the option -E to flag enables Extended Regular Expressions and where R1 and R2 successfully synchronized in the correct RF orientation.

2. awk 'NR>1 {sum+=\$7} END {print "Sample A avg depth: ", sum/16}' ...Sample_A.coverage.txt This code calculates the mean depth across the whole genome. The actions of the code are: 1) ignore the first row, which correspond to the header; 2) Target the column seven and sum all the values that belong to all 16 chromosomes; 3) Finally, divides the total by the 16 chromosomes of S. cerevisiae.

Alignment summary

Sample	% mapped	% properly paired	Mean depth
A	100	100	19.99
B	100	100	4.00
C	100	100	16.44
D	100	100	19.13

The results show that every read processed through the mapping pipeline resolved with a high overall alignment rate. Additionally, the data indicated that filtering out poly-G sequences and adapters successfully preserved paired-end coordination. Finally, samples A and D achieved deep ~20X analytical coverage targets, whereas sample B presented coverage constraints (~4X depth).