Genome Mapping Software: What to Use to Align Your Own Genome

If you have your own whole-genome sequencing data and you want to map it, the answer is short. Align Illumina short reads with bwa-mem2 against the GRCh38 no-alt analysis set, sort and compress to CRAM, then call variants with DeepVariant. DRAGMAP and Sentieon are alternatives if you want speed and can pay for it. For long reads, use minimap2 with the preset matched to your chemistry. Everything else described as “genome mapping software” is either a visualization tool such as IGV, UCSC Genome Browser or JBrowse, or a different laboratory technology entirely. Optical genome mapping on a Bionano Saphyr is one such technology, and it does not sequence bases at all. Knowing which of the three you need saves a great deal of wasted reading.

Three different things are called genome mapping

The oldest meaning is physical and genetic mapping: ordering markers and clones along a chromosome without knowing the sequence between them. That was the framework of the Human Genome Project’s first decade, and the reason the project was structured as mapping first and sequencing second ¹. It is history now for human work, because we have a finished reference. The informatics problems that project surfaced are still the same ones you will meet on your laptop, namely data volume, format standardization and distributed analysis ².

The second meaning, and the one that matters if you hold a FASTQ file, is read mapping: finding the position in a reference genome where each short DNA fragment most likely originated. This is the first computational step of every resequencing pipeline, and the quality of everything downstream depends on it ³.

The third meaning is optical genome mapping, in which very long DNA molecules are fluorescently labeled at a sequence motif and imaged in nanochannels. The output is a pattern of label distances, with no base calls. It resolves large structural rearrangements well and single nucleotide variants not at all. Unless you specifically ordered that assay, it is not your data.

The aligner we would use, and why

For paired-end Illumina reads, bwa-mem2 is what we run. It is a reimplementation of BWA-MEM with the same alignment output and roughly two to three times the throughput, at the cost of a much larger index (about 80 GB of RAM for a human reference build, which is the main reason people fall back to original bwa mem on small machines). A minimal invocation:

bwa-mem2 mem -t 32 -K 100000000 -Y \
  -R '@RG\tID:L001\tSM:sample01\tPL:ILLUMINA\tLB:lib1' \
  GRCh38_no_alt.fa R1.fastq.gz R2.fastq.gz \
| samtools sort -@ 8 -m 2G -o sample01.bam -
samtools index -@ 8 sample01.bam

The flags matter. -K 100000000 fixes the number of bases processed per batch so results are deterministic regardless of thread count, which you want if you ever intend to reproduce a run. -Y uses soft clipping for supplementary alignments, which preserves the full read sequence and helps structural variant callers that parse clipped segments. The -R read group is not optional in practice: GATK, DeepVariant, and duplicate markers all expect it, and adding it later means rewriting the file.

The alternative worth knowing is Sentieon’s BWA implementation or Illumina’s DRAGEN, which push a 30x genome from FASTQ to VCF in well under an hour on the right hardware, versus roughly four to eight hours of wall time for bwa-mem2 on 32 cores. Both are commercial. DRAGMAP is the open-source reimplementation of the DRAGEN mapper and is a middle path. It is what the current GATK best-practice reference pipelines use. We stay on bwa-mem2 for personal work because its failure modes are extremely well characterized across a decade of population-scale use, and because matching the aligner used to build the large public variant catalogs makes your call set easier to compare against them ⁴.

For PacBio HiFi or Oxford Nanopore data, use minimap2 -ax map-hifi or -ax lr:hq respectively, then sort the same way. Long reads change the calculus substantially: mappability in segmental duplications improves, and structural variants become tractable with sniffles2 or pbsv.

Reference choice is the decision that compounds

Use the GRCh38 analysis set with no ALT contigs, plus the decoy sequences and HLA contigs, unless you have a specific reason to do otherwise. ALT contigs represent alternate haplotypes for polymorphic regions, and an aligner that is not ALT-aware will distribute reads between the primary and alternate copies, dropping mapping quality to zero in exactly the regions you were curious about. If you do include ALTs, you must run bwa-postalt.js from bwakit to lift alignments back, and most downstream callers still will not handle them correctly.

T2T-CHM13v2.0, the first truly complete human assembly, is better than GRCh38 in acrocentric arms, segmental duplications, and centromeric regions. It also removes several thousand collapsed or erroneous bases that generate persistent false positive calls. The tradeoff is annotation ecosystem: clinical variant databases, gnomAD frequencies, and transcript coordinate systems are still primarily GRCh38. Our recommendation is to make GRCh38 the primary alignment and run T2T as a second pass if a region of interest looks suspicious. Comparing the two is often more informative than either alone.

One detail people miss: for a sample with two X chromosomes, leave chrY in the reference but rely on the analysis set’s hard-masked pseudoautosomal regions on Y, so PAR reads go to X rather than being split. Getting this wrong produces phantom heterozygous calls on Y.

Files, sizes, and the storage decision

A 30x human genome in paired 2x150 FASTQ is roughly 90 to 110 GB gzipped. The sorted BAM is about 60 to 80 GB. Converted to CRAM 3.1 against the same reference, it is roughly 15 to 25 GB, and lossless:

samtools view -@ 8 -T GRCh38_no_alt.fa -C -o sample01.cram sample01.bam
samtools index sample01.cram

