Module 2.3 bwa mem2

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 THE FILES FROM 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_fastp/*trimmed.fastq.gz .

Now lets map one read set to see how it works.

bwa-mem2 mem \
Fraxinus_pennsylvanica_PE-57.fasta \
 21_FP_OH_FSB_F_10_R1.trimmed.fastq.gz \
 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. (this takes about 7 minutes)

bwa-mem2 mem \
-t 2 \
Fraxinus_pennsylvanica_PE-57.fasta \
21_FP_OH_FSB_F_10_R1.trimmed.fastq.gz \
21_FP_OH_FSB_F_10_R2.trimmed.fastq.gz \
| \
samtools sort \
-@ 2 \
-o 21_FP_OH_FSB_F_10.sorted.bam \
-

Build an index (.bai)

samtools index 21_FP_OH_FSB_F_10.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 got from AI. I normally use sed, this seems easier).

# Loop over every trimmed R1 file.
for R1 in ../02_fastp/*_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"

    bwa-mem2 mem \
        Fraxinus_pennsylvanica_PE-57.fasta \
        -t 2 \
        $R1 \
        $R2 \
      | 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"

# Where the finished BAM files should go.
OUTDIR="."

# 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"
    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"

Stats

We need to know if our mapping worked. We are going to use two different samtools commands then aggregate with multiqc. Put this in a script and run it

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

Take a look at the three types of output files. We can aggregate them with multiqc.

multiqc .