Bio Variant Normalization

Version Compatibility

Reference examples tested with: bcftools 1.19+, cyvcf2 0.30+

Before using code patterns, verify installed versions match. If versions differ:

• Python: pip show <package> then help(module.function) to check signatures
• CLI: <tool> --version then <tool> --help to confirm flags

The --atomize/--old-rec-tag flags require bcftools 1.12+ (--keep-sum requires 1.11+). Earlier versions require vt decompose_blocksub as an alternative.

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

Variant Normalization

Left-align indels, decompose MNPs, and split multiallelic sites using bcftools norm.

"Put my VCF in canonical form before comparing, annotating, or merging" -> Enforce one representation per biological event so tuple-keyed operations recognize the same variant across sources.

• CLI: bcftools norm -m-any -f ref.fa (split multiallelic, then left-align + parsimony each biallelic)
• Alternatives: vt normalize + vt decompose; GATK LeftAlignAndTrimVariants

The Normalized Form (governing principle)

A variant is normalized if and only if it is BOTH (Tan, Abecasis & Kang 2015):

1. Parsimonious (minimal representation). No base can be trimmed from either end of REF/ALT without creating a zero-length allele or changing the variant. Indels keep exactly ONE anchor base (VCF forbids empty alleles), so a deletion of A is written REF=CA ALT=C, not REF=A ALT=.
2. Left-aligned. POS cannot be shifted further left while keeping the allele lengths and the reference-relative meaning unchanged.

Both are required: parsimony alone leaves an indel positionally ambiguous inside a repeat; left-alignment alone can leave redundant shared bases. Together they define a single canonical (CHROM, POS, REF, ALT) key, and normalization is idempotent -- normalizing an already-normalized VCF changes nothing. That idempotence is what makes the key safe to join, merge, and match on.

The algorithm (as in vt normalize and bcftools norm):

• Right-trim: while all alleles end in the same base AND all have length >= 2, drop the last base of each.
• Left-extend: while any allele has length 0, prepend the reference base immediately to the left to every allele and decrement POS.
• Left-trim: while all alleles start with the same base AND all have length >= 2, drop the first base of each and increment POS.

The right-trim-then-left-roll loop is exactly what shifts an indel through a repeat to its leftmost equivalent position.

When Normalization is Mandatory

Not normalizing before certain operations leads to missed matches and false discordance. Normalization is required:

• Before comparing variants from different callers. Each caller may represent the same indel at different positions or encode MNPs differently. Without normalization, identical variants appear discordant.
• Before database annotation. dbSNP, ClinVar, and gnomAD store variants in canonical left-aligned, parsimonious representation. A right-aligned or non-parsimonious indel will fail to match its database entry.
• Before merging VCF files from different sources. bcftools merge matches on CHROM/POS/REF/ALT; different representations of the same variant produce duplicate entries instead of a single merged record.
• Before any variant set operations. Intersection (bcftools isec), complement, and union operations all rely on exact positional matching. Non-normalized variants silently fall through set comparisons.

Normalization is generally safe to skip only when a single caller produced all variants and no cross-file comparison or database lookup is needed.

Why Left-Alignment Matters in Repeats (the silent miss)

The same variant can be written multiple ways:

chr1  100  ATCG  A      (right-aligned)
chr1  100  ATC   A      (left-aligned, parsimonious -- the canonical form)
chr1  101  TCG   T      (shifted position, different anchor base)

The ambiguity is worst in homopolymers and tandem repeats. In reference ...AAAAAAA..., a single-base A deletion is positionally ambiguous: deleting ANY one of the seven A's yields the identical alternate haplotype, so different callers/aligners emit it at different POS. Left-alignment defines the canonical (leftmost) position, giving the one biological event one representation.

The load-bearing consequence: a one-base-off (non-left-aligned) indel in a homopolymer is a valid VCF line that produces a different (CHROM, POS, REF, ALT) tuple. Annotation and matching keyed on that tuple then silently miss the dbSNP/ClinVar/gnomAD entry that sits at the canonical coordinate -- a clinically actionable variant is reported as absent or novel, with no error thrown. The VCF is structurally valid; the numbers are wrong. bcftools merge of a normalized and an un-normalized cohort likewise emits duplicate rows for the same event, inflating counts and splitting allele frequencies.

Decision: always left-align + parsimony (bcftools norm -f ref.fa) against the EXACT downstream reference before annotation, database matching, set operations, or merging -- never rely on callers to emit canonical form.

Recommended Normalization Pipeline

The order of operations matters. Performing these steps out of order can produce incorrect results (e.g., left-aligning a multiallelic record may normalize differently than splitting first, then left-aligning each biallelic record independently).

The correct order:

1. Decompose MNPs into atomic SNPs (--atomize)
2. Split multiallelic sites into biallelic records (-m-)
3. Left-align and trim against the reference (-f reference.fa)

Combined as a piped pipeline:

bcftools norm --atomize input.vcf.gz | \
    bcftools norm -m- | \
    bcftools norm -f reference.fa -Oz -o normalized.vcf.gz
bcftools index normalized.vcf.gz

For VCFs without MNPs (e.g., GATK HaplotypeCaller output, which does not emit MNPs), the atomize step can be skipped:

