Aligning Your Whole-Genome FASTQs with BWA-MEM

By the end of this guide you will have taken a pair of gzipped FASTQ files from a 30x human whole-genome run and produced a coordinate-sorted, duplicate-marked CRAM aligned to GRCh38, along with the index and the QC numbers you need to decide whether the alignment is any good. You need roughly 64 GB of RAM, 16 or more cores, and about 300 GB of free disk for intermediates. BWA’s index for the human genome occupies about 5 GB resident, but sorting and duplicate marking are the real memory consumers. The final CRAM will be around 15 to 20 GB. Software: bwa 0.7.17 or newer, samtools 1.19 or newer, and either samtools markdup or GATK’s MarkDuplicates. Everything below assumes a Linux box with a real filesystem, not a network mount, because alignment is I/O-bound at the write end more often than people expect.

1. Understand what BWA-MEM is doing before you run it

BWA stands for Burrows-Wheeler Aligner, named after the Burrows-Wheeler transform, a reversible reordering of the reference sequence that makes it possible to build an FM-index supporting fast substring queries in space proportional to the genome itself ¹. The package contains three algorithms. BWA-backtrack was designed for Illumina reads up to about 100 bp and is now obsolete ¹. BWA-SW extended the approach to longer, more divergent reads using a Smith-Waterman-like extension off seeds ². BWA-MEM is the one you want. It finds maximal exact matches (MEMs) between the read and the reference using the FM-index. It then chains the seeds that are colinear and extends the best chain with a banded affine-gap alignment.

Two properties of BWA-MEM matter for interpreting your output. First, it performs local alignment. A read that only partially matches the reference will be soft-clipped rather than forced into a full-length alignment. That covers reads spanning a structural breakpoint or carrying adapter at the 3’ end. Second, a single read can produce one primary alignment plus supplementary alignments for other segments. Those are flagged with SAM flag 0x800, which is how split reads that support structural variants are represented. Any downstream tool that silently drops supplementary records will lose that evidence.

On the question of whether BWA-MEM is still the best choice: for short paired-end Illumina reads against a human reference, it is the aligner whose behavior downstream variant callers have been tuned against, and comparative benchmarks of next-generation sequencing aligners have consistently placed BWA among the most accurate on realistic error and indel models ³. The main alternatives worth knowing are Bowtie2 and Dragen-style hardware or the open-source bwa-mem2. Bowtie2 is comparable in accuracy but handles indels less gracefully. bwa-mem2 reimplements the same algorithm with SIMD-optimized banded extension and produces near-identical output roughly two to three times faster. There is also BWA-MEME, which replaces the FM-index seeding step with a learned index and reports substantial speedups on the seeding phase while preserving BWA-MEM’s output ⁴. If you have the RAM (bwa-mem2’s index needs about 28 GB for GRCh38, BWA-MEME more still), use bwa-mem2. The command line is otherwise identical to what follows.

2. Get the right reference and index it

The reference choice affects your variant calls more than the aligner flags do. Use the GRCh38 analysis set that includes the decoy contigs and the ALT contigs but has the pseudoautosomal regions on chrY hard-masked. The canonical file is GCA_000001405.15_GRCh38_no_alt_analysis_set.fna if you want to skip ALT-aware handling, or the full analysis set with .alt file if you want it. We recommend the no-ALT analysis set with decoys for a first pass, because ALT-aware postprocessing (bwa-postalt.js) adds a step that most downstream germline pipelines do not expect, while the decoy contigs soak up reads from unplaced sequence that would otherwise pile up as false variants in the real chromosomes.

REF=GCA_000001405.15_GRCh38_no_alt_analysis_set.fna
bwa index -a bwtsw "$REF"
samtools faidx "$REF"
samtools dict "$REF" -o "${REF%.fna}.dict"

bwa index on a human genome takes 50 to 70 minutes on a single core and cannot be parallelized. The -a bwtsw flag selects the BWT construction algorithm suited to large genomes. The alternative is is, which fails above 2 GB. You will end up with five files totaling about 5.5 GB. Their extensions are .amb, .ann, .bwt, .pac and .sa. Index once, keep it on fast local storage, and reuse it. If you run this repeatedly across samples, the index load time (about 30 seconds from warm page cache, several minutes cold) is worth avoiding by batching samples in one process rather than one container per sample.

3. Check your FASTQs before you align

Alignment will happily consume garbage and give you a plausible-looking BAM. Spend ten minutes first. Run fastqc or, faster, seqkit stats:

seqkit stats -a sample_R1.fastq.gz sample_R2.fastq.gz

Confirm that read counts match between R1 and R2, that read length is what you expected (150 bp for most current WGS), and that the total base count divided by 3.1e9 gives you the coverage you paid for. A 30x human genome at 2x150 is roughly 310 million read pairs. Check the quality encoding is Phred+33, which everything modern uses, and glance at the adapter content. BWA-MEM soft-clips adapter reasonably well, so pre-trimming is usually unnecessary and occasionally harmful, because trimming tools sometimes truncate reads in a way that biases indel representation. Trim only if adapter content exceeds a few percent, and if you do, use fastp --detect_adapter_for_pe rather than a fixed adapter string.