Keep the CRAM and the reference FASTA together, permanently. CRAM is reference-differential, so a CRAM without the exact reference (matched by MD5 of each sequence) is close to unreadable. Keep the FASTQ only until you trust the CRAM, then decide: FASTQ lets you re-align to any future reference, CRAM lets you re-call variants but locks in the alignment. We keep FASTQ in cold object storage and CRAM in warm storage.

Where mapping quietly fails

Between five and ten percent of reads in a typical human genome land with mapping quality zero, meaning the aligner found two or more equally good positions. Those reads are concentrated in the regions where interesting biology hides, and they are not noise. The regions include SMN1/SMN2, CYP2D6 and its pseudogene, PMS2/PMS2CL, GBA/GBAP1, and the HLA locus. Short-read mapping is structurally unable to resolve some of these, and no amount of extra coverage fixes it. Assays targeting those genes specifically, or long-read data, are the remedy. This limitation of personal short-read genomes has been evident since the earliest individual genome papers ⁵, and careful work has shown that coverage depth and platform choice both change which fraction of the genome is confidently callable ⁶.

Before trusting any call set, run the checks. samtools flagstat for duplicate and properly-paired rates (a duplicate rate above roughly 15 to 20 percent suggests a low-complexity library). mosdepth for coverage uniformity and the fraction of the callable genome above 10x. verifybamid2 for cross-sample contamination, where anything above about 1 to 2 percent will start producing spurious heterozygous calls. somalier to confirm the sample is who you think it is.

What the compute costs

The sequencing is the expense. The mapping is nearly free. On a cloud spot instance with 32 to 64 vCPUs and enough RAM for the bwa-mem2 index, alignment plus sorting plus duplicate marking for a 30x genome runs a few dollars of compute, and DeepVariant adds a few more (much less on a GPU instance). Storage of the CRAM in standard object storage is well under a dollar a month. A used workstation with 128 GB of RAM and fast NVMe will do the same job overnight for the cost of electricity. The economics of sequencing itself have collapsed by orders of magnitude since the first drafts, and the projections made twenty years ago about a thousand-dollar genome have been met and passed ⁷. Analysis cost followed the same curve, which is why the bottleneck is now interpretation rather than either sequencing or alignment ⁸.

A mapped genome and a called VCF tell you what bases you carry. They do not tell you what those bases mean for your health, and the gap between the two is wide. If a variant in your data appears in a medically actionable gene, the correct next step is a genetics clinician and confirmatory testing in an accredited laboratory. A deeper dive into your own BAM file is not the next step.

Questions people also ask

How much does genome mapping cost?

If you mean the compute to align reads you already have, roughly five to twenty dollars of cloud time per 30x genome, or a night on your own machine. If you mean the sequencing that produces the reads, a 30x whole genome from a commercial provider currently runs from a few hundred dollars for data-only service to several thousand for a clinical-grade report. Optical genome mapping is a separate laboratory assay, priced per sample through the vendor, and is not something you can run on existing sequence data.

Is the UCSC Genome Browser a genome mapping tool?

It is a visualization and annotation browser. You use it to look at a region of the reference alongside gene models, conservation tracks, and variant databases. You can also load your own BAM or VCF as a custom track. It does not align reads. For inspecting your own alignments read by read, IGV running locally against your CRAM is faster and more informative.

Do I need to align reads myself if my provider gave me a VCF?

Not necessarily. The VCF is still a lossy summary of someone else’s parameter choices. Having the CRAM lets you check whether a surprising call is supported by well-mapped reads or by three reads at the edge of a repeat, which is the single most useful thing you can do with your own genome. Population-scale projects made the same point institutionally: the value of a variant map depends on being able to return to the underlying alignments ⁴.

How much coverage do I need?

Thirty-fold average is the working standard for germline short-read variant calling, and it is what the large reference cohorts used ³. Higher depth improves sensitivity in low-coverage regions and for detecting mosaicism, but it does not rescue regions where mapping quality is zero for structural reasons.

Footnotes

1. Lloyd Smith, Leroy Hood. Mapping and Sequencing the Human Genome: How to Proceed. Nature Biotechnology, 1987. https://doi.org/10.1038/nbt0987-933 ↩

2. Karen A. Frenkel. The human genome project and informatics. Communications of the ACM, 1991. https://doi.org/10.1145/125490.125492 ↩

3. Rachel L Goldfeder, Dennis P Wall, Muin J Khoury, et al. Human Genome Sequencing at the Population Scale: A Primer on High-Throughput DNA Sequencing and Analysis. American Journal of Epidemiology, 2017. https://doi.org/10.1093/aje/kww224 ↩ ↩²

4. The 1000 Genomes Project Consortium. A map of human genome variation from population-scale sequencing. Nature, 2010. https://doi.org/10.1038/nature09534 ↩ ↩²

5. Michael Snyder, Jiang Du, Mark Gerstein. Personal genome sequencing: current approaches and challenges. Genes & Development, 2010. https://doi.org/10.1101/gad.1864110 ↩

6. Subramanian S. Ajay, Stephen C.J. Parker, Hatice Ozel Abaan, et al. Accurate and comprehensive sequencing of personal genomes. Genome Research, 2011. https://doi.org/10.1101/gr.123638.111 ↩

7. Simon T Bennett, Colin Barnes, Anthony Cox, et al. Toward the $1000 Human Genome. Pharmacogenomics, 2005. https://doi.org/10.1517/14622416.6.4.373 ↩

8. Leroy Hood, Lee Rowen. The human genome project: big science transforms biology and medicine. Genome Medicine, 2013. https://doi.org/10.1186/gm483 ↩