Read Alignment with BWA-MEM

Before you start

You've done Trimming and Filtering Reads, so you have clean FASTQ reads to align.
Comfortable running a few command-line lines. A free Google account to run the notebook.
Any term new? The Glossary has it.

Learning objectives

By the end of this lesson you will be able to: explain what read alignment does and why a reference index is needed, build an index and align reads with a tool like BWA-MEM, read the basic fields of the resulting SAM file, and run a quick sanity check on how many reads mapped.

What alignment actually does

A sequencer hands you millions of short reads, but it does not tell you where in the genome each one came from. Alignment (also called mapping) answers exactly that question: for every read, it finds the spot on a known reference genome that the read matches best, and records that position.

Think of the reference as the picture on a jigsaw box and your reads as the loose pieces. Alignment is the act of laying each piece onto the spot it fits. Once every read has a coordinate, you can stack them up, see how many reads cover each position, and spot the places where your sample differs from the reference. Those differences are variants, and finding them is what the rest of this track builds toward.

Decode the jargonReference genome

A reference genome is a finished, agreed-upon sequence for a species (for humans, the GRCh38 assembly) that the community uses as a common coordinate system. It is not a "perfect" or "correct" genome, just a shared map so everyone's positions mean the same thing. You align your reads against it to describe where they sit and how they differ.

Decode the jargonBWA-MEM

BWA (Burrows-Wheeler Aligner) is a widely used program for mapping reads to a reference, and MEM is its recommended algorithm for reads longer than about 70 bases. You'll also hear of Bowtie2 (a close alternative for DNA), HISAT2 and STAR (used for RNA-seq, because they handle reads that span across spliced introns). For DNA mapping, BWA-MEM is the dependable default.

▶ Run it in Google Colab

Colab is a free Linux machine in your browser, so BWA runs with nothing to install on your computer.

Open the notebook in Colab →

1Install the tools and make a tiny reference

We install BWA and samtools (we'll use samtools to convert and inspect the result). The notebook then writes a small reference genome and a handful of reads taken from it, so the whole thing runs in seconds.

!apt-get -qq install -y bwa samtools
!bwa 2>&1 | head -3

2Index the reference

Before BWA can search the reference, it builds an index, a set of helper files that let it look up sequences fast (the same idea as the index at the back of a book). You do this once per reference.

!bwa index reference.fasta

Decode the jargonIndex

An index is a pre-built lookup structure. Searching the raw genome letter by letter for every read would be hopelessly slow, so the aligner first reorganizes the reference into a form it can query in an instant. Indexing takes a little time up front and is reused for every read afterward.

3Align the reads

Now the main event. bwa mem takes the indexed reference and your reads and prints alignments in SAM format. We pipe that straight into samtools to save it as a compact BAM file.

!bwa mem reference.fasta reads.trimmed.fastq > aligned.sam
# or, save space by going straight to BAM:
!bwa mem reference.fasta reads.trimmed.fastq | samtools view -b -o aligned.bam

Reading the SAM file

SAM (Sequence Alignment/Map) is just a text table, one line per read. After a header (lines starting with @), each alignment line has eleven required columns. You rarely read them by hand, but knowing the key ones demystifies everything downstream.

Column	Name	What it tells you
1	`QNAME`	The read's name.
2	`FLAG`	A bundle of yes/no facts about the read (mapped? reverse strand? paired?).
3	`RNAME`	Which reference sequence (chromosome) it mapped to.
4	`POS`	The position on that chromosome where the alignment starts.
5	`MAPQ`	Mapping quality: how confident the aligner is about this position.
6	`CIGAR`	A compact code for how the read matches (matches, insertions, deletions, clips).
10	`SEQ`	The read's actual bases.
11	`QUAL`	The per-base quality scores.

Decode the jargonBAM vs SAM

SAM is the human-readable text version. BAM is the exact same information compressed into a binary file, much smaller and faster for programs to read. In real projects you keep BAM on disk and only convert to SAM when you want to peek inside. Tools read BAM directly.

Decode the jargonMAPQ (mapping quality)

MAPQ is a score for how sure the aligner is that a read landed in the right place. A high MAPQ means the read fits one spot clearly. A low MAPQ (near 0) means the read could fit several places about equally well, often because it falls in a repetitive region. Variant callers lean on MAPQ to decide which reads to trust.

4Did it map? A quick sanity check

Before moving on, confirm that your reads actually mapped. samtools flagstat prints a one-glance summary, including the all-important percentage of mapped reads.

!samtools flagstat aligned.bam

For a clean sample against the right reference, you expect a high mapping rate (often well above 90%). A low rate is a red flag: wrong reference, leftover adapters, or contaminated reads.

⚠️ Align to the right reference

Reads only map well to the genome they actually came from. Aligning human reads to a mouse genome, or using the wrong build, produces a low mapping rate and garbage downstream. Always match your reference to your organism (and note which build you used, so your work is reproducible).

Where the mapped reads go next

You now have a BAM file: every read, with a coordinate. But it's in the order the reads came off the machine, which is not useful yet. The next lesson is the samtools survival kit, where you sort and index the BAM so tools can jump to any region instantly, the step that turns a raw alignment into something you can actually analyze.

🔶 Level up: read groups and paired-end

Real pipelines align paired-end reads (two files) and tag each BAM with a read group (-R) that records the sample and sequencing run. Variant callers like GATK require read groups. The full command grows to something like bwa mem -R '@RG\tID:run1\tSM:sample1' ref.fasta R1.fastq R2.fastq, but the core idea is exactly what you just did.

Check your understanding

Why do you index the reference before aligning?

Searching the raw genome for every one of millions of reads would be far too slow. The index reorganizes the reference into a lookup structure the aligner can query almost instantly. You build it once and reuse it for all your reads.

What does the POS column in a SAM file tell you?

It's the position on the reference (the coordinate on a given chromosome) where the read's alignment begins. Combined with RNAME (which chromosome), it pins the read to an exact spot on the genome.

Your flagstat shows only 40% of reads mapped. What might be wrong?

A low mapping rate usually means a mismatch between reads and reference: the wrong species or genome build, heavy adapter or quality contamination that trimming missed, or a contaminated sample. It's worth investigating before trusting anything downstream.

What does a low MAPQ (close to 0) usually indicate about a read?

MAPQ measures how sure the aligner is about a read's position. A value near 0 means the read maps about equally well to several spots, commonly because it falls in a repetitive region, so callers trust it less.

Which aligners are designed for RNA-seq reads that span across spliced introns?

HISAT2 and STAR are splice-aware aligners built for RNA-seq, where reads can cross intron boundaries. For plain DNA mapping, BWA-MEM and Bowtie2 are the usual choices.

Next in Track 2

Read Alignment: Putting Reads Back on the Genome