Also extract the read group information now. You need it at alignment time because adding it afterward means rewriting the whole file.

zcat sample_R1.fastq.gz | head -1
# @A00123:456:HXXXXDSXX:1:1101:1181:1000 1:N:0:ATCACGTT

The instrument, run number, flowcell, and lane give you a platform unit. Construct the read group string:

RG='@RG\tID:HXXXXDSXX.1\tSM:sample01\tLB:lib1\tPU:HXXXXDSXX.1.ATCACGTT\tPL:ILLUMINA'

SM is the sample name and is what most downstream tools key on. ID must be unique across every read group in the final file. If your sample was sequenced across four lanes, align each lane separately with its own ID and PU, then merge. This preserves the ability to model lane-specific error behavior in base quality recalibration.

4. Run the alignment

The core command pipes BWA’s SAM output straight into samtools so you never write an uncompressed SAM to disk.

bwa mem -t 16 -K 100000000 -Y \
    -R "$RG" \
    "$REF" \
    sample_R1.fastq.gz sample_R2.fastq.gz \
  | samtools fixmate -m -u -@ 4 - - \
  | samtools sort -u -@ 8 -m 3G -T /scratch/sort_tmp - \
  | samtools markdup -@ 8 --reference "$REF" --output-fmt CRAM - sample01.cram

samtools index sample01.cram

The flags that matter:

-t 16 sets alignment threads. BWA-MEM scales close to linearly up to about 16 threads and then flattens as the single-threaded output writer becomes the bottleneck. Work on parallelizing BWA’s single-end path has shown the same pattern, where thread scaling is limited by stages outside the alignment kernel itself ⁵. If you have 64 cores, you will get more throughput running four samples at -t 16 than one at -t 64.

-K 100000000 fixes the number of bases loaded per processing batch. Without it, BWA estimates the insert size distribution from whatever chunk happens to be in memory, which makes output depend on thread count and batch boundaries. Setting -K makes the run deterministic, which you want if you ever need to reproduce a result.

-Y uses soft clipping rather than hard clipping for supplementary alignments, so the full read sequence is retained in every record. Structural variant callers such as Manta and GRIDSS expect this.

samtools fixmate -m adds the mate score tag that markdup needs to pick which duplicate to keep. The -u flags on fixmate and sort request uncompressed BAM between pipeline stages, which saves considerable CPU that would otherwise go to compressing and decompressing intermediate data.

-m 3G on sort is per-thread memory, so eight threads will use about 24 GB. Set -T to a directory on a fast local disk with at least 200 GB free. Sort writes temporary chunks there, and the default of the output directory often is not what you want.

Writing CRAM instead of BAM costs a little CPU and saves roughly 40 percent of the disk, and CRAM is reference-based, so keep the exact reference FASTA you aligned against forever. Losing it means losing the ability to decode the file.

Expect 6 to 10 hours of wall clock for a 30x genome at -t 16 on a modern server, most of it in BWA itself. If you are running many samples and the cost matters, distributed implementations exist: SparkBWA wraps BWA-MEM in Spark and reports near-linear scaling across a cluster without changing the alignment results ⁶. For a single personal genome, the added complexity is not worth it.

5. Mark duplicates correctly

The command above uses samtools markdup, which is fast and adequate for PCR-free libraries. If your library was PCR-amplified, or you care about optical duplicates specifically, use GATK:

gatk MarkDuplicates \
    -I sample01.bam \
    -O sample01.md.bam \
    -M sample01.dup_metrics.txt \
    --OPTICAL_DUPLICATE_PIXEL_DISTANCE 2500 \
    --CREATE_INDEX true

The pixel distance of 2500 is the right value for patterned flowcells (NovaSeq, HiSeq X). The default of 100 was tuned for unpatterned ones and will undercount optical duplicates badly on a NovaSeq. Read the metrics file: for a PCR-free 30x library you should see a duplicate fraction under 5 percent, often under 2 percent. Ten percent or higher means the library was over-amplified relative to its complexity, and your effective coverage is lower than the raw base count suggests.

6. Verify the alignment

Never move to variant calling without looking at the numbers. Three commands give you most of what you need.

samtools flagstat -@ 8 sample01.cram
samtools stats -@ 8 sample01.cram > sample01.stats
mosdepth --by 1000 --fast-mode -t 4 sample01 sample01.cram

From flagstat, expect above 99 percent of reads mapped and above 95 percent properly paired for a good human WGS library. Mapped fractions in the 80s usually mean contamination or the wrong reference. From samtools stats, look at the insert size distribution: a PCR-free library should show a roughly unimodal peak between 300 and 500 bp with a standard deviation under about 100. A long right tail or a second mode suggests chimeric molecules from ligation artifacts. The error rate line should be around 0.2 to 0.5 percent for Illumina.

From mosdepth, check mean coverage. More importantly, check the fraction of the genome at 10x or more. For a nominal 30x genome, you want above 90 percent of the callable genome at 10x. Coverage uniformity, not mean depth, determines how much of your genome you can call confidently.

