Whole Genome Comparison: What You Can Compare, and With Which Tools

“Whole genome comparison” means four different jobs, and the tool depends on which one you want. Comparing your genome to the reference is variant calling (DeepVariant or GATK, output is a VCF). Comparing two callsets of the same sample is concordance evaluation (rtg vcfeval or hap.py). Comparing your genome to another person’s is relatedness and shared-segment analysis (somalier, KING, IBD segment callers). Comparing human sequence to mouse or chicken sequence is cross-species alignment, and as an individual you consume the output rather than produce it: precomputed phyloP and phastCons tracks from the UCSC multiple alignments. If you have your own WGS data and you are trying to learn something about yourself, you will spend almost all your time in the first three.

Comparing your genome to the reference

This is the comparison that produces your data. A 30x PCR-free Illumina genome aligned to GRCh38 yields roughly 4.5 to 5 million variant sites per person: about 3.5 million SNVs and 500,000 to 1,000,000 short indels, plus a few thousand structural variants that short reads can see. High-coverage sequencing of individual genomes has been used since the early 2010s to characterize this full spectrum of variation, including the indel and CNV classes that array genotyping misses entirely ¹.

The pipeline we would run:

bwa-mem2 mem -t 32 -K 100000000 -Y GRCh38_full_analysis_set.fa r1.fq.gz r2.fq.gz \
  | samtools sort -@ 8 -o sample.bam -
samtools index sample.bam

run_deepvariant --model_type=WGS --ref=GRCh38_full_analysis_set.fa \
  --reads=sample.bam --output_vcf=sample.vcf.gz --output_gvcf=sample.g.vcf.gz \
  --num_shards=32

Use the analysis set with ALT contigs and the decoy. Do not use a bare primary assembly. Reads from segmental duplications and from sequences absent in GRCh38 will otherwise pile onto the closest paralog and generate confident-looking false heterozygotes. Mapping is the step where most downstream differences are born, and different mappers make different tradeoffs between sensitivity in repeats and speed, which propagates directly into the variant set ².

Then normalize before you compare anything to anything:

bcftools norm -f GRCh38_full_analysis_set.fa -m -any -Oz -o sample.norm.vcf.gz sample.vcf.gz
bcftools index -t sample.norm.vcf.gz

-m -any splits multi-allelic records, and left-alignment against the reference fixes indel representation. Skip this and half your “differences” between two callsets will be the same indel written two ways.

Comparing two callsets of the same person

You will do this more often than you expect: two labs, two pipelines, GRCh37 versus GRCh38, or DeepVariant versus GATK on the same BAM. Do not use bcftools isec as your primary tool here. The command matches on position and allele string, so a 5 bp deletion in a homopolymer that one caller placed at 12345 and another at 12347 counts as two differences.

Use RTG Tools, which compares by haplotype equivalence:

rtg format -o GRCh38.sdf GRCh38_full_analysis_set.fa
rtg vcfeval -b callsetA.vcf.gz -c callsetB.vcf.gz -t GRCh38.sdf \
  --evaluation-regions=confident.bed --ref-overlap --vcf-score-field=GQ -o eval/

hap.py --engine=vcfeval wraps the same comparison and gives you a stratified summary (SNV versus INDEL, transition/transversion, by genomic context). Restrict to high-confidence regions with --evaluation-regions. Genome-wide, two good pipelines on the same 30x BAM will agree on more than 99% of SNVs and noticeably less on indels, with the disagreement concentrated in homopolymers and short tandem repeats. That concentration is the signal: if your two callsets disagree outside low-complexity sequence, something is wrong with one of them.

For a quick sanity check that two files are the same person at all, somalier extract then somalier relate takes about a minute and compares genotypes at a few thousand informative sites.

Comparing your genome to another person’s

A parent, a sibling, a child. The comparison that is well defined here is kinship and identity by descent, not “how many variants differ,” which is a number dominated by common polymorphism.

