BMS5021 · Introduction to Bioinformatics
Algorithms for Alignment & Read Mapping
Week 6 of Monash University BMS5021 is the algorithmic heart of Topic 2: pairwise sequence alignment by dynamic programming (global Needleman-Wunsch and local Smith-Waterman), scoring matrices and gap penalties, and fast read mapping to a reference using the Burrows-Wheeler transform (Bowtie2/HISAT2), followed by SAM/BAM post-processing with samtools. These are the steps that turn raw reads into mapped, countable data for the RNA-seq report - so the alignment maths and the mapping stats here are directly assessable.
What this chapter covers
- 01Alignment = find regions of similarity between a query and a reference, maximising matched positions
- 02Global (Needleman-Wunsch) aligns end-to-end; local (Smith-Waterman) finds the best matching sub-region; both use dynamic programming
- 03Dynamic programming (slow but OPTIMAL) vs heuristic methods (fast but approximate, e.g. BLAST, short-read aligners)
- 04DP recurrence: F(i,j) = max( F(i-1,j-1)+s(xi,yj), F(i-1,j)-d, F(i,j-1)-d ); three steps: initialise, fill, traceback
- 05Identity scoring: Id(x,y) = 1 if x=y, 0 if x!=y; alignment score = sum of aligned-pair scores
- 06Scoring matrices and gap penalties; affine gap = large opening penalty + smaller extension penalty
- 07Read mapping and indexing: hash tables / spaced seeds and the Burrows-Wheeler transform (BWT) + FM-index (Bowtie2, BWA, HISAT2, STAR)
- 08Post-processing with samtools: SAM<->BAM convert, sort THEN index (.bai), flagstat, remove PCR duplicates
Identity score and percent identity of an alignment (Week 6, has_math)
- +1Score each column with Id(x,y). Column 1 G/G = 1; column 2 A/A = 1; column 3 T/C = 0; column 4 gap/T = 0 (a gap can never match); column 5 A/A = 1; column 6 C/C = 1. The alignment identity score is the sum: 1+1+0+0+1+1 = 4.
- +1Compute percent identity. There are 4 identical positions across an alignment 6 columns long, so percent identity = 4/6 = 0.667 = about 66.7% (the gap column is part of the alignment length and counts against identity).
- +1Explain optimality. Needleman-Wunsch is a dynamic-programming method: it fills a full scoring matrix using the recurrence F(i,j) = max(diagonal + match/mismatch, up + gap, left + gap) and tracebacks the best path, so it is guaranteed to return the highest-scoring (optimal) global alignment - at the cost of speed. A heuristic aligner (e.g. BLAST or a short-read mapper) uses shortcuts like seeds to run fast on huge data, so it is approximate and can miss the true optimum.
Key terms
- Needleman-Wunsch (global alignment)
- A dynamic-programming algorithm that aligns two sequences end-to-end, guaranteeing the optimal global alignment. It initialises the first row/column with cumulative gap penalties, fills each cell F(i,j) = max(F(i-1,j-1)+s(xi,yj), F(i-1,j)-d, F(i,j-1)-d), then tracebacks.
- Smith-Waterman (local alignment)
- The dynamic-programming algorithm for the optimal LOCAL alignment - the best-matching sub-region of two sequences - using the same fill/traceback idea but resetting negative scores to zero.
- Identity score
- The simplest scoring scheme: Id(x,y) = 1 if the two aligned symbols are identical and 0 if not; the alignment score is the sum of the aligned-pair scores.
- Affine gap penalty
- A gap-scoring scheme with a large opening penalty plus a smaller extension penalty, so one long gap is penalised less than many separate short gaps - biologically realistic for a single insertion/deletion event.
- Burrows-Wheeler transform (BWT) / FM-index
- A reversible, compressible transform of a reference sequence that, with an FM-index, lets short-read aligners (Bowtie2, BWA, HISAT2) map billions of reads to the genome quickly using little memory.
- samtools sort/index
- Tools to manipulate SAM/BAM files: you sort a BAM by coordinate and THEN index it (.bai) so downstream tools can jump directly to any genomic position; flagstat summarises mapping flags.
Algorithms for Alignment & Read Mapping FAQ
How does this chapter appear in the assessment?
Week 6 is the mapping stage of the Topic 2 RNA-seq report (35%): you align trimmed reads to the reference with Bowtie2 or HISAT2, convert/sort/index with samtools, and check the alignment statistics before counting. Understanding the alignment algorithms and being able to interpret mapping stats is exactly what the report tests at this stage. Confirm the tool and dataset details on Moodle.
What is the difference between global and local alignment?
Global alignment (Needleman-Wunsch) aligns two sequences over their entire length, best when the sequences are of similar length and expected to match end-to-end. Local alignment (Smith-Waterman) finds the single best-matching sub-region, best when you only expect a shared domain or motif. Both use dynamic programming and are optimal; they differ in initialisation and in resetting negative scores to zero (local).
Why use a Burrows-Wheeler index instead of just running Needleman-Wunsch on every read?
Dynamic programming is optimal but far too slow for billions of short reads against a whole genome. The Burrows-Wheeler transform plus an FM-index lets aligners like Bowtie2 and HISAT2 find where each read maps very fast and with a small memory footprint. It is a heuristic-style speed-up that makes genome-scale mapping feasible; exact DP is reserved for aligning short pairs of sequences.
Why sort a BAM before indexing it?
Because the index lets a tool jump straight to a genomic coordinate, like using the tabs in a dictionary - and that only works if the records are in coordinate order. So the workflow is convert SAM to BAM, sort by position, then index (.bai). Skipping the sort makes the index meaningless, which is a common practical error in the Topic 2 pipeline.
Exam move
This is the most algorithmic week, so practise by hand and by tool. On paper, fill a small dynamic-programming grid for two short sequences with a match/mismatch/gap scheme and traceback the alignment, and score given alignments with the identity rule - this makes the Needleman-Wunsch recurrence F(i,j) = max(diagonal, up-gap, left-gap) concrete. Keep the one-liners straight: DP = optimal but slow, heuristic (BLAST, short-read aligners) = fast but approximate; global = end-to-end, local = best sub-region; sort THEN index. In the workshop, run the real mapping step (Bowtie2/HISAT2 then samtools) and read the alignment statistics. Ask Sia to generate fresh grids and alignments to score and to check your traceback. Confirm the Topic 2 tools and expectations on Moodle.
Working through Algorithms for Alignment & Read Mapping in BMS5021? Sia is AskSia’s AI Biology tutor — ask any BMS5021 Algorithms for Alignment & Read Mapping question and get a clear, step-by-step explanation grounded in how BMS5021 is taught and assessed. Read this chapter free, then take your hardest questions to Sia.