One last check worth doing: mapping quality distribution. BWA-MEM assigns MAPQ 0 to reads that map equally well to multiple places, and about 5 to 8 percent of a human genome falls in such regions. If your MAPQ 0 fraction is much higher, either the reference lacks decoys (so junk reads scatter across real chromosomes) or you are aligning to the wrong build.

Common problems

The alignment runs but almost nothing maps. Nearly always the wrong reference, a corrupted index, or FASTQs that are not what you think they are. Take a thousand reads and BLAST a few by hand to see what organism and what region they come from. Local alignment search tools remain the quickest way to identify an unknown sequence when you need a sanity check rather than a full alignment ⁷.

Reads map but the properly-paired rate is low. Check that R1 and R2 are in the same order and the same length. If a FASTQ was truncated during transfer, BWA will still run and will silently pair the wrong reads once the files desynchronize. seqkit stats catching a read-count mismatch is the cheap way to catch this. Also confirm you did not pass the same file twice.

samtools sort dies with “No space left on device”. The -T directory filled up. Sorting a 30x genome needs about 100 to 150 GB of temporary space even when the final output is 15 GB, because the intermediates are uncompressed or lightly compressed.

markdup reports zero duplicates. You skipped samtools fixmate -m, or you ran markdup on a name-sorted rather than coordinate-sorted file. The order of the pipeline stages is not optional.

Output differs between runs on the same input. You omitted -K. BWA’s insert size estimation is batch-dependent, so different thread counts produce slightly different MAPQ values and occasionally different alignments for ambiguous reads. Fix -K and the run becomes reproducible.

You are aligning something other than short Illumina reads. BWA-MEM handles reads up to around a few kilobases and was extended from the original short-read design specifically to cope with longer, more divergent sequences ². For PacBio HiFi or Oxford Nanopore, use minimap2 instead. BWA-MEM will be slow and will misplace reads with high indel rates. And if what you have is protein or a set of homologous coding sequences to compare across species, BWA is the wrong tool category entirely. That is multiple sequence alignment, where tools work on evolutionary similarity rather than mapping to a reference. Aligning protein-coding DNA is often better done by translating to amino acids, aligning there, and back-translating ⁸. Benchmarking work on noncoding alignment shows how differently these tools behave from read mappers and from each other ⁹.

Two final notes on interpretation. A BAM or CRAM is a hypothesis about where each read came from, not a fact, and MAPQ is the aligner’s own estimate of how likely it is to be wrong. Anything you conclude about a specific variant in a low-complexity or segmentally duplicated region deserves manual inspection in IGV before you believe it. And nothing in an alignment pipeline is a clinical result. If a variant you find looks medically meaningful, confirm it with an orthogonal assay in an accredited laboratory and discuss it with a clinician or genetic counselor before acting on it in any way.

Footnotes

1. Heng Li, Richard Durbin. Fast and accurate short read alignment with Burrows–Wheeler transform. Bioinformatics, 2009. https://doi.org/10.1093/bioinformatics/btp324 ↩ ↩²

2. Heng Li, Richard Durbin. Fast and accurate long-read alignment with Burrows–Wheeler transform. Bioinformatics, 2010. https://doi.org/10.1093/bioinformatics/btp698 ↩ ↩²

3. Matthew Ruffalo, Thomas LaFramboise, Mehmet Koyutürk. Comparative analysis of algorithms for next-generation sequencing read alignment. Bioinformatics, 2011. https://doi.org/10.1093/bioinformatics/btr477 ↩

4. Youngmok Jung, Dongsu Han. BWA-MEME: BWA-MEM emulated with a machine learning approach. Bioinformatics, 2022. https://doi.org/10.1093/bioinformatics/btac137 ↩

5. Heeseung Jo, Gunhwan Koh. Faster single-end alignment generation utilizing multi-thread for BWA. Bio-Medical Materials and Engineering, 2015. https://doi.org/10.3233/bme-151480 ↩

6. José M. Abuín, Juan C. Pichel, Tomás F. Pena, et al. SparkBWA: Speeding Up the Alignment of High-Throughput DNA Sequencing Data. PLOS ONE, 2016. https://doi.org/10.1371/journal.pone.0155461 ↩

7. S. Donkor Eric, T. K. D. Dayie Nicholas, K Adiku Theophilus. Bioinformatics with basic local alignment search tool (BLAST) and fast alignment (FASTA). Journal of Bioinformatics and Sequence Analysis, 2014. https://doi.org/10.5897/ijbc2013.0086 ↩

8. Olaf RP Bininda-Emonds. transAlign: using amino acids to facilitate the multiple alignment of protein-coding DNA sequences. BMC Bioinformatics, 2005. https://doi.org/10.1186/1471-2105-6-156 ↩

9. Daniel A Pollard, Casey M Bergman, Jens Stoye, et al. Benchmarking tools for the alignment of functional noncoding DNA. BMC Bioinformatics, 2004. https://doi.org/10.1186/1471-2105-5-6 ↩