Genome comparison tools: which one to use for which comparison

Oak
genomics bioinformatics variant-calling pangenome microbiome
Two glowing DNA-like strands held on parallel rails in a dark lab instrument, linked by light where they match and flaring where they differ.

There is no single best genome comparison tool because “compare two genomes” means at least four different computations. If you have variant calls (VCF/gVCF) and want to know how two samples differ, use bcftools isec after normalization, and hap.py/rtg vcfeval when one side is a truth set. If you have two assemblies (FASTA), use minimap2 -cx asm5 plus paftools.js call, or MUMmer’s nucmer/dnadiff when you want a one-page identity report. If you are comparing many bacterial genomes from a metagenome, start with k-mer distances (skani, mash, sourmash), then move to a pangenome tool (panaroo) or a clustering tool (PopPUNK) once you know which genomes are close enough to align. If you want to compare across species to find conserved non-coding sequence, that is an alignment-and-visualization problem, and VISTA-style pipelines still do it well 1.

Pick the layer first. Most wasted effort in comparative genomics comes from running an assembly aligner on data that should have stayed in VCF space, or from comparing VCFs that were never normalized against the same reference.

Comparing your variants to another person, or to a truth set

Two VCFs from the same pipeline are rarely directly comparable, and two VCFs from different pipelines never are. An indel can be left-aligned or not, a MNP can be one record or three, multi-allelic sites can be split or merged. Normalize both sides first:

bcftools norm -m -any -f GRCh38.fa -Ou in.vcf.gz \
  | bcftools norm -d exact -Oz -o norm.vcf.gz
bcftools index norm.vcf.gz

Then set differences are cheap:

bcftools isec -p isec_out -n=1 -c none a.norm.vcf.gz b.norm.vcf.gz

-c none forces exact allele matching. -c some is more forgiving and will hide real representation differences, so use it only when you know why.

For anything resembling accuracy measurement, do not use set intersection. Use hap.py with the vcfeval engine, which does haplotype-aware comparison and will correctly score a 4 bp deletion written two different ways as a match:

hap.py truth.vcf.gz query.vcf.gz \
  -f HG002_GRCh38_benchmark.bed \
  -r GRCh38.fa --engine=vcfeval -o eval

The benchmark material worth knowing is Genome in a Bottle, built by sequencing seven genomes across many platforms and reconciling the calls into a high-confidence variant set plus a BED of regions where the truth is trusted 2. That BED is the important half. It deliberately excludes segmental duplications, most of the MHC, and other hard regions, so a 99.5% SNV F-score inside it says nothing about the genes that live outside it. If you are benchmarking your own pipeline, run the same reference material through it and compare in truth-set space rather than trusting reported Q-scores.

Two practical constraints. Short-read sequencing is good at SNVs and small indels and systematically weak at structural variants, repeat expansions, and anything in low-complexity or highly similar sequence 3. And reference build mismatch (hg19 vs GRCh38 vs T2T-CHM13) invalidates a comparison silently: coordinates shift, chromosome names change (chr1 vs 1), and CrossMap/liftOver drops a non-trivial fraction of variants at build boundaries. Check bcftools view -h for the ##reference line and ##contig lengths before doing anything else.

Comparing assemblies, not variant calls

If you have long-read data and an actual assembly, alignment beats variant calling for large-scale comparison. Our default:

minimap2 -cx asm5 --cs -t 16 ref.fa query.fa > aln.paf
paftools.js call -L 5000 -f ref.fa aln.paf > var.vcf

asm5 assumes under ~5% divergence, which is right for two human haplotypes. Use asm10 or asm20 for cross-species. --cs emits the difference string that paftools.js needs. -L 5000 restricts calling to alignment blocks of at least 5 kb, which suppresses most spurious calls from repeat-induced misalignment.

MUMmer is still the fastest way to get a summary:

nucmer --maxmatch -c 100 -p out ref.fa query.fa
dnadiff -d out.delta -p diff

diff.report gives average identity, aligned fraction of each sequence, and counts of insertions, relocations, and inversions in one screen. For inversions and translocations specifically, run SyRI on the alignment and plot with plotsr. For a quick visual check, D-GENIES will render a dotplot from a PAF in the browser, which is the fastest way to notice that your assembly has a chimeric contig.

Human assembly comparison increasingly means graph comparison rather than pairwise alignment. A single linear reference misrepresents the large fraction of human sequence that is polymorphic in presence or absence, and pangenome graphs built from many assemblies capture variation that pairwise alignment to GRCh38 hides 4. minigraph-cactus and pggb build the graph, vg aligns reads to it. This is worth the trouble mainly if you are chasing structural variation or genes with high copy-number variability.

Comparing bacterial genomes, including ones from your own microbiome

If your comparison targets come from a stool or oral metagenome, the first step is removing human reads. Do this before you assemble, before you upload anything, and before you share. Methods differ substantially in how many human reads they leave behind and how many microbial reads they wrongly discard, and alignment-based filtering against a complete human reference is the safer default 5. Map with bowtie2 or minimap2 against T2T-CHM13 and keep the unmapped pairs.

