Sanity check · Variant Calling
How to Catch a Genome Build Mismatch in Variant Calling
A BED file from the wrong build gives you confident, wrong overlaps, and the header check that catches it takes thirty seconds.
By Ming "Tommy" Tang, Director of Bioinformatics in Big Pharma · Reviewed September 2026 · 5 min read
You inherited a pipeline, or you're merging a capture panel BED from three years ago with a freshly aligned hg38 cohort. The variant caller runs clean, the VCF has no red errors in the log, and you move straight into annotation and interpretation. Nothing about the run tells you that half your files are on a different genome build than the other half, both hg19 and hg38 are valid FASTA files, and every tool in the chain will happily call variants against whichever one you point it at.
The cost shows up later, and it doesn't look like a crash. hg19 and hg38 coordinates diverge by anywhere from kilobases to megabases across the genome, not by a fixed offset. A truly pathogenic variant gets flagged as "outside target region" because your BED file is on the wrong build. A benign polymorphism gets reported as sitting inside a hotspot it never touched. In a clinical or translational pipeline, that's not a bioinformatics inconvenience, it's a wrong call that reaches a report.
This page gives you the checks to run in the next hour: header comparisons, a mitochondrial length check, and a landmark-gene spot check, in order of how cheap they are to run. None of them require re-aligning anything. They just tell you, for every file already sitting in your pipeline, which build it's actually on.
What it looks like when it's happening
- GATK throws an explicit error naming two input files and reporting the same contig with two different lengths
- IGV loads a BAM or VCF and a gene you expect to see lands in the wrong place, or in an intergenic desert, relative to the genome build selected in the browser
- A capture panel BED file's coordinates overlap introns or intergenic regions instead of the exons the panel was designed to target
- The VCF header shows ##contig=<ID=chr1,length=249250621> in one file and length=248956422 in another file from the same sample
- A well-known common variant comes back with allele frequency 0 or absent from gnomAD/ClinVar annotation, even though it should be a documented hit
- Mitochondrial variant positions are consistently off by roughly 2 base pairs relative to published haplogroup marker positions
- One file's header uses chr-prefixed contig names (chr1) while another in the same pipeline uses bare numbers (1), and a join or overlap between them silently returns zero or nonsensical results
- A liftover run produces an unusually large .unmap file, with entries tagged Fail(KeyError) pointing at contigs the target reference doesn't recognize
Why it happens
GRCh38 fixed thousands of misassemblies present in GRCh37, filled gaps, and added alternate-haplotype and decoy sequences that GRCh37 never had. Every contig's coordinate system shifted as a result, but not by a uniform offset you could add or subtract, the shift is idiosyncratic per region, which is exactly why liftover requires a chain file and a dedicated tool rather than simple arithmetic. There is no shortcut that lets you eyeball a coordinate and know which build it belongs to; you have to check the file's actual header or a known landmark.
Naming convention makes this worse because it travels independently of the build itself. hg19 uses chr-prefixed contig names (chr1), while b37, which carries the identical GRCh37 coordinates, uses bare numbers (1). hg38's canonical release again uses chr-prefixed names. Pipelines built around 1000 Genomes resources tend to inherit b37's bare naming; pipelines built around UCSC or more recent GATK bundles tend to use chr-prefixed names. A file can look mismatched purely on naming while actually being on the same coordinates, or it can look consistent on naming while sitting on a genuinely different build. You cannot tell which from the contig name alone; you need the length.
In practice the mismatch enters through the joints between tools, not through any single tool's internals. A shared reference directory on a cluster gets updated for new projects, but an inherited script still hardcodes a path to the old FASTA. A vendor's capture panel BED was designed against hg19 and never relifted when the group standardized on hg38. An annotation database like ClinVar or gnomAD gets pulled in whatever build is the current release, while a legacy germline call set from years ago sits untouched in an older one. Each individual file is internally consistent; the mismatch only exists in the comparison between files.
GATK does check contig header lengths for the files you feed it directly at call time and will fail loudly if two of its own inputs disagree. That check does not extend to your BED file, your annotation source, or the human doing manual review in IGV. Those are exactly the places a mismatch survives, because nothing in the software stack enforces build consistency across the whole workflow, only within the narrow set of files a given tool invocation touches directly.
The checks
Run them in order. Each one tells you what healthy looks like and what the problem looks like.
0/7 checked · saved in this browser
Pull the @SQ lines from the BAM header and the ##contig lines from the VCF header, and diff the lengths for the same chromosome name across BAM, VCF, and reference FASTA/.dict.
bashsamtools view -H sample.bam | grep "^@SQ" bcftools view -h sample.vcf.gz | grep "^##contig" grep "^>" reference.fasta | head -5- Healthy
- chr1 (or 1) reports the identical length in the BAM header, the VCF header, and the reference FASTA: 249250621 for hg19/GRCh37, 248956422 for hg38/GRCh38.
- Red flag
- The same contig name reports different lengths in two of the three files, this is the literal condition GATK's own "incompatible contigs" error is built to detect.
Check whether contig names are chr-prefixed (chr1) or bare (1) in each file that will be joined or overlapped with another.
bashsamtools view -H sample.bam | grep "^@SQ" | head -1 zcat panel.bed.gz | head -1 bcftools view -h sample.vcf.gz | grep "^##contig" | head -1- Healthy
- All files use the same naming convention, whatever it is.
- Red flag
- BAM/VCF use chr-prefixed names while the BED panel uses bare numbers (or vice versa), bcftools and BED-based overlap tools either error out or silently return zero matches depending on how strict they are, and a silent zero looks identical to "this sample has no variants in the panel."
Run the variant-calling step (HaplotypeCaller, Mutect2, or joint genotyping) and read the log output before assuming a clean exit means the files matched.
- Healthy
- The tool proceeds without complaint, meaning every file it directly read shares the reference's contig dictionary.
- Red flag
- The run halts with an error naming two specific input files and reporting the same contig with two different lengths, this check is free because GATK already performs it, but it only covers the files passed directly into that command, not your BED or annotation sources.
Extract the header line for the MT or chrM contig from the BAM or VCF and read off its declared length in base pairs.
bashsamtools view -H sample.bam | grep -E "SN:(MT|chrM)\b"- Healthy
- 16569 bp if the file uses the revised Cambridge Reference Sequence (rCRS), the MT sequence used by b37 and GRCh38.
- Red flag
- 16571 bp (GRCh37's MT, which is 2 bp longer than rCRS) turning up in a file you believe is hg38, or MT variant positions consistently off by about 2 bp from published haplogroup marker coordinates.
Pick a gene you know well, look up its coordinate range on UCSC or Ensembl for the specific build you believe your file is on, and confirm a variant you expect inside that gene actually falls within that range in your VCF.
- Healthy
- The variant's position falls inside the gene body at the coordinate range the browser reports for your intended build.
- Red flag
- The position falls outside the gene entirely, or lands inside a completely different gene, a fast tell that the file is on a different build than the one you looked up.
Count and inspect the contig names in the reference and BAM header for ALT-suffixed or decoy sequences.
bashgrep "^>" reference.fasta | wc -l samtools view -H sample.bam | grep "^@SQ" | grep -Ei "_alt|decoy"- Healthy
- For GRCh37 variant calling, hs37d5 with its decoy sequences present, per the recommended flavor for better variant-calling sensitivity. For GRCh38, the no_alt_analysis_set with no ALT-suffixed or unplaced-random contigs, which avoids mapper confusion.
- Red flag
- A GRCh38 BAM header full of *_alt or *_random contigs (the wrong flavor for a mapper-friendly analysis), or a GRCh37 pipeline running against a decoy-free FASTA when the project assumed hs37d5-level sensitivity.
Run the liftover and inspect the resulting .unmap file size and content relative to the input.
bashCrossMap.py vcf hg19ToHg38.over.chain.gz input.vcf reference_hg38.fa output.vcf wc -l output.vcf.unmap- Healthy
- A failure rate in the same range as CrossMap's own reported benchmark, about 0.19% of variants (1,567,838 of 1,570,748 converted successfully going GRCh37 to GRCh38).
- Red flag
- A .unmap file with a meaningfully higher fraction of variants failing, especially with Fail(KeyError) tags, usually the wrong chain file, or an input VCF that was already partly on the target build.
What to do about it
Fix a naming-only mismatch with reheader
When: The contig lengths already agree across files (same build, same coordinates) and the only difference is chr-prefixed vs bare naming.
Use samtools reheader to swap the header's SN: lines to match the naming convention the rest of your pipeline expects.
Caveat: Reheader replaces header text only; it does not verify or reorder coordinates. Confirm lengths match first (check #1) or you will relabel a genuine build mismatch as a cosmetic naming fix and hide the real problem.
Fix a contig order mismatch with Picard
When: Same build, same contig names, but the sort order differs between a BAM and the reference dictionary, common after merging files produced by different pipelines.
Run Picard ReorderSam on the BAM against the target reference, and Picard SortVcf with the reference's .dict for the VCF.
Caveat: Requires a matching .dict file for the target reference alongside the FASTA. This does not fix actual coordinate differences between builds, only ordering, it is the wrong tool if the contigs disagree in length.
Liftover with CrossMap for a genuine build mismatch
When: Variants were called natively on hg19/GRCh37 but need to sit alongside an hg38-based cohort or annotation database.
Run CrossMap.py vcf with the correct chain file and target reference FASTA, then inspect the .unmap output (check #7) before trusting the converted file, and spot-check a handful of converted positions in IGV against the new build.
Caveat: Liftover is lossy near breakpoints: indels, structural variants, and MT variants (2 bp offset) often fail to convert cleanly or convert to slightly wrong positions. Prefer re-calling from FASTQ against the correct reference for anything you plan to keep using long-term; reserve liftover for one-off comparisons.
Re-download the recommended reference flavor and re-align
When: You're nominally on the right build but the wrong flavor of it, b37 instead of hs37d5 for GRCh37, or a reference with ALT contigs instead of no_alt_analysis_set for GRCh38.
Pull hs37d5 for GRCh37 work or the no_alt_analysis_set FASTA for GRCh38 work, and re-align from FASTQ against it.
Caveat: This is a full re-alignment, not a header fix, budget for the compute and storage cost across a whole cohort before promising a quick turnaround. Reheadering an existing BAM does not relocate reads that would have mapped differently against decoy or ALT sequences.
Liftover the BED file, not just the VCF
When: A capture panel, blacklist, or gene-list BED predates or postdates your team's switch to a new build.
Liftover the BED itself with CrossMap or liftOver, then re-run the landmark gene check (check #5) on the lifted BED before using it to filter or annotate calls.
Caveat: BED files use 0-based coordinates while GFF/GTF and VCF use 1-based coordinates. A naive coordinate copy across that boundary introduces a new off-by-one error even after the build itself is correctly converted.
When not to "fix" it
Not every apparent mismatch is a bug. A legacy hg19 germline call set kept for continuity with a clinical trial's original annotation, sitting next to a new hg38 tumor cohort for a different study, does not need reconciling, each stays valid on the build it was called on, and any comparison between them should go through an explicit, logged liftover step, not a forced retroactive conversion of the historical set just to make filenames look consistent. Forcing an old cohort onto a new build purely for cosmetic uniformity risks introducing exactly the liftover artifacts (indel shifts, MT offset, failed conversions) that a stable legacy file set doesn't have. If your analysis never leaves one build internally, don't spend cycles lifting over a reference panel that will never be compared across builds.
Five things experienced analysts do here
- Put the build name in the filename itself, not just in documentation, sample.hg38.bam rather than sample.bam, because headers get checked deliberately but filenames get trusted at a glance, and that glance should be correct.
- Bake the contig-length header check into your pipeline as an assertion that runs before variant calling starts, so a reference downloaded to the wrong path fails loudly in seconds instead of silently miscalling a cohort.
- Standardize on one reference flavor per build across the whole team, hs37d5 for GRCh37, no_alt_analysis_set for GRCh38, and store the exact FASTA and .dict path in a shared config file, not a wiki page someone forgets to update.
- When you inherit someone else's VCF, read the ##contig header lines before you load it in IGV, IGV will display it against whatever genome build you have selected and show confident, wrong tracks if you don't set the build first.
- Treat a liftover's .unmap output as a QC metric you log every time, not a nuisance file you delete, investigate any run with a meaningfully higher failure rate than the tool's documented benchmark.
Questions people ask
- How do I tell if my BAM file is hg19 or hg38?
Run
samtools view -H sample.bam | grep "^@SQ"and read the length for chr1. A length of 249250621 means hg19/GRCh37; 248956422 means hg38/GRCh38. Do this before you trust any downstream annotation, because the FASTA header carries no obvious build label otherwise.- Can I just rename chr1 to 1 to fix a build mismatch?
Only if the underlying coordinates already agree and the only difference is naming convention. Renaming a contig with
samtools reheaderchanges the label, not the sequence positions, so if you rename without first confirming the contig lengths match, you paper over a real coordinate mismatch instead of fixing one.- What's the difference between hg19 and b37 reference FASTAs?
They carry the same GRCh37 coordinates, but hg19 uses chr-prefixed names (chr1) while b37 uses bare numbers (1). The one substantive difference beyond naming is the mitochondrial contig: GRCh37's MT is 2 base pairs longer than the revised Cambridge Reference Sequence (rCRS) used in b37/hg38, so MT variant positions can shift by that amount if you mix them.
- Do I need to liftover my whole cohort to hg38?
Only if you need it to sit alongside an hg38-native cohort or annotation database. If you're not merging or comparing across builds, leave a completed hg19 call set alone. When you do need to move to hg38, re-calling from FASTQ against the correct reference is more accurate than lifting over an existing VCF, because liftover cannot recover reads that mapped differently under the new assembly.
- How accurate is liftover with CrossMap?
CrossMap reports about 99.81% successful conversion for GRCh37 to GRCh38 variant coordinates in its own benchmark (1,567,838 of 1,570,748 variants). The remainder lands in the
.unmapoutput file, which you should inspect rather than discard, since a spike in failures usually means you used the wrong chain file or the input was already partly on the target build.
Related pages
- Guide · How to Catch a Sample Swap in Variant Calling
- Guide · How to Detect Contamination in Variant Calling
- Glossary · Variant calling
Related reading on the blog
Sources
- The Most Common Stupid Mistakes In Bioinformatics — 0-based vs 1-based coordinate mistakes compounding build mismatches
- Errors about input files having missing or incompatible contigs, GATK — GATK's own contig-length error behavior
- Human genome reference builds - GRCh38 or hg38 - b37 - hg19, GATK — GATK build support notes
- CrossMap Documentation — CrossMap.py vcf usage and .unmap output
- CrossMap GitHub Repository — 99.81% GRCh37-to-GRCh38 conversion accuracy figure
- VCF 4.1 Specification — ##contig header line format
Part of the Genome build mismatches series.