Mapping Reads to a Reference Genome

A DNA sequencer does not read a genome from start to finish. Instead it produces millions of short strings of letters called reads, each usually 50 to 150 bases long, and it does not tell you where in the genome each read came from. Read mapping, also called alignment, is the step that figures this out: it places each read at its most likely position on a reference genome. A reference genome is an already-assembled sequence for the species, stored as plain text in a FASTA file, that works like the picture on a jigsaw puzzle box you match your pieces against.

Common Aligners

The program that does the placing is called an aligner. For short Illumina reads the two most common choices are BWA-MEM and Bowtie2; both are fast and accurate for reads up to a few hundred bases. For long, error-prone reads from Oxford Nanopore or PacBio machines, minimap2 is the standard. In this lesson we use BWA-MEM, but all three take the same inputs, a reference genome and your reads, and produce the same kind of output.

# Build an index of the reference so the aligner can search it quickly.
bwa index reference.fa

# Align paired-end reads to the reference; the alignment goes to aln.sam.
bwa mem reference.fa reads_R1.fastq.gz reads_R2.fastq.gz > aln.sam

[M::bwa_idx_load_from_disk] read 0 ALT contigs
[M::process] read 20000 sequences (1000000 bp)...
[M::mem_pestat] analyzing insert size distribution for orientation FR...
[M::mem_pestat] (25, 50, 75) percentile: (297, 342, 397)
[M::mem_process_seqs] Processed 20000 reads in 1.982 CPU sec, 2.030 real sec
[main] Version: 0.7.17-r1188
[main] CMD: bwa mem reference.fa reads_R1.fastq.gz reads_R2.fastq.gz
[main] Real time: 2.590 sec; CPU: 2.431 sec

SAM And BAM Files

The aligner writes its results in SAM format (Sequence Alignment/Map), a plain-text file with one line per read describing where and how it aligned. SAM files are human-readable but large, so they are almost always converted to BAM, the compressed binary version of the same data. The everyday tool for handling them is samtools: it converts SAM to BAM, sorts the alignments by their position along the genome, and builds an index so downstream tools can jump straight to any region without scanning the whole file.

# Convert the text SAM file to compressed binary BAM (-b means output BAM).
samtools view -b aln.sam > aln.bam

# Sort the alignments by their coordinate on the reference genome.
samtools sort aln.bam -o aln.sorted.bam

# Build an index (creates aln.sorted.bam.bai) for fast random access.
samtools index aln.sorted.bam

# Print just the first alignment record so we can inspect it.
samtools view aln.sorted.bam | head -n 1

SRR1770413.1	99	NC_000913.3	1	60	50M	=	312	361	AGCTTTTCATTCTGACTGCAACGGGCAATATGTCTCTGTGTGGATTAAAA	FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF	NM:i:0	MD:Z:50	MC:Z:50M	AS:i:50	XS:i:0

Reading A Record

Each SAM record is one line split into tab-separated columns. Column 1 (QNAME) is the read name, a unique label for that sequencing read. Column 3 (RNAME) names the reference sequence the read landed on, and column 4 (POS) is the 1-based position of the read's leftmost base on that sequence. Column 5 (MAPQ) is the mapping quality, a Phred-scaled score of how confident the aligner is that the position is correct; with BWA-MEM a value of 60 means highly confident and 0 means the read fits equally well in more than one place. Column 6 (CIGAR) summarizes the alignment shape, where 50M means 50 consecutive bases were aligned to the reference and letters like I, D, and S would mark insertions, deletions, and soft-clipped ends.