4: Reference Genome Architecture, Indexing & Read Alignment

Chapter 4: Reference Genome Architecture, Indexing & Read Alignment

Read alignment is the computational core of resequencing genomics. High-throughput sequencers generate hundreds of millions of short read fragments per sample. The aligner must rapidly map each read to its true chromosomal origin within a 3.1-billion base reference genome in the presence of sequencing errors and biological polymorphisms.

1. The Burrows-Wheeler Transform & FM-Index

Linear scanning through a 3.1 Gb genome for each 150 bp read would take centuries. Modern aligners index the reference genome using the Burrows-Wheeler Transform (BWT) coupled with an auxiliary Ferragina-Manzini (FM) index:

• BWT Compression: Reversibly permutes characters of the reference string based on lexicographical sorting of all cyclic suffixes. Characters following identical prefixes are clustered together, creating long runs of repeated nucleotides that compress efficiently into memory.
• Backward Search: Queries exact seed matches character-by-character from right-to-left in time proportional strictly to the read length \(O(m)\), completely independent of the genome size \(G\).

2. Production Read Alignment with BWA-MEM2

bwa-mem2 utilizes AVX-512 and AVX2 vector instructions to accelerate alignment 2–3× over classical BWA-MEM:

# 1. Index the reference FASTA (only run once)
bwa-mem2 index -p GRCh38_noalt_decoy GRCh38_noalt_decoy.fasta

# 2. Align paired-end reads and attach Read Group (RG) headers
# Note: GATK and Picard REQUIRE @RG metadata (ID, SM, PL, LB)
bwa-mem2 mem \
    -t 16 \
    -R '@RG\tID:SAMPLE1_L1\tSM:SAMPLE1\tPL:ILLUMINA\tLB:LIB1\tPU:FLOWCELL_L1' \
    GRCh38_noalt_decoy \
    sample_trimmed_R1.fastq.gz \
    sample_trimmed_R2.fastq.gz \
    | samtools sort -@ 8 -m 4G -o sample1_aligned.sorted.bam -

# 3. Create BAM index (.bai)
samtools index -@ 8 sample1_aligned.sorted.bam

3. Long-Read Alignment with minimap2

For PacBio HiFi or Oxford Nanopore reads, minimap2 uses approximate seed matching via minimizers:

# Align PacBio HiFi reads
minimap2 -ax map-hifi -t 16 GRCh38.fasta pacbio_hifi.fastq.gz | samtools sort -o pacbio.bam -

# Align Oxford Nanopore genomic reads
minimap2 -ax map-ont -t 16 GRCh38.fasta ont_reads.fastq.gz | samtools sort -o ont.bam -