bcftools norm -m- input.vcf.gz | \
    bcftools norm -f reference.fa -Oz -o normalized.vcf.gz

A single-pass bcftools norm -f ref.fa -m-any is acceptable for basic use cases but does not control the decomposition order and skips MNP atomization.

Tool Discordance: bcftools vs vt vs GATK

All three implement the same Tan-2015 left-align + parsimony core and agree on simple biallelic indels. They diverge on decomposition:

The discordance that burns people: vt decompose_blocksub splits MNPs/block substitutions into SNPs by default, bcftools norm does not (it needs --atomize). Running the same normalization with the two tools yields a different variant count (bcftools keeps an MNP as one record; vt emits separate SNPs). If cohort A is atomized with vt and cohort B is left un-atomized with bcftools, every MNP is a systematic representation mismatch, manufacturing spurious cohort-private "variants" in any cross-cohort comparison.

Decision: standardize on ONE normalization tool + exact flag set across every cohort intended for comparison, and record the command. bcftools norm is the de-facto production standard (htslib-maintained; integrates split, atomize, and left-align in one pass). Note GATK's left-alignment window (--max-leading-bases) means indels inside repeats longer than the window may not fully left-align -- verify the installed default with gatk LeftAlignAndTrimVariants --help before trusting STR-embedded indels.

Left-Alignment

"Normalize my VCF before comparing callers" -> Left-align indel representations and split multiallelic sites for consistent variant comparison.

bcftools norm -f reference.fa input.vcf.gz -Oz -o normalized.vcf.gz

Left-alignment cannot roll an indel leftward without the reference bases to its left, so -f/--fasta-ref is non-negotiable. Two traps make the reference choice load-bearing:

• The FASTA must be the exact reference the VCF was called against. A different build patch, chr1-vs-1 contig naming, or a masked-vs-unmasked sequence yields wrong left-aligned coordinates. On a REF mismatch bcftools norm errors and exits non-zero by default (-c e) rather than warning -- do not paper over it with -c w; a single off-by-one in a contig shifts every indel.
• It must be the SAME reference used downstream (the annotation-database build, the other cohort). Normalizing to GRCh38 and matching against a GRCh37 dbSNP is a guaranteed miss even when local left-alignment is perfect.

Check for Normalization Issues

bcftools norm -f reference.fa -c w input.vcf.gz > /dev/null

Check modes (-c):

• e - Error and exit on mismatch (default)
• w - Warn on mismatch and continue (use this to enumerate all mismatches)
• x - Exclude mismatches
• s - Set correct REF from reference

Multiallelic Splitting

Split Multiallelic to Biallelic

bcftools norm -m-any input.vcf.gz -Oz -o split.vcf.gz

Before:

chr1  100  .  A  G,T  30  PASS  .  GT  1/2

After:

chr1  100  .  A  G  30  PASS  .  GT  1/0
chr1  100  .  A  T  30  PASS  .  GT  0/1

Splitting Caveats

Splitting creates artificial missing information. A sample with genotype 1/2 (compound heterozygous for two different ALT alleles) becomes 0/1 in both split records. The information that both alleles were present at the same site in the same individual is lost. This has consequences for:

• Phasing and compound heterozygosity detection. Clinical pipelines that identify compound hets (two damaging variants on different alleles of the same gene) can misinterpret split records as independent heterozygous calls rather than co-occurring alleles at one site.
• Allele depth (AD) interpretation. AD values are retained per allele in each split record, but the genotype relationship between alleles at the same site is gone.
• Population allele frequency estimation. Splitting followed by naive frequency calculation can double-count samples at multiallelic sites.

Decision guidance:

When a downstream tool does not require splitting, prefer keeping multiallelic sites intact to preserve genotype relationships.

Split Options

Join Biallelic to Multiallelic

bcftools norm -m+any input.vcf.gz -Oz -o merged.vcf.gz

Rejoining after analysis can restore compound heterozygosity context, but only if the split records were not independently filtered (removing one allele of a 1/2 site makes the remaining record misleading).

Field Reapportionment and Spanning Deletions

Splitting a multiallelic is NOT information-lossless. Per-allele fields must be reapportioned according to their header Number code, and bcftools norm -m- uses those codes to subset correctly:

• Number=A (one value per ALT, e.g. AF, AC) -> take the k-th element for the k-th ALT.
• Number=R (one per allele including REF, e.g. AD -> ref depth first) -> off-by-one relative to A.
• Number=G (one per genotype, e.g. PL, GL) -> subsetting needs the genotype-index formula.
• Number=. -> variable count; parsers cannot auto-subset, so these fields are silently carried whole onto every split record.

The trap: a custom INFO/FORMAT field mis-declared Number=. when it is really A is NOT subset on split, so every biallelic record keeps the full multiallelic vector and downstream tools read the wrong allele's value. Joining back (-m+) cannot always reconstruct the original PLs exactly. Rule: split once, early, and stay biallelic; join only for final delivery if a consumer requires it.

Spanning-deletion * alleles (VCF: "allele missing due to overlapping deletion") are meaningful only relative to the overlapping deletion recorded on another line. Splitting can strand a * allele from the deletion it references; bcftools handles this, but naive third-party splitters corrupt it. Use --keep-sum AD when the summed allele depth must be preserved across split records.