Variant Normalization¶

How gbcms prepares variants before counting — validation, left-alignment, and homopolymer decomposition.

Overview¶

Before counting any reads, every variant passes through the preparation pipeline in prepare_variants(). This ensures consistent, biologically correct coordinates regardless of how the input was generated.

flowchart LR
    Input(["📄 Input Variant"]):::start --> S1["Step 1: MAF Anchor"]
    S1 --> S2["Step 2: REF Validation"]
    S2 -->|"REF_MISMATCH"| Fail(["❌ Rejected"]):::fail
    S2 -->|"PASS*"| S3["Step 3: Left-Alignment"]
    S3 --> S4["Step 4: ref_context Fetch"]
    S4 --> S5["Step 5: Homopolymer Detection"]
    S5 --> Output(["✅ PreparedVariant"]):::pass

    classDef start fill:#9b59b6,color:#fff,stroke:#7d3c98,stroke-width:2px;
    classDef step fill:#3498db,color:#fff,stroke:#2471a3,stroke-width:2px;
    classDef fail fill:#e74c3c,color:#fff,stroke:#c0392b,stroke-width:2px;
    classDef pass fill:#27ae60,color:#fff,stroke:#1e8449,stroke-width:2px;
    class S1,S2,S3,S4,S5 step;

Each input variant produces exactly one PreparedVariant, even if validation fails — this guarantees output row count always matches input.

Step 1: MAF Anchor Resolution¶

MAF files represent indels using - dashes. gbcms converts them to VCF-style anchor-based coordinates, which requires fetching the anchor base from the reference FASTA.

Step 2: REF Validation¶

The REF allele is compared against the reference genome at the stated position. An exact match passes immediately. If the exact match fails, a similarity score is computed — variants with ≥90% base match are corrected to the FASTA REF and proceed to counting.

flowchart TD
    Fetch["Fetch ref bases at pos"] --> Exact{"Exact match?"}
    Exact -->|"Yes"| Pass(["✅ PASS"]):::pass
    Exact -->|"No"| Sim["Compute similarity\n(matching bases / max length)"]
    Sim --> Thresh{"≥ 90%?"}
    Thresh -->|"Yes"| Corrected(["⚠️ PASS + WARN_REF_CORRECTED\nFASTA REF used downstream"]):::warn
    Thresh -->|"No"| Fail(["❌ FAIL + REF_MISMATCH\n→ excluded from counting"]):::fail

    classDef pass fill:#27ae60,color:#fff,stroke:#1e8449,stroke-width:2px;
    classDef warn fill:#f39c12,color:#fff,stroke:#d68910,stroke-width:2px;
    classDef fail fill:#e74c3c,color:#fff,stroke:#c0392b,stroke-width:2px;

When a variant gets a WARN_REF_CORRECTED reason (verdict stays PASS), the MAF's REF allele is replaced with the FASTA sequence for all downstream steps (left-alignment, haplotype construction). The original MAF REF is preserved in original_ref for auditing.

Step 3: Left-Alignment¶

For indels and complex variants, gbcms applies bcftools-style left-alignment to ensure consistent positioning in repetitive regions.

