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

Commands

mkdir bwamem2_3
cd bwamem2_3/

ln -s /lustre/isaac24/proj/UTK0505/rsalani/exam2/rawdata/ref/* .
ln -s ../fastp2_withWARNING/*trimmed* .

#make a bash script for bwamem
nano loop_for_bwamems.sh

------------------------------------------------------------------------------
# Loop over every trimmed R1 file in the current directory
for R1 in *_R1.trimmed.fastq.gz
do
    # Get matching R2 filename
    R2="${R1/_R1/_R2}"

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

    echo "Aligning sample: $SAMPLE"
    echo "R1: $R1"
    echo "R2: $R2"

    # Check that the matching R2 file exists
    if [ ! -f "$R2" ](/rsalani13/EPP622-2026-SNPCalling/wiki/-!--f-"$R2"-); then
        echo "ERROR: Missing R2 file for $SAMPLE"
        continue
    fi

    # Align reads, sort alignments, and output BAM
    bwa-mem2 mem \
        -t 2 \
        sacCer3.fa \
        "$R1" \
        "$R2" \
      | samtools sort \
        -@ 2 \
        -o "${SAMPLE}.sorted.bam" \
        -

    # Stop if alignment or sorting failed
    if [ ${PIPESTATUS[0]} -ne 0 ](/rsalani13/EPP622-2026-SNPCalling/wiki/-${PIPESTATUS[0]}--ne-0-); then
        echo "ERROR: Alignment or sorting failed for $SAMPLE"
        continue
    fi

    # Index the sorted BAM
    samtools index "${SAMPLE}.sorted.bam"

    echo "Finished: $SAMPLE"
done
------------------------------------------------------------------------------
#run it
bash loop_for_bwamems.sh

#bash script for samtools
nano did_bwamem_work.sh
----------------------------------------------------------------------------------------------
for bamfile in *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
----------------------------------------------------------------------------------------------
nano did_bwamem_work.sh

#run multiqc again!
multiqc .

#rename the multiqc again!
mv multiqc_report.html multiqc_report_for_bwamem.html

Results

MultiQC report path and name: /lustre/isaac24/proj/UTK0505/rsalani/exam2/analyses/bwamem2_3/multiqc_report_for_bwamem.html

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