Then work down the resolution ladder:

  1. Distance. skani dist --qi --ri computes ANI in seconds per pair and is accurate down to roughly 80% ANI. mash sketch -s 10000 plus mash triangle is faster still for thousands of genomes. Above 95% ANI you are usually within a species.
  2. Quality and identity. CheckM2 for completeness and contamination, GTDB-Tk classify_wf for taxonomy against a consistent tree. Bins under ~50% complete are not worth comparing.
  3. Gene content. bakta or prokka to annotate, then panaroo -i *.gff3 -o pan --clean-mode strict -a core for the core/accessory split. Panaroo handles annotation error better than Roary, which matters when your input is MAGs rather than isolates. The core versus accessory question is the central one in bacterial comparative genomics, and it only became tractable once genome counts got large enough to estimate a pangenome curve 6.
  4. Strain-level. snippy or ska2 for SNP distances within a species, PopPUNK for clustering across a large collection. Typing schemes (cgMLST, wgMLST, SNP-based) each carry different assumptions about what counts as the same strain, and the choice changes your clusters 7.

Reference choice matters here more than tool choice. Cultured-isolate collections with matched genomes improved metagenomic classification substantially over relying on whatever was in public databases 8, and the accuracy of everything downstream tracks the reference set you classify against 9.

Cross-species comparison

For finding conserved non-coding elements, the pipeline is whole-genome alignment (LASTZ or Cactus) followed by conservation scoring (phastCons/phyloP) and a percent-identity plot. VISTA packaged this into a usable form and its output format (identity over a sliding window, with conserved regions shaded) is still the clearest way to look at a locus 1. Human-mouse comparison in particular identified regulatory elements that sequence inspection of the human genome alone did not reveal, which is the argument for doing it at all 10.

RefSeq or GenBank

GenBank is the submitter’s record: whatever was deposited, including duplicates, inconsistent annotation, and assemblies nobody curated. RefSeq is NCBI’s curated, non-redundant derivative with uniform annotation. For comparative work use RefSeq (GCF_ accessions) so annotation differences reflect biology rather than submitter habits. Use GenBank (GCA_) when the assembly you need has no RefSeq counterpart, which is common for recent MAGs, or when you specifically want the submitter’s original annotation.

One caution on interpretation: comparison tools tell you how sequences differ. They do not tell you what a difference means for a person. Anything you find in your own genome that looks clinically relevant needs a clinical geneticist and an orthogonally validated test, not a PAF file.

Questions people also ask

What is the best tool for comparing bacterial genomes? There isn’t one, but the sequence we’d use is skani for ANI, GTDB-Tk for taxonomy, panaroo for pangenome, and PopPUNK or ska2 for strain clustering. Use pairwise viewers like ACT or pyGenomeViz only when you are down to a handful of genomes and want to look at synteny directly.

What tools can I use to compare genomes? By layer: bcftools isec and hap.py/vcfeval for variants, minimap2 plus paftools.js or nucmer/dnadiff for assemblies, SyRI for structural rearrangements, pggb/minigraph-cactus/vg for pangenome graphs, VISTA or phastCons for cross-species conservation.

Is RefSeq or GenBank better? RefSeq for consistency, GenBank for coverage. If a genome exists in both, the GCF_ and GCA_ accessions share sequence and differ in annotation.

How do I compare my own genome to someone else’s? Normalize both VCFs against the same reference build with bcftools norm -m -any -f, then bcftools isec -c none. Expect roughly 4–5 million variant sites per genome against GRCh38, most of them common and shared, so the interesting output is the private set, not the count.

Why do two pipelines disagree on the same sample? Representation differences in indels, different filtering thresholds, different reference builds, and genuine ambiguity in hard regions. Comparing inside a benchmark high-confidence BED separates pipeline disagreement from region difficulty 2.

Oak builds longitudinal molecular profiles of individuals: whole-genome sequencing, RNA sequencing, proteomics, blood biomarkers, and continuous glucose data, integrated into one model of you. Build your profile.

Footnotes

  1. K. A. Frazer, L. Pachter, A. Poliakov, et al. VISTA: computational tools for comparative genomics. Nucleic Acids Research, 2004. https://doi.org/10.1093/nar/gkh458 ↩ ↩2

  2. Justin M. Zook, David Catoe, Jennifer McDaniel, et al. Extensive sequencing of seven human genomes to characterize benchmark reference materials. 2015. https://doi.org/10.1101/026468 ↩ ↩2

  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. Rachel M. Sherman, Steven L. Salzberg. Pan-genomics in the human genome era. Nature Reviews Genetics, 2020. https://doi.org/10.1038/s41576-020-0210-7 ↩

  5. Stephen J. Bush, Thomas R. Connor, Tim E.A. Peto, et al. Evaluation of methods for detecting human reads in microbial sequencing datasets. Microbial Genomics, 2020. https://doi.org/10.1099/mgen.0.000393 ↩

  6. Dawn Field, Gareth Wilson, Christopher van der Gast. How do we compare hundreds of bacterial genomes?. Current Opinion in Microbiology, 2006. https://doi.org/10.1016/j.mib.2006.08.008 ↩

  7. Laura Uelze, Josephine Grützke, Maria Borowiak, et al. Typing methods based on whole genome sequencing data. One Health Outlook, 2020. https://doi.org/10.1186/s42522-020-0010-1 ↩

  8. Samuel C. Forster, Nitin Kumar, Blessing O. Anonye, et al. A human gut bacterial genome and culture collection for improved metagenomic analyses. Nature Biotechnology, 2019. https://doi.org/10.1038/s41587-018-0009-7 ↩

  9. George M. Weinstock. Genomic approaches to studying the human microbiota. Nature, 2012. https://doi.org/10.1038/nature11553 ↩

  10. Len A. Pennacchio. Insights from human/mouse genome comparisons. Mammalian Genome, 2003. https://doi.org/10.1007/s00335-002-4001-1 ↩