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

Step 3: Alignment (bwa-mem2 + samtools)

Commands

BWA-MEM2 was used to align the trimmed paired end sequencing reads to the reference genome. Samtools was then used to sort, index, and evaluate the resulting BAM files. MultiQC combined the alignment statistics to compare mapping performance and sequencing depth across samples.

Created the alignment directory and linked the reference genome and trimmed reads:

mkdir 03_bwamem2
cd 03_bwamem2
ln -s /lustre/isaac24/proj/UTK0505/aduckwo4/622_test2/ref
ln -s ../02_fastp/fastp_trimmed/*trimmed.fastq.gz .
nano run_bwa.sh

Created run_bwa.sh:

for R1 in *_R1.trimmed.fastq.gz
do
    R2="${R1/_R1/_R2}"

    SAMPLE="${R1##*/}"
    SAMPLE="${SAMPLE%_R1.trimmed.fastq.gz}"

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

    bwa-mem2 mem \
        -t 2 \
        ref/sacCer3.fa \
        "$R1" \
        "$R2" \
    | samtools sort \
        -@ 2 \
        -o "$SAMPLE.sorted.bam" \
        -

    samtools index "$SAMPLE.sorted.bam"

done

# R2: identifies the matching reverse read file by replacing _R1 with _R2.
# SAMPLE: extracts the sample name by removing the file path and R1 suffix.
# echo: prints the filenames and sample name.
# bwa-mem2 mem: aligns paired end sequencing reads to the reference genome.
# -t 2: uses two threads for alignment.
# ref/sacCer3.fa: specifies the reference genome.
# "$R1" and "$R2": specify the paired end input files.
# |: pipes alignment output directly into samtools sort.
# samtools sort: sorts alignments by genomic position.
# -@ 2: uses two additional threads for sorting.
# -o: specifies the output sorted BAM filename.
# -: reads alignment data from standard input.
# samtools index: creates a BAM index

Ran the alignment script and created run_sam.sh:

bash run_bwa.sh
nano run_sam.sh

Created run_sam.sh:

for bamfile in *.sorted.bam
do
    # Extract sample name
    SAMPLE="${bamfile%.sorted.bam}"

    echo "$SAMPLE"

    # Generate alignment statistics
    samtools stats -@ 2 "$bamfile" > "$SAMPLE.stats.txt"

    samtools flagstat -@ 2 "$bamfile" > "$SAMPLE.flagstat.txt"

    samtools coverage "$bamfile" > "$SAMPLE.coverage.txt"

done
# SAMPLE: removes the .sorted.bam suffix to identify the sample.
# echo: prints the current sample name.
# samtools stats: generates detailed alignment statistics.
# -@ 2: uses two additional threads.
# samtools flagstat: summarizes mapped, paired, and properly paired reads.
# samtools coverage: calculates reference coverage and mean sequencing depth.
# >: redirects command output into the specified text file.

Ran Samtools and generated the MultiQC report:

bash run_sam.sh
multiqc .
cp multiqc_report.html multiqc_bam_report.html
scp [email protected]:/lustre/isaac24/proj/UTK0505/aduckwo4/622_test2/03_bwamem2/multiqc_bam_report.html ~/Downloads/
# bash: executes the Samtools statistics script.
# multiqc .: combines alignment statistics into one report.
# cp: copies the MultiQC report to a new filename.

Results

MultiQC report path and name: /lustre/isaac24/proj/UTK0505/aduckwo4/622_test2/03_bwamem2/multiqc_bam_report.html

Sample	% mapped	% properly paired	Mean depth
A	100%	100%	20.0 x
B	100%	100%	4.0 x
C	100%	100%	16.4 x
D	100%	100%	19.1 x