somalier extract -d extracted/ -s sites.GRCh38.vcf.gz -f GRCh38.fa sampleA.bam sampleB.bam
somalier relate --infer extracted/*.somalier

KING-robust kinship coefficients fall into well-separated bands: above 0.354 means duplicate or monozygotic twin, 0.177 to 0.354 first degree (parent-child or full sibling), 0.0884 to 0.177 second degree. Parent-child and full siblings share the same expected kinship but differ in the fraction of sites with zero shared alleles, which is near zero for parent-child and about 0.25 for siblings. That statistic is what separates them.

If you want population-level context (ancestry, admixture proportions and deep ancestry inference), whole-genome data supports much finer reconstruction than array genotypes. The reason is that the rare and low-frequency variants that carry most of the demographic signal are only visible in sequencing ³. Run PCA against 1000 Genomes reference panels with plink2 --pca after LD pruning, or use ADMIXTURE with a supervised reference set.

Phasing matters if you are asking whether two variants in the same gene sit on the same haplotype. Read-backed phasing from short reads (WhatsHap, HapCUT2) gives you accurate but short blocks. Statistical phasing against a large reference panel (SHAPEIT, Eagle, Beagle) gives you chromosome-scale blocks with switch errors scattered through them. Systematic comparison of these strategies on whole human genomes shows trio-based and read-based approaches each dominate on different metrics, and combining read data with a population panel outperforms either alone ⁴. If you have parental sequence, trio phasing is the accurate option and you should use it.

Comparing human to other species

You are not going to run a multiple genome alignment yourself, and you do not need to. The useful artifact is a per-base conservation score. Download hg38.phyloP100way.bw and hg38.phastCons100way.bw from UCSC and query them with bigWigAverageOverBed against a BED of your variant positions. Cross-species sequence comparison identifies functional non-coding elements that human sequence alone does not distinguish from background, which is the main reason to bother with it ⁵. The comparison remains a filter rather than an answer: the algorithmic problems of aligning large eukaryotic genomes, with rearrangement, duplication, and lineage-specific insertion, were laid out clearly two decades ago and the hard parts have not disappeared ⁶.

Practically: a variant of unknown significance sitting at phyloP > 4 in a deeply conserved element is more interesting than one at phyloP ≈ 0. This is a prioritization heuristic for your own reading. It is not a clinical interpretation. Any variant you are inclined to act on needs confirmation in a CLIA-certified laboratory and discussion with a clinical geneticist or genetic counselor.

Assembly versus assembly

If you have a long-read de novo assembly (HiFi plus Hi-C, or ONT ultra-long), compare it to the reference directly rather than through a VCF:

minimap2 -cx asm5 --cs -t 32 GRCh38.fa sample.asm.fa > aln.paf
paftools.js call -L 5000 -f GRCh38.fa aln.paf > asm_calls.vcf

For a diploid assembly use dipcall, which handles both haplotypes and emits a confident-region BED. For dotplot-level views of rearrangement, nucmer --maxmatch followed by dnadiff from MUMmer4 is still the fastest way to see structure. Assembly-based comparison is how you see the ~1-2% of the genome that short reads never resolve: large insertions, VNTR expansions, the acrocentric arms.

Failure modes that produce fake differences

Contig naming (chr1 versus 1), which silently produces zero overlaps in bedtools. Mixed builds, since GRCh37 and GRCh38 coordinates differ by megabases in places. Convert with CrossMap.py vcf and a real chain file, and expect a small percentage of sites to fail to lift. Unnormalized indels. Comparing a filtered VCF against an unfiltered one without checking the FILTER column. Regions of the reference with MAPQ 0 pileups, where every caller produces garbage: exclude them with the ENCODE blacklist or restrict to GIAB high-confidence regions. And coverage: below about 20x, heterozygous sites start dropping out, and a comparison between a 15x and a 40x genome will show a deficit of heterozygotes in the shallower one that has nothing to do with biology. Long-read and short-read platforms also differ systematically in which variant classes they recover, so cross-platform comparisons need stratification before they mean anything ⁷.

Questions people also ask

What tool can I use to compare two genomes? For two human samples represented as VCFs, rtg vcfeval for haplotype-aware concordance and somalier for relatedness. For two assemblies, minimap2 -cx asm5 plus paftools.js call, or dipcall for diploid assemblies. For two bacterial or viral genomes, MUMmer4 (nucmer/dnadiff) or Mash for a quick distance estimate.

What are the most common tools in comparative genomics? The aligners are minimap2, LASTZ, MUMmer and Cactus for multi-species. Alongside them are the UCSC chain/net toolchain, conservation scoring with phyloP and phastCons via PHAST, and orthology tools like OrthoFinder. For human-to-human work the variant-centric stack is bcftools, RTG Tools, hap.py, PLINK and KING.

Can you give an example of comparative genomics? Scoring every base of the human genome by conservation across 100 vertebrates and using that score to pick out non-coding regulatory elements that show no signal in human sequence alone ⁵.

What are the applications of comparative genomics? Annotating functional elements, prioritizing variants of unknown significance, reconstructing population history and demography from shared segments and rare variant sharing ³, and tracking pathogen evolution.

How many differences are there between two human genomes? Roughly 4 to 5 million variant sites relative to GRCh38 per individual, of which about 3.5 million are SNVs, with the remainder indels and larger structural variants ¹. Two unrelated people differ at a comparable order of magnitude, and most of those differences are common polymorphisms shared with millions of others.

Footnotes

1. Hui Shen, Jian Li, Jigang Zhang, et al. Comprehensive Characterization of Human Genome Variation by High Coverage Whole-Genome Sequencing of Forty Four Caucasians. PLoS ONE, 2013. https://doi.org/10.1371/journal.pone.0059494 ↩ ↩²

2. Sophie Schbath, Véronique Martin, Matthias Zytnicki, et al. Mapping Reads on a Genomic Sequence: An Algorithmic Overview and a Practical Comparative Analysis. Journal of Computational Biology, 2012. https://doi.org/10.1089/cmb.2012.0022 ↩

3. Krishna R. Veeramah, Michael F. Hammer. The impact of whole-genome sequencing on the reconstruction of human population history. Nature Reviews Genetics, 2014. https://doi.org/10.1038/nrg3625 ↩ ↩²

4. Yongwook Choi, Agnes P. Chan, Ewen Kirkness, et al. Comparison of phasing strategies for whole human genomes. PLOS Genetics, 2018. https://doi.org/10.1371/journal.pgen.1007308 ↩

5. Len A. Pennacchio, Edward M. Rubin. Comparative genomic tools and databases: providing insights into the human genome. Journal of Clinical Investigation, 2003. https://doi.org/10.1172/jci17842 ↩ ↩²

6. Webb Miller. Comparison of genomic DNA sequences: solved and unsolved problems. Bioinformatics, 2001. https://doi.org/10.1093/bioinformatics/17.5.391 ↩

7. Jianhua Zhao, Struan F.A. Grant. Advances in Whole Genome Sequencing Technology. Current Pharmaceutical Biotechnology, 2011. https://doi.org/10.2174/138920111794295729 ↩