Every reference window prep fetches (the left-align window, ref_context, the adaptive repeat scan, an indel's shift region and the event reference) is clamped to the contig: near a contig end it holds the bases that exist. The FASTA reader rejects a window that passes the end. So before 6.6.0, each window failed within its own reach of a contig end: left-alignment within about 100bp, ref_context within its padding, the shift region within 256bp and the event reference within 60bp. Exact fetches (REF validation, a MAF anchor) are not clamped: a REF running past the end fails validation.

Every variant that reaches REF validation gets a type label derived from its final alleles (after anchor resolution, REF correction and alignment), with the same rule the readers use. A row rejected before that (EMPTY_ALLELE, NON_SEQUENCE_ALLELE, ALT_EQUALS_REF, START_BEFORE_CONTIG: see QC Flags; a failed MAF anchor fetch) keeps the reader's label.

Counting never reads this label — reads are classified by the alleles — so it is what gbcms normalize reports, not an input to counting.

Step 4: ref_context Fetch (Adaptive)¶

For indels and complex variants, a flanking reference sequence is fetched around the (possibly normalized) position. This context serves two purposes:

1. Windowed indel detection — Safeguard 3 uses ref_context to verify that shifted indels delete/insert the expected reference bases
2. Phase 3 alignment — the context is used to build REF and ALT haplotypes for dual-haplotype alignment (WFA triage → PairHMM). See Phase 3: Alignment Fallback.

Adaptive Context Padding¶

By default, gbcms dynamically adjusts the context padding based on nearby tandem repeats. This prevents ambiguous SW alignments in repetitive regions where a fixed padding might consist entirely of the repeat motif.

Multi-Anchor Footprint Sweep¶

Historically, scanning only the variant's strict anchor base caused miscalculations for left-aligned INDELs, as the trailing INDEL mutations generated by BWA mapping often fell outside the default 5bp padding window.

The engine now performs a Multi-Anchor Footprint Sweep. It scans across the entire variant footprint (pos through pos + ref_len) looking for tandem repeats, and extracts the maximum necessary padding to clear the repeat entirely.

flowchart LR
    SweepFootprint["Sweep variant footprint\n(pos → pos + ref_len)"] --> RepFound{"Tandem repeat\nfound?"}
    RepFound -->|"Yes (span = S)"| Formula["padding = max(5, S/2 + 3)\ncapped at 50bp"]
    RepFound -->|"No"| Default["padding = 5bp (default)"]
    Formula --> Fetch["Fetch ref_context"]
    Default --> Fetch

Example: Poly-A Insertion¶

Variant: chr7:100 A→AAAAAAAAAAAAA (12bp poly-A insertion)
Footprint: chr7:100-101
Flanking: ...AAAAAAAAAA[anchor]AAAAAAAAAA...

Sweep Detected: homopolymer max span=20bp → padding = max(5, 20/2+3) = 13bp

Fixed padding=5:     AAAAA + REF/ALT + AAAAA  (all A's — ambiguous, BWA alignments drop)
Adaptive padding=13: GCTTAAAAA... + REF/ALT + ...AAAAATTGAC  (anchored)

Supported Repeat Types¶

Step 5: Homopolymer Decomposition Detection¶

Some variant callers merge nearby events in a homopolymer run into one complex variant with an inflated deletion, e.g. CCCCCC→T where the reads show a smaller change at the run's end. --rescue-homopolymer is deprecated in 6.7.0 and removed in 6.8.0: it is accepted, warns and is ignored, and gbcms builds and counts no corrected allele. Until 6.6.x, with the flag, gbcms also counted a corrected allele for such calls and reports whichever of the two more reads support.

Detection Criteria¶

All four conditions must be met:

flowchart LR
    C1{"C1: REF is homopolymer?\n(all same base, ≥3bp)"}
    C2{"C2: Net deletion?\n(ref_len > alt_len)"}
    C3{"C3: ALT last base ==\nnext ref base after span?"}
    C4{"C4: ALT last base ≠\nhomopolymer base?"}
    Detected(["🔍 Decomposition detected"]):::detect
    Skip(["⏭️ Skip — not a decomposition"]):::skip

    C1 -->|Yes| C2
    C2 -->|Yes| C3
    C3 -->|Yes| C4
    C4 -->|Yes| Detected

    C1 -->|No| Skip
    C2 -->|No| Skip
    C3 -->|No| Skip
    C4 -->|No| Skip

    classDef detect fill:#f39c12,color:#fff,stroke:#d68910,stroke-width:2px;
    classDef skip fill:#95a5a6,color:#fff,stroke:#7f8c8d,stroke-width:2px;

Corrected Allele Construction¶

When detected, a corrected ALT allele is built by keeping every base of the run but the last and replacing that last base with the base that follows the run. The corrected allele has the same length as REF:

Example:  REF = CCCCCC, ALT = T, next_ref_base = T

Step 1: homopolymer base = C
Step 2: Confirmed: alt_last (T) == next_ref_base (T), T ≠ C
Step 3: corrected_alt = C × (ref_len - 1) + alt_last = CCCCC + T = CCCCCT

Dual-Counting Flow¶

Historical (the flag is ignored in 6.7.0, removed in 6.8.0): with --rescue-homopolymer, the corrected variant is stored in PreparedVariant.decomposed_variant; without it, prep builds no twin. During counting, both the original and corrected alleles are counted independently. The result with the higher alt_count wins:

flowchart TD
    Input(["⚡ Decomposition detected"]):::warn
    Input -->|"in parallel"| CountOrig
    Input -->|"in parallel"| CountDecomp

    CountOrig["count (original: CCCCCC→T)"]
    CountDecomp["count (corrected: CCCCCC→CCCCCT)"]

    CountOrig --> Compare
    CountDecomp --> Compare
    Compare{"corrected AD > original AD?"}

    Compare -->|"Yes"| UseDecomp["Return corrected counts"] --> FlagWarn["reason +=\nWARN_HOMOPOLYMER_DECOMP"]:::warn
    Compare -->|"No"| UseOrig["Return original counts"] --> FlagPass["gbcms_status = PASS"]:::pass

    classDef warn fill:#f39c12,color:#fff,stroke:#d68910,stroke-width:2px;
    classDef pass fill:#27ae60,color:#fff,stroke:#1e8449,stroke-width:2px;

The corrected allele is scored as unique sequence (repeat_span 0). In RNA with a GTF it takes the original's gene strand, so --enforce-strandedness applies to it as to the original. WARN_HOMOPOLYMER_DECOMP is set per sample: a row carries it only when that sample's corrected allele won.

Validation Status Reference¶

Status is carried in two fields: a verdict, gbcms_status (exactly PASS or FAIL), and a |-separated reason list, gbcms_status_reason (empty for a clean PASS). Both appear in the MAF (two columns) and the VCF INFO (GS = verdict, GSR = reasons). | is used (never ;/,, which are VCF-unsafe), so the reason string is byte-identical in the MAF and the VCF.

Reasons stack: a PASS variant can carry WARN_REF_CORRECTED|WARN_HOMOPOLYMER_DECOMP.

Limitations¶

1. Left-alignment only — gbcms left-aligns but does not right-align. If your input is right-aligned, the windowed scan (±5bp) may still catch it, but left-alignment is recommended.
2. Dynamic normalization window — Starting at 100bp, the window automatically doubles (up to 2500bp) when a variant shifts to the window edge. In practice, this covers even variants in centromeric/telomeric repeats. Variants shifted beyond 2500bp won't be corrected.
3. Homopolymer detection is conservative — only triggers for homopolymer REF ≥3bp with a clear D(n)+SNV pattern. Other miscollapsed events are not detected.

Related¶

• gbcms normalize — CLI command for standalone normalization
• Input Formats — VCF/MAF coordinate conventions
• Allele Classification — How reads are classified after normalization