Alignment: putting reads back on the genome
The sequencer produces FASTQ: sequences and quality values with no position attached. Alignment puts each read back where it most likely came from on the reference genome, producing a BAM file. Two clean-up steps follow: marking PCR duplicates and recalibrating base qualities. The variant caller reads this cleaned BAM.
Alignment: finding each read’s place
The GATK Best Practices align reads to the reference genome with BWA-MEM. New projects mostly use GRCh38, older ones GRCh37; the TCGA-LAML MAF is in GRCh37 coordinates. Paired reads are aligned together: facing each other at a sensible distance, they form a properly paired read pair. The result is written as SAM/BAM (or the more compact CRAM), sorted by coordinate and indexed, so that the reads of any region can be fetched quickly.
Every read carries a mapping quality (MAPQ), the confidence that its position is right. Repeats, pseudogenes and gene families give a short read several equally good places, and BWA assigns such reads MAPQ 0. Callers ignore reads with very low MAPQ, so these regions are effectively invisible in WES.
Marking duplicates and recalibrating base qualities
Duplicate marking. MarkDuplicates finds read pairs whose two ends start at the same positions in the same orientation, treats them as copies of one original molecule, keeps the best pair and flags the rest (without deleting them). Optical duplicates are recognised here too. This is where the duplicate rate is computed. Libraries with UMIs can tell apart different molecules that happen to share a start; standard WES mostly has none.
Base quality recalibration. The base qualities a sequencer reports are systematically off. BQSR counts the actual mismatches away from known variant sites, models them by read group, cycle and neighbouring bases, and rewrites each base quality. Mutect2 weighs each read’s evidence by base quality, so recalibrated values are more trustworthy.
Older papers also mention local indel realignment. GATK4 no longer needs it, because Mutect2 reassembles each active region locally.
What the BAM records, and which QC metric it feeds
| In the BAM | Produced by | QC metric |
|---|---|---|
| Position and CIGAR | BWA-MEM alignment | Mapping rate, soft-clip fraction |
| MAPQ | Uniqueness of the position | Fraction of low-MAPQ reads |
| TLEN (template length) | Distance between the paired reads | Insert-size distribution |
| Duplicate flag (1024) | MarkDuplicates | Duplicate rate |
| Base qualities | Sequencer, recalibrated by BQSR | Fraction of bases at ≥Q30 |
| Read group (RG) | Sample, library and lane | Separates batches for duplicate marking and BQSR |
A high soft-clip fraction usually means inserts so short that the reads ran into the adapter; it can also mark structural-variant breakpoints.
Common pitfalls
- Mixing reference builds. GRCh37 and GRCh38 have different coordinates, and chromosome names may or may not carry a “chr” prefix. BAM, BED, panel of normals and annotation must all use the same build.
- Processing tumour and normal differently. Both need the same reference, aligner and settings; otherwise the pipeline difference itself creates false “somatic” signal.
- Assuming low-MAPQ regions have no mutations. They are just unseen. Genes with close pseudogenes, such as PMS2, need dedicated methods.