Zubair Khalid

Virologist/Molecular Biologist | Veterinarian | Bioinformatician

Conventional & Molecular Virology • Vaccine Development • Computational Biology

Dr. Zubair Khalid is a veterinarian and virologist specializing in conventional and molecular virology, vaccine development, and computational biology. Dedicated to advancing animal health through innovative research and multi-omics approaches.

Dr. Zubair Khalid - Veterinarian, Virologist, and Vaccine Development Researcher specializing in Computational Biology, Multi-omics, Animal Health, and Infectious Disease Research

Category: Guides

Bwa Genomics

BWA (Burrows,Wheeler Aligner) is a high‑performance software package that maps short sequencing reads to a large reference genome. It is the de facto tool for alignment in resequencing, variant calling, and many metagenomic pipelines. This guide is for bioinformaticians, molecular biologists, and clinical researchers who need a practical, source‑grounded understanding of BWA , how to choose the correct algorithm, execute a reliable workflow, spot common errors, and recognise where the tool stops being appropriate.

At a Glance

Feature Detail
Algorithm family Burrows‑Wheeler Transform, FM‑index
Main programs bwa aln (backtrack), bwa bwasw, bwa mem
Best suited for Short reads (70,1000 bp), Illumina, paired‑end or single‑end
Typical output SAM (Sequence Alignment/Map format)
Key strength Speed and memory efficiency for mammalian‑size genomes
Key limitation Not designed for long reads (PacBio, ONT) or bisulfite‑converted DNA

Core Concepts

BWA indexes the reference genome using the Burrows‑Wheeler Transform, then quickly identifies positions where each read most likely originated. The three algorithms serve different read lengths and error profiles. bwa aln (backtrack) works well for reads shorter than 100 bp with low divergence. bwa bwasw handles longer reads and gapped alignment. Most modern workflows use bwa mem, the default since BWA 0.7, because it is faster, supports paired‑end data natively, and tolerates higher divergence. The official documentation on the NCBI Bookshelf discusses the FM‑index principles behind these algorithms, and the EMBL‑EBI Training materials provide a beginner‑friendly introduction to sequence alignment concepts.

The essential input is a reference genome in FASTA format and one or two FASTQ files containing reads. The output is a SAM file (or compressed BAM/CRAM) that records mapping positions, quality scores, and pairing information. Downstream tools such as Samtools, Picard, and GATK then build on this alignment for variant calling, coverage analysis, or metagenomic profiling.

Decision Points

Before running BWA, evaluate your data to pick the correct algorithm and indexing strategy.

Read length and type

  • Reads <70 bp (e.g., ancient DNA or some ChIP‑seq): use bwa aln with -l to adjust the seed length. Alternatively, newer BWA versions still offer aln, but mem can handle very short reads with reduced sensitivity.
  • Reads 70,1000 bp, standard Illumina: use bwa mem. This covers whole‑genome resequencing, exomes, RNA‑seq (though STAR or HISAT2 are more appropriate for spliced alignments), and targeted amplicons.
  • Reads >1 kb or with high error rate (PacBio CCS, ONT): BWA is not optimal. Consider minimap2 or ngmlr.

Paired‑end, mate‑pair, or single‑end

  • bwa mem automatically recognises paired FASTQ files when you supply -p for interleaved or two separate files. It correctly flags proper pairs and inserts sizes.
  • bwa aln requires separate runs for each read and then bwa sample (deprecated) or bwa sampe to produce paired output. Use mem to simplify.

Reference genome size and repeats

  • For a human genome (3 Gb), default index parameters are adequate. For highly repetitive genomes or large plant genomes (e.g., wheat, 16 Gb), increase the -b band width or use bwa index with a larger block size. The Galaxy Training Network includes workflows that demonstrate how to adjust BWA parameters for non‑standard genomes.

Specialised applications

  • Methylation‑seq (bisulfite conversion): BWA does not handle C‑to‑T asymmetry. Use Bismark or BSMAP.
  • Metagenomic shotgun data: BWA‑MEM is commonly used, see the metagenomic dataset described in a recent study shotgun metagenomic dataset of leaf endophytic microbiome where BWA‑MEM aligned reads from Salvia officinalis endophytes. However, for very high diversity or viral samples, consider kraken2 or centrifuge.
  • Amplicon or targeted panels: BWA works, but consider bwa aln with strict mismatch thresholds to reduce off‑target alignments.

Practical Workflow

The following steps form a robust implementation sequence. All commands assume BWA 0.7.17 or later and Samtools 1.10+.

1. Prepare the reference index

bwa index reference.fa

This creates five files (.amb, .ann, .bwt, .pac, .sa). Index only once per reference version. Many public genomes (e.g., hg38, GRCh38) have pre‑built indexes available on the NCBI Sequence Read Archive or UCSC download sites. Using a pre‑indexed reference saves time.

2. Align reads with bwa mem

Basic paired‑end command:

bwa mem -t 8 reference.fa read1.fastq.gz read2.fastq.gz > aligned.sam
  • -t specifies threads. Use number of CPU cores minus one.
  • For single‑end: bwa mem reference.fa reads.fq > aligned.sam.
  • For interleaved paired FASTQ: bwa mem -p reference.fa interleaved.fq.

3. Convert SAM to sorted BAM

Sorting and indexing are essential for efficient downstream access.

samtools view -bS aligned.sam | samtools sort -o sorted.bam
samtools index sorted.bam

4. Mark optical/ PCR duplicates (optional but standard for DNA‑seq)

Use Picard MarkDuplicates or samtools markdup. Duplicate marking improves downstream variant calling accuracy.

5. Quality filtering (optional)

Remove reads with low mapping quality, unmapped pairs, or supplementary alignments.

samtools view -b -q 30 -F 0x04 sorted.bam > filtered.bam

6. Generate alignment statistics

samtools flagstat sorted.bam

This output reports total reads, mapped reads, properly paired reads, and secondary alignments. The Bioconductor ecosystem provides R packages (Rsamtools, GenomicAlignments) to compute coverage, insert‑size distributions, and enrichment.

Quality Checks and Validation

Three critical numbers must be inspected before any biological interpretation.

Mapping rate , At least 90% for high‑quality Illumina data against a complete reference. Lower rates suggest contamination, poor library quality, or a mismatch between organism and reference. In a recent study of Aspergillus fumigatus mitochondrial genomes, methods for large‑scale intraspecific SNP analyses used BWA and reported mapping rates above 95% after quality filtering.

Properly paired percentage , For paired‑end data, >80% of reads should align in the correct orientation and insert‑size range. A drop indicates possible chimera reads, adapter contamination, or inconsistent insert lengths.

Coverage uniformity , Use samtools depth or mosdepth to examine mean depth and standard deviation. Excessively high variance may point to amplification bias or collapsed repeats.

A real‑world example: in a metagenomic diagnostic study for febrile urinary tract infections, researchers used BWA‑MEM to align metagenomic reads and reported mapping statistics to confirm bacterial DNA enrichment metagenomic sequencing as a diagnostic tool for urine culture negative febrile urinary tract infection. Always include such metrics in publications to demonstrate data quality.

Common Mistakes

1. Incorrect algorithm for read length , Using bwa aln on 150‑bp reads yields slower performance and fewer properly paired alignments. Use bwa mem for reads 70 bp or longer.

2. Forgetting to index the reference , BWA will complain with an unhelpful error. Always check that all five index files exist.

3. Mixing FASTQ read‑pair order , Supplying read2.fastq before read1 will cause BWA to report many improper pairs. Use consistent naming (e.g., _R1 and _R2), and verify with a small subset.

4. Ignoring read groups , Variant callers (GATK, FreeBayes) require read group tags (@RG) in the SAM header. Add them during alignment:

bwa mem -R '@RG\tID:sample1\tSM:sample1\tPL:ILLUMINA' ...

Omitting this will force you to re‑run or use expensive BAM manipulation later.

5. Over‑trimming reads , Aggressive quality trimming before alignment can remove bases necessary for unique mapping. Light adapter trimming (e.g., with fastp or cutadapt) is fine, but avoid hard clipping all bases below Q20 unless the data are exceptionally poor.

Limits of Interpretation

BWA is not a universal aligner. Understand its boundaries to avoid over‑interpreting results.

Structural variant detection , BWA‑MEM can hint at deletions or duplications through discordant read pairs or split reads, but dedicated tools (DELLY, Manta, Lumpy) are required for confident calls. Relying only on BWA output for structural variant detection leads to high false‑positive and false‑negative rates.

Highly divergent sequences , BWA assumes low sequence divergence (typically <5%). If your organism diverges significantly from the reference (e.g., non‑model species, cross‑species alignment), consider a more tolerant aligner like minimap2 or LAST.

RNA‑seq , BWA does not handle splice junctions. Using it for RNA‑seq will map reads across exons incorrectly. Use STAR, HISAT2, or Salmon.

Methylation or bisulfite data , BWA will destroy C‑to‑T conversion information. Specialised aligners perform the necessary strand‑aware alignment.

Unmapped or multimapped reads , Reads that map to multiple positions (often repeats) are assigned a random location by BWA. The SAM field XS:i: reports the second‑best score, but handling multimapped reads requires careful filtering or rescuing tools (RScript using Bioconductor packages). Multimapped reads can inflate copy‑number estimates or introduce false SNP calls.

In a study of inherited variations in autism spectrum disorder families, researchers used BWA for alignment and then explicitly filtered out multimapped reads before their genotype‑phenotype correlation analysis in silico genotype‑phenotype correlation analysis. This is a necessary precaution.

Frequently Asked Questions

1. How do I choose between bwa aln and bwa mem for a short‑read project? For reads 70 bp or longer, always use bwa mem. It is faster, supports paired ends natively, and produces longer alignments. Only revert to bwa aln for read lengths below 70 bp or for legacy pipeline compatibility.

2. Why does my alignment output contain unmapped reads with high quality scores? Unmapped reads can arise from adapter contamination, low‑complexity sequences, or a reference that lacks the corresponding genomic region. Check adapter content with FastQC and consider trimming. Also verify the reference genome matches your sample species.

3. Can I use BWA for PacBio HiFi reads if I lower the seed length? BWA‑MEM works poorly with long reads ( >1 kb) even with reduced seed lengths. PacBio HiFi (CCS) reads have high accuracy but are long, minimap2 is the standard choice. BWA will produce many small alignments and a high number of supplementary records.

4. How do I measure mapping accuracy without a ground truth? You can inspect per‑base quality scores, the distribution of insert sizes (for paired ends), and the proportion of reads that map uniquely versus multi‑mapping. Simulating reads from your reference with tools like wgsim or art_illumina and then aligning them back provides an accuracy benchmark.

References and Further Reading

Related Articles