Sequence Alignment

intermediate

Alignment Pipeline

BWA-MEM2 Alignment Steps

🔨

Build BWA Index (once)

bwa-mem2 index hg38.fa - takes ~1 hour; stored permanently

🗺

Align Paired Reads

bwa-mem2 mem -t 16 -R "@RG ID:sample..." hg38.fa R1.fastq.gz R2.fastq.gz

📋

Sort BAM

samtools sort -@ 8 -o sample.sorted.bam - coordinate sort required for GATK

📍

Index BAM

samtools index sample.sorted.bam - creates .bai file for random access

📊

Alignment QC

samtools flagstat; mosdepth coverage; VerifyBamID2 contamination check

Reference Genomes

GRCh38/hg38 (2013) - current standard; better repeat resolution, alternative loci, recommend for new projects
GRCh37/hg19 (2009) - still used for legacy datasets and ClinVar/HGMD compatibility
T2T-CHM13 (2022) - first truly complete human genome; 200 Mb of new sequence resolved
Never mix reference versions - variants called on hg19 cannot be compared with hg38 without liftover
GATK resource bundle: curated reference files for hg38 (FASTA, dbSNP, Mills, 1000G) - download from Broad
Decoy sequences: add hs38d1 decoy to reduce false alignments to highly repetitive regions

BWA-MEM2 (Recommended)

BWA-MEM2 is 2× faster than BWA-MEM with identical results. The standard aligner for Illumina short reads in clinical pipelines.

code

# Build index (once, ~1 hour for hg38)
bwa-mem2 index hg38.fa

# Align paired-end reads - pipe directly to samtools for efficiency
bwa-mem2 mem \
  -t 16 \
  -R "@RG\tID:sample1\tSM:patient001\tPL:ILLUMINA\tLB:lib1\tPU:run1" \
  hg38.fa \
  R1_clean.fastq.gz R2_clean.fastq.gz \
  | samtools sort -@ 8 -o sample.sorted.bam

# Index BAM
samtools index sample.sorted.bam

# Quick alignment stats
samtools flagstat sample.sorted.bam
samtools idxstats sample.sorted.bam | head

Read Groups (Critical)

Read groups (@RG) are mandatory for GATK. They track which reads came from which flow cell, lane, and library - essential for detecting batch effects and running BQSR correctly.

ID: unique identifier for this read group (e.g., flowcell.lane)
SM: sample name - must match across RGs from same patient
PL: platform (ILLUMINA, PACBIO, NANOPORE)
LB: library preparation ID - identifies PCR duplicates across lanes
PU: platform unit (flowcell_barcode.lane.sample_barcode) - most specific
CN: sequencing centre name

CIGAR Strings Explained

M - alignment Match (can be match or mismatch)
I - Insertion relative to reference
D - Deletion relative to reference
N - Skipped region (intron in RNA-seq)
S - Soft clipping (bases present in read, not aligned)
H - Hard clipping (bases not in read sequence)
Example: 50M2D30M1I20M = 50 aligned, 2 deleted from ref, 30 aligned, 1 inserted, 20 aligned

Read Trimming Post-Alignment Processing