Allele Classification¶

What call each read shape gets, and the decision behind it, is specified in Read Judgment; this page describes how the engine computes it.

How gbcms classifies each read as supporting the reference allele, the alternate allele, or neither.

Dispatch¶

After passing read filters, each read is dispatched to a type-specific allele checker based on the variant's allele lengths (ref_len × alt_len), not on the variant_type string. This makes dispatch robust against inconsistent upstream type labels:

flowchart LR
    Fetch(["📥 Fetch Reads"]):::fetch --> Filter(["🔍 Apply Filters"]):::filter
    Filter --> Dispatch{{"Dispatch\nref_len × alt_len?"}}

    Dispatch -->|"1 × 1"| SNP["check_snp"]
    Dispatch -->|"N × N"| MNP["check_mnp"]
    Dispatch -->|"1 × N"| INS["check_insertion"]
    Dispatch -->|"N × 1, alt[0]==ref[0]"| DEL["check_deletion"]
    Dispatch -->|"N × 1, alt[0]≠ref[0]\nor complex"| CPX["check_complex"]

    MNP -.->|"inconclusive"| CPX
    INS -.->|"Phase 3 fallback"| CPX
    DEL -.->|"del_len ≥ 5, S3 fail"| CPX

    SNP --> Count(["📊 Update Counts"]):::count
    INS --> Count
    DEL --> Count
    MNP --> Count
    CPX --> Count

    classDef fetch fill:#4a90d9,color:#fff,stroke:#2c6fad,stroke-width:2px;
    classDef filter fill:#e67e22,color:#fff,stroke:#bf6516,stroke-width:2px;
    classDef count fill:#27ae60,color:#fff,stroke:#1e8449,stroke-width:2px;

SNP (Single Nucleotide Polymorphism)¶

A single base substitution — the simplest and most common variant type.

Algorithm¶

flowchart TD
    Start(["🧬 SNP Check"]):::start --> Walk["Walk CIGAR to variant pos"]
    Walk --> CigarOp{"CIGAR op at pos\nis M or X?"}
    CigarOp -->|"No (I/D/N/S)"| Neither0(["Neither — no consumable base"]):::neither
    CigarOp -->|"Yes"| Found{"Position found?"}
    Found -->|No| Neither1(["Neither"]):::neither
    Found -->|Yes| BQ{"Base quality ≥ min_baseq?"}
    BQ -->|No| Neither2(["Neither"]):::neither
    BQ -->|Yes| NCheck{"Base is N?"}
    NCheck -->|Yes| NeitherN(["Neither (has_n_base=true)\n→ n_count++"]):::nbase
    NCheck -->|No| Compare["Compare base to REF and ALT"]
    Compare --> IsRef{"base == REF?"}
    IsRef -->|Yes| Ref(["✅ REF"]):::ref
    IsRef -->|No| IsAlt{"base == ALT?"}
    IsAlt -->|Yes| Alt(["🔴 ALT"]):::alt
    IsAlt -->|No| Neither3(["Neither"]):::neither

    classDef start fill:#9b59b6,color:#fff,stroke:#7d3c98,stroke-width:2px;
    classDef ref fill:#27ae60,color:#fff,stroke:#1e8449,stroke-width:2px;
    classDef alt fill:#e74c3c,color:#fff,stroke:#c0392b,stroke-width:2px;
    classDef neither fill:#95a5a6,color:#fff,stroke:#7f8c8d,stroke-width:2px;
    classDef nbase fill:#e67e22,color:#fff,stroke:#d35400,stroke-width:2px;

Visual Example¶

Variant: chr1:100 C→T (SNP)

Reference: 5'─ ...G  A  T  C  G  A  T  C  G  A... ─3'
                          98 99 100 101
                               ▲
                           variant pos

Read 1:    5'─ ...G  A  T [T] G  A  T  C... ─3'   → ALT ✅
                               ↑
                          base=T matches ALT

Read 2:    5'─ ...G  A  T [C] G  A  T  C... ─3'   → REF ✅
                               ↑
                          base=C matches REF

Read 3:    5'─ ...G  A  T [A] G  A  T  C... ─3'   → Neither
                               ↑
                          base=A ≠ REF or ALT

Insertion¶

Bases inserted after an anchor position. The anchor is the last reference base before the inserted sequence.

Algorithm¶

The insertion check uses a single CIGAR walk with four detection strategies:

Part A — CIGAR-Based Strict Detection¶

flowchart LR
    Start(["🧬 Insertion Check"]):::start --> Walk["Walk CIGAR → left to right"]
    Walk --> MatchBlk{"Match block contains anchor?"}
    MatchBlk -->|No| WinCheck(["→ Windowed Scan"]):::next
    MatchBlk -->|Yes| Bwd{"Anchor at block start\nAND prev op = Ins?"}
    Bwd -->|Yes| BwdMatch{"Same haplotype?\n(quality-masked)"}
    BwdMatch -->|Yes| BWAlt(["🔴 ALT — backward"]):::alt
    BwdMatch -->|No| Fwd
    Bwd -->|No| Fwd{"Anchor at block end?"}
    Fwd -->|No| RefCov["Mark ref coverage"]
    Fwd -->|Yes| Next{"Next op = Ins?"}
    Next -->|No| RefCov
    Next -->|Yes| LenQ{"Length matches?"}
    LenQ -->|Yes| SeqQ{"Seq matches?\n(quality-masked)"}
    SeqQ -->|Yes| StrictAlt(["🔴 ALT — strict"]):::alt
    SeqQ -->|"Confident mismatch"| Third(["⚪ Neither + partial\n(third allele)"]):::partial
    SeqQ -->|"No base readable\n(N or < min_baseq)\nor insert past read end"| FlagURA["Flag has_unreadable_insert"]:::partialflag
    LenQ -->|No| Trunc{"Truncation of expected insert?\n(≥4bp, ≥90% identity,\nboth non-low-complexity)"}
    Trunc -->|Yes| TruncAlt(["🔴 ALT — truncated same event"]):::alt
    Trunc -->|No| WLPartial(["⚪ Neither + partial\n(distinct allele)"]):::partial
    FlagURA --> RefCov
    RefCov --> WinCheck

    classDef start fill:#9b59b6,color:#fff,stroke:#7d3c98,stroke-width:2px;
    classDef next fill:#3498db,color:#fff,stroke:#2471a3,stroke-width:2px;
    classDef alt fill:#e74c3c,color:#fff,stroke:#c0392b,stroke-width:2px;
    classDef partial fill:#8e44ad,color:#fff,stroke:#6c3483,stroke-width:2px;
    classDef fallback fill:#f39c12,color:#fff,stroke:#d68910,stroke-width:2px;

Part B — Windowed Scan + S1/S2/S3 Safeguards¶

flowchart TD
    WinIn(["→ from CIGAR walk"]):::entry --> CheckWin{"Ins within window?\n(±max(5, repeat_span+2),\nor in the shift region; or a readable\nsame-length insert with a placement,\nbase by base or scored, inside the\ndiscrimination window: RJ-21/RJ-22)"}
    CheckWin -->|No| Continue["Continue CIGAR walk"]
    CheckWin -->|Yes| SameLen{"Same length?"}
    SameLen -->|Yes| Readable{"Any of its own\ninserted bases readable?\n(RJ-20)"}
    Readable -->|"No, inside the\ndiscrimination window"| FlagUR["Flag has_unreadable_insert"]:::partialflag
    Readable -->|"No, outside it,\nanother indel in the window"| FlagWL
    Readable -->|"No, outside it:\nseparate event"| Continue
    FlagUR --> Continue
    Readable -->|Yes| Places["Every placement it can take\nas well (slide: the base placed\nonto the reference must fit it)"]
    Places --> AnyAlt{"A placement shows the ALT\nby read bases? (at the junction\nas the strict path; or the variant's\nhaplotype, the only change)"}
    AnyAlt -->|Yes| S2["S2: Track closest match"]
    S2 --> Continue
    AnyAlt -->|No| Change{"The variant written with\nanother change, or another\nindel in the window?"}
    Change -->|Yes| FlagWL
    Change -->|No| Spells{"The read's bases spell\nthe ALT across the window?"}
    Spells -->|Yes| S2
    Spells -->|No| InDW3{"A placement inside the\ndiscrimination window?"}
    InDW3 -->|Yes| FlagAA["Flag another_allele_in_window"]:::partialflag
    FlagAA --> Continue
    InDW3 -->|No| Scored{"Every junction scored by its\nreadable mismatches: a fewest-\nmismatch placement inside? (RJ-22)"}
    Scored -->|Yes| FlagAA
    Scored -->|"No: separate event"| Continue
    SameLen -->|No| FlagWL["Flag has_distinct_allele_nearby\n(any size)"]:::partialflag
    FlagWL --> Continue
    Continue --> MoreOps{"More CIGAR ops?"}
    MoreOps -->|Yes| CheckWin
    MoreOps -->|No| Eval{"Windowed candidate found?"}
    Eval -->|Yes| WinAlt(["🔴 ALT — windowed"]):::alt
    Eval -->|No| URCheck{"has_unreadable_insert?"}
    URCheck -->|Yes| URPartial(["⚪ Neither + partial\n(length, not sequence)"]):::partial
    URCheck -->|No| AACheck{"another_allele_in_window?"}
    AACheck -->|Yes| AAPartial(["⚪ Neither + partial\n(another allele, RJ-21/RJ-22)"]):::partial
    AACheck -->|No| WLCheck{"has_distinct_allele_nearby\nAND ref coverage?"}
    WLCheck -->|"Yes, the read changes\nthe discrimination window"| WLPartial(["⚪ Neither + partial\n(another allele, RJ-7)"]):::partial
    WLCheck -->|"Yes, outside the window"| RefSep(["✅ REF\n(a separate event, RJ-8)"]):::ref
    WLCheck -->|No| HasRef{"Anchor covered by M?"}
    HasRef -->|Yes| Ref(["✅ REF"]):::ref
    HasRef -->|No| AnchorSpan{"Read spans anchor?\n(e.g. soft-clip at anchor)"}
    AnchorSpan -->|Yes| CPXA(["🔄 Phase 3"]):::fallback
    AnchorSpan -->|No| Neither(["⧯ Neither"]):::neither

    classDef entry fill:#3498db,color:#fff,stroke:#2471a3,stroke-width:2px;
    classDef ref fill:#27ae60,color:#fff,stroke:#1e8449,stroke-width:2px;
    classDef alt fill:#e74c3c,color:#fff,stroke:#c0392b,stroke-width:2px;
    classDef neither fill:#95a5a6,color:#fff,stroke:#7f8c8d,stroke-width:2px;
    classDef partial fill:#8e44ad,color:#fff,stroke:#6c3483,stroke-width:2px;
    classDef partialflag fill:#8e44ad,color:#fff,stroke:#6c3483,stroke-width:2px;
    classDef fallback fill:#f39c12,color:#fff,stroke:#d68910,stroke-width:2px;

Windowed Scan Safeguards¶

Three layers of validation prevent false-positive windowed matches:

Visual Example¶

Variant: chr1:100 G→GTG (insertion of TG after anchor G, inside a TG repeat)

Reference:     5'─ ...A  T  G  T  G  T  G  C... ─3'
                          99 100 101 102 103 104 105
                              ▲
                          anchor pos (repeat context — shifted placements stay equivalent)

Read 1 (ALT, strict):  CIGAR = 5M 2I 5M
               5'─ ...A  T  G [T  G] T  G  T  G  C... ─3'
                               └──┘
                          inserted bases at anchor → ALT ✅ (strict)

Read 2 (ALT, windowed): CIGAR = 7M 2I 3M
               5'─ ...A  T  G  T  G [T  G] T  G  C... ─3'
                                         └──┘
                    insertion shifted +2bp (one repeat unit): the same haplotype
                    (TG + TG = TG + TG) → ALT ✅ (windowed); at +1bp it reads GT, a
                    rotation, and is the same haplotype too

Read 3 (REF):  CIGAR = 10M
               5'─ ...A  T  G  T  G  T  G  C... ─3'
                              ↑
                          no insertion after anchor → REF ✅

Deletion¶

Bases deleted after an anchor position. Mirrors insertion but looks for Del CIGAR operations.

Deletions of 50 bp or more¶

A prepared pure deletion of 50 bp or more is judged by its bases at the POS breakpoint, however the aligner wrote it: a D op, a soft clip, or a split read (RJ-24, since 6.7.0). Short-read aligners without realignment write most of these carriers as clips (on HG002, BWA-MEM wrote 1,958 of the junction carriers that way), so a CIGAR-only rule read VAF 0.000 at 287 of 300 heterozygous deletions.

• One window per allele, of equal length, read from the left flank: 4 clean flank bases ending at the anchor, the microhomology (the bases the deletion can slide right through, shared by both alleles), then 4 clean bases of the allele's own continuation: the deleted bases for REF, the right flank for ALT. A REF and an ALT molecule hold their windows from the same read starts.
• Nothing past the window is read, for either allele: checking REF through every deleted base it covers would expose it to errors the ALT reading never meets (measured: HG002 103 vs 100 of 166 deletions in band, the panel 40 vs 29 of 60 within 10% of sign-out).
• Placement: the read is read from its aligned first window base. One whose left flank is not aligned there (in its soft clip, under its own deletion, or written over the deleted bases' end as mismatches) is placed from the right flank, which both alleles share, as the ALT's window sits there; REF's window ends inside the deleted bases, which no ALT molecule holds. The placement is a hypothesis, so the read's confident bases (at or above --min-baseq, not N) must favour the left flank over the deleted bases REF has before the right flank by one base read at --min-baseq, and by more than nothing (a guard-masked flank base weighs nothing for ALT, the deleted base it is weighed against stays read): a REF molecule starting inside the deleted span behind a poor clip shows the same right flank as a carrier, and nothing of the left flank. It is ALT where it fits, else another allele; written over the deleted bases' end, only where it fits. Hard-clipped bases are absent, so a hard-clipped record holds no window on its own.
• Outcomes: REF or ALT when the read fits exactly one window (bases below --min-baseq, and N, fit anything; ALT also needs one minimum-quality base's evidence); fitting neither, another allele (partial_alt when closer to ALT); not holding the windows, or crossing a splice inside them, depth only. A read fitting REF whose own alignment deletes or inserts 5+ bases inside the deleted span is another allele (C26).
• Depth: a decided read whose windows read soft-clipped bases inside its fragment counts toward depth in DNA, as for complex variants.
• A change recurring at one window offset on an allele's reads (3+ fragments, 20% of that allele's reads) is a separate event, masked at that offset in both windows: in the flank and microhomology, which both alleles share, one the REF reads carry too (a germline SNP or MNP, a run included, on 3 or more REF reads); one of the ALT reads' own with no own change at either neighbouring base (the base either side of the windows included); and past the junction, one on the REF reads with none recurring beside it (a run there, on reads fitting REF within a base, is another breakpoint: another allele). A run of changes on the ALT reads alone is an indel at the junction: another allele, not masked. The isolated changes on the ALT reads alone name the larger allele in gbcms_diagnostic (the deletion with the changes); beside a run none is named. Masks never take every offset where REF and ALT differ: those stay read, and a read changed there is neither.
• Another unit count beside a run (#297, RJ-24 amended 2026-10-30 and 2026-10-31, rule X): where the bases past the deletion repeat with the period of another allele at its breakpoint (one more or one fewer unit of a tandem run whose deleted copy differs from it, so the windows stop at the copy's differing base), that allele reads like the ALT through the windows and on to its first difference. Every such allele is computed exactly from the reference (the other deletions at the deletion's breakpoint, fewer bases or more up to twice the given, and the listed rows; a reference N tells nothing), never by a run detector. An insertion at the breakpoint (a tandem copy of the bases before it) is not looked for: a REF read inside the copy fits it as well as REF. A listed one is always a candidate (rule X, operator 2026-10-31: the input is taken as correct, so no read need show the count, and a listed count is held wherever it first differs, past every read's reach too: the reads that cannot reach the difference are partial evidence); one nobody listed counts only where it is present at the site, spelled in full by its own bases on 3 or more fragments the aligner places here (XS below AS; MAPQ above 0 where it writes no XS) and 5% of those telling the ALT from it: read on 24 bases past its first difference, fitting it but for one base past the windows, with its base at two or more of the bases where it and the ALT differ and the ALT's at none, and fitting REF worse by two) and the longest read at the site can reach that difference, an ALT call stands only where the read's confident bases where the two differ favour the ALT by one base read at --min-baseq (a base fitting neither costs both alike, an error; one base of the other's against three of the ALT's is forgiven); otherwise the read is the ALT or that allele: neither, with partial evidence (RJ-25), and a read whose own bases spell an unlisted count in full is never ALT whatever the presence floor (partial evidence; depth only where the aligner doubts its placement), the floor deciding only the reads that cannot tell. So the given count's reads that end before the first difference stay partial evidence at a listed count that no molecule carries (the cost of rule X: the S1 rows in the listed form); the same reads with the count unlisted and unshown are ALT. A read that fits the deletion with a listed row applied past the junction (an indel or a substitution the input names) as well, but for one base and whose confident bases show that row's change, is the deletion with that neighbour: it spells no other allele and its ALT call stands; a row no read carries makes no read a neighbour's. The windows do not grow: REF, told from both at the junction, keeps its reads and the sibling claim its window. A site with no other count keeps its reads. Where an unlisted count's difference lies past what a read can spell (the first difference and 24 bases), nothing is held, the counts stand and gbcms_diagnostic names the row, UNIT_COUNT_UNSPANNED(change:bases/longest). Under the sequence judge the alleles found so are site states, so its reads tie the same way. The allele that holds the calls is named, OBSERVED_ALLELE or COEXISTING_ALLELE at the deletion's anchor, listed or not.
• RNA: the same, after splice triage (a deletion written as a splice N stays excluded).

The length band that used to credit a deletion of another length within 3 bases (and split D+M+D representations) is retired: it credited no real read on the panel or HG002. A deletion the prepared reference cannot hold the windows of, or an unprepared one, keeps the CIGAR checks below.

Algorithm (CIGAR checks)¶

Same single-walk strategy as insertion, with three additional features:

1. Wrong-length rule — A wrong-length D at the anchor is a distinct allele in the same tract → neither + partial_alt, never Phase 3. See Wrong-Length Pure Indels.

2. A windowed deletion that fails S3 — an exact-length deletion elsewhere is ALT when the read's bases spell the ALT despite where its gap sits (compensating mismatches); otherwise, at 5bp or more, another allele (has_distinct_allele_nearby). Short deletions (<5bp) that fail S3 remain CIGAR-definitive REF — they are almost certainly unrelated spurious deletions. See Case 3 in the Complex Indels guide.

3. Haplotype fallback — A read that spans the anchor with no aligned base on it and no deletion over it (in practice, its anchor inside a splice N) goes to Phase 3 for haplotype-based comparison. Reads mapping entirely inside a large deleted span carry no information about the variant and are classified neither.

Part A — CIGAR-Based Strict Detection¶

flowchart LR
    Start(["🧬 Deletion Check"]):::start --> Walk["Walk CIGAR → left to right"]
    Walk --> MatchBlk{"Match block contains anchor?"}
    MatchBlk -->|No| WinCheck(["→ Windowed Scan"]):::next
    MatchBlk -->|Yes| AtEnd{"Anchor at block end?"}
    AtEnd -->|No| RefCov["Mark ref coverage"]
    AtEnd -->|Yes| NextDel{"Next op = Del?"}
    NextDel -->|No| RefCov
    NextDel -->|Yes| LenMatch{"Length matches exactly?"}
    LenMatch -->|Yes| StrictAlt(["🔴 ALT — strict"]):::alt
    LenMatch -->|No| WLPartial(["⚪ Neither + partial\n(distinct allele)"]):::partial
    RefCov --> WinCheck

    classDef start fill:#9b59b6,color:#fff,stroke:#7d3c98,stroke-width:2px;
    classDef next fill:#3498db,color:#fff,stroke:#2471a3,stroke-width:2px;
    classDef alt fill:#e74c3c,color:#fff,stroke:#c0392b,stroke-width:2px;
    classDef partial fill:#8e44ad,color:#fff,stroke:#6c3483,stroke-width:2px;

Part B — Windowed Scan + S1/S2/S3 Safeguards¶

flowchart TD
    WinIn(["→ from CIGAR walk"]):::entry --> CheckWin{"Del within window?\n(±max(5, repeat_span+2),\nor a placement in the shift region)"}
    CheckWin -->|No| Continue["Continue walk"]
    CheckWin -->|Yes| S1{"S1: Length check"}
    S1 -->|"Exact match"| S3{"S3: Same haplotype?"}
    S1 -->|"Wrong length, ≥5bp"| FlagWL["Flag has_distinct_allele_nearby"]:::partialflag
    S1 -->|"Wrong length, <5bp"| Continue
    S3 -->|Yes| Only{"The read's only change\nacross the window?"}
    Only -->|Yes| S2["S2: Track closest"]
    Only -->|No| FlagWL
    S3 -->|No| DSpells{"Exact length and the read's\nbases spell the ALT?"}
    DSpells -->|Yes| S2
    DSpells -->|"No, exact length ≥ 5"| FlagWL
    DSpells -->|"No, < 5bp"| Continue
    S2 --> Continue
    FlagWL --> Continue
    Continue --> MoreOps{"More ops?"}
    MoreOps -->|Yes| CheckWin
    MoreOps -->|No| Eval(["→ Evaluation"]):::next

    classDef entry fill:#3498db,color:#fff,stroke:#2471a3,stroke-width:2px;
    classDef next fill:#3498db,color:#fff,stroke:#2471a3,stroke-width:2px;
    classDef partialflag fill:#8e44ad,color:#fff,stroke:#6c3483,stroke-width:2px;

Part C — Evaluation¶

flowchart TD
    EvalIn(["→ from windowed scan"]):::entry --> Win{"Windowed match found?"}
    Win -->|Yes| WinAlt(["🔴 ALT — windowed"]):::alt
    Win -->|No| Spans{"Ref coverage at anchor?"}
    Spans -->|No| AnchorSpan{"Read spans anchor?\n(its anchor in a splice N)"}
    AnchorSpan -->|Yes| CPX(["🔄 Phase 3"]):::fallback
    AnchorSpan -->|No| Neither(["⧯ Neither\n(no variant information)"]):::neither
    Spans -->|Yes| WLCheck{"has_distinct_allele_nearby?"}
    WLCheck -->|"Yes, the read changes\nthe discrimination window"| WLPartial(["⚪ Neither + partial\n(another allele, RJ-7)"]):::partial
    WLCheck -->|"Yes, outside the window"| RefSep(["✅ REF\n(a separate event, RJ-8)"]):::ref
    WLCheck -->|No| Ref(["✅ REF"]):::ref

    classDef entry fill:#3498db,color:#fff,stroke:#2471a3,stroke-width:2px;
    classDef ref fill:#27ae60,color:#fff,stroke:#1e8449,stroke-width:2px;
    classDef alt fill:#e74c3c,color:#fff,stroke:#c0392b,stroke-width:2px;
    classDef neither fill:#95a5a6,color:#fff,stroke:#7f8c8d,stroke-width:2px;
    classDef partial fill:#8e44ad,color:#fff,stroke:#6c3483,stroke-width:2px;
    classDef fallback fill:#f39c12,color:#fff,stroke:#d68910,stroke-width:2px;

Windowed Scan Safeguards¶

Three layers of validation prevent false-positive windowed matches:

Visual Example¶

Variant: chr1:100 GTG→G (deletion of TG after anchor G, inside a TG repeat)

Reference:     5'─ ...A  T  G  T  G  T  G  C... ─3'
                          99 100 101 102 103 104 105
                              ▲
                          anchor pos (repeat context — shifted placements stay equivalent)

Read 1 (ALT, strict):  CIGAR = 5M 2D 5M
               5'─ ...A  T  G  ──  ──  T  G  C... ─3'
                               └─────┘
                          2bp deletion at anchor → ALT ✅ (strict)

Read 2 (ALT, windowed): CIGAR = 7M 2D 3M
               5'─ ...A  T  G  T  G  ──  ──  C... ─3'
                                        └─────┘
                  deletion shifted +2bp (one repeat unit): S1 length matches,
                  S3 ref bases at shifted position (TG) match expected TG → ALT ✅ (windowed)

Read 3 (REF):  CIGAR = 12M
               5'─ ...A  T  G  T  G  T  G  C... ─3'
                              ↑
                          no deletion after anchor → REF ✅

Wrong-Length Pure Indels (partial_alt)¶

A read whose CIGAR proves a pure indel of a different length at the variant anchor is evidence of a different allele — never definitive REF, and never the queried ALT. This is the load-bearing rule for repeat-tract loci, where coexisting distinct-length populations are distinct slippage alleles (a 1bp germline slippage allele coexisting with an annotated 2bp somatic deletion; a trinucleotide-repeat deletion ladder). Counting every tract-touching indel as the annotated event inflated VAF several-fold at such loci.

The rule¶

Phase 3 deliberately does not arbitrate the definitive wrong-length/wrong-sequence cases: its haplotype window is length-blind inside repeat tracts, and alignment scoring promotes a wrong-sequence insert to ALT because it still beats the gapped REF alignment. Delins/complex variants are unaffected — they route to check_complex, whose Phase-3 realignment correctly resolves split and mismatch-absorbed representations of one event.

Deletions of 50 bp or more¶

A deletion of another length is another allele at any size. From 50 bp a prepared deletion is judged by its bases at the POS breakpoint (RJ-24); the 3-base length band that used to credit another length there, on the premise that a large D op is never an artifact, is retired: it credited no real read on the panel or HG002.

Reading the output¶

At an affected locus, alt_count reflects only exact-support reads, and the distinct-allele evidence appears in partial_alt/any_alt. When partial evidence dominates (partial_alt > alt_count), the variant is flagged PARTIAL_DOMINANT in gbcms_diagnostic — at pure-indel loci this usually means the locus carries a coexisting allele the annotation does not describe (at ≥50bp deletion loci, a real different large event). Per-read indel lengths are visible with --trace.

Splice-Aware Evidence (RNA)¶

Spliced reads (CIGAR N) pass through a triage before any of the per-type checkers: a read testifies only through aligned bases (or a D op) at the discriminating positions, and an N spanning all of them means the read observes nothing there — it classifies neither AND is excluded from DP/fragment depth (samtools-pileup semantics). D is deletion evidence; N never is, so the call cannot flip on the aligner's D-vs-N representation choice. Indel ops directly after a splice N (M-N-D-M) get the same anchor/windowed inspection as ops after an M block, and Phase 3 never scores across a splice (window extraction and check_complex reconstruction both refuse N-crossing windows rather than stitch exon arms into a junction-chimeric sequence). At a pure-deletion locus whose anchor is spliced out, a junction read covering the entire deleted span with aligned bases counts REF (span-aligned REF testimony — the deletion is demonstrably absent); partial coverage or competing indel evidence keeps the read on its existing path. Details: RNA Splice-Junction Handling.

MNP (Multi-Nucleotide Polymorphism)¶

Multiple adjacent bases substituted simultaneously.

Algorithm¶

The MNP checker uses masked per-position evaluation — each discriminating position is independently assessed, and only unmasked positions vote for the final classification.

flowchart TD
    Start(["🧬 MNP Check"]):::start --> Find["Find read position of first base"]
    Find --> Found{"Position found?"}
    Found -->|No| CPX
    Found -->|Yes| Cover{"Read covers entire MNP region?"}
    Cover -->|No| CPX
    Cover -->|Yes| Contig{"No indels within MNP block?"}
    Contig -->|"Indel found"| CPX
    Contig -->|"Contiguous"| PerPos["Per-position masked evaluation"]

    PerPos --> Loop{"For each discriminating pos\n(where REF ≠ ALT)"}
    Loop --> Mask{"BQ < min_baseq\nOR base is N?"}
    Mask -->|Yes| Masked["Position masked\n(cannot vote)"]
    Mask -->|No| Vote{"Unmasked: compare\nbase to REF/ALT"}
    Vote --> AltVote["ALT vote"]
    Vote --> RefVote["REF vote"]
    Vote --> OtherVote["Other vote"]

    Masked --> NextPos( )
    AltVote --> NextPos
    RefVote --> NextPos
    OtherVote --> NextPos
    NextPos --> MorePos{"More positions?"}
    MorePos -->|Yes| Loop
    MorePos -->|No| Classify

    Classify{"All unmasked = 0?"}
    Classify -->|Yes| LowQ(["⬜ LowQuality\nNeither + partial tracked"]):::neither
    Classify -->|No| AllAlt{"All unmasked\nmatch ALT?"}
    AllAlt -->|Yes| Alt(["🔴 ALT"]):::alt
    AllAlt -->|No| AllRef{"All unmasked\nmatch REF?"}
    AllRef -->|Yes| Ref(["✅ REF"]):::ref
    AllRef -->|No| Third(["⬜ ThirdAllele\n(partial_alt tracked)"]):::neither

    CPX(["🔄 check_complex\n(Phase 3 fallback)"]):::fallback

    classDef start fill:#9b59b6,color:#fff,stroke:#7d3c98,stroke-width:2px;
    classDef ref fill:#27ae60,color:#fff,stroke:#1e8449,stroke-width:2px;
    classDef alt fill:#e74c3c,color:#fff,stroke:#c0392b,stroke-width:2px;
    classDef neither fill:#95a5a6,color:#fff,stroke:#7f8c8d,stroke-width:2px;
    classDef fallback fill:#f39c12,color:#fff,stroke:#d68910,stroke-width:2px;

Classification Outcomes¶

Contiguity Check¶

The contiguity check is performed first (before quality or sequence comparison) as a fail-fast for structural issues. gbcms compares the read positions of the first and last MNP base — if the distance doesn't equal len - 1, an indel exists within the block and the read is routed to check_complex for haplotype-based resolution.

Complex (Indel + Substitution)¶

Variants where REF and ALT differ in both sequence and length. It also covers complex Del+SNV variants (N×1 with a substituted anchor), complex Ins+SNV variants (1×N with a substituted anchor, since 6.6.0) and MNP reads that carry an indel inside the block.

The exact-carrier rule¶

A read is REF or ALT for a complex variant only when its own bases carry that allele across the whole event. Nothing is inferred: gbcms counts the allele it is given.

• The windows. Each allele's window is its event plus two reference bases of flank on each side.
• The event is every base where REF and ALT can differ over all equivalent placements: the union of the minimal difference trimmed from the left first and from the right first, grown through any tandem repeat touching either end on either allele (unit 1–6 bases, or the length change). An aligner may place an indel anywhere in such a repeat, including a run the ALT's own bases continue: for TA>ACC after AAA, the inserted A can sit anywhere in the run. So a delins in or beside a repeat is judged over its whole ambiguity.
• One pair of equal-length windows at one flank (RJ-23, since 6.7.0). Both alleles are read inward from the same flank: the shorter allele's whole window, and as many of the longer allele's bases. A read holding them is read on through the rest of the longer allele's window as far as it reaches, and a later base that contradicts an allele rules that allele out (one A more than the ALT's run is not the ALT). Windows of equal length at one flank give a REF and an ALT molecule the same number of read starts that can judge them, so neither allele is favoured by where reads start. The flank is the POS (left) flank; an RNA variant whose left flank an exon edge cuts reads from the right, and a read judged there counts although it starts past POS (inside the longer allele), as the other allele's molecule starting at the same place would: otherwise the POS gate would drop only the longer allele's reads. Before 6.7.0 the shorter window was padded with flank to the longer one's length, which cost reads (a third of them on long delins), and long events were read at both junctions, which let the longer allele count at either end. The windows run on past the event through any bases the shorter allele's flank shares with the longer allele's next bases, and on through two bases from where they first differ, as a long deletion's run through its microhomology: a read carrying only a substitution of the event's bases then shows the longer allele's bases where the shorter allele's differ, and is neither allele.
• The window is read where it sits. It is read from the read position aligned to the window's first reference base, and back from the one aligned to its last; the closer reading counts. The flank bases are outside the event, so they align the same way on either allele. Inside the event the read's bases are taken as one stretch, so where the aligner put an indel or a soft clip does not matter, and an indel it placed past one flank leaves the other in place. A copy of the window elsewhere in the read is never a match.
• Placed from the other flank. A read whose reading flank lies in its soft clip, under its own deletion, or aligned at the event's other copy (an aligner placing a carrier by its far part, its flank over the event's last bases as mismatches) is placed from the other flank as the allele its bases favour would sit there. Each allele's window is placed from that flank at its own length, and the read's confident bases (at or above --min-baseq, not N) must mismatch one allele by at least one base read at --min-baseq more than the other; bases below it weigh nothing, so a REF molecule behind a poor clip, whose clip happens to fit the ALT, shows no placement and is depth only. The read counts for that allele where it fits; otherwise it carries another allele (partial_alt when the favoured allele is the ALT). One placed through aligned bases counts only where its allele fits, and counts although it starts past POS, as the other allele's molecule starting there would.
• Masked bases. Bases below --min-baseq (and N) match anything. This is the only tolerance, and it is the same quality rule every backend applies. n_count counts a read only when an N sits among the compared bases.
• Outcome. A read decides only when it holds both the REF and the ALT window.

A read one confident base off the given ALT carries a different allele. A read that ends before holding the pair cannot show either, like a pure-indel read ending inside its tract; a read that holds it and ends inside the longer allele is judged by the bases it has. - Long events. When both alleles' windows exceed 50 bases, no read can hold the shorter one whole. The pair then ends one base past the first base where the alleles differ, read from the flank nearer it (an event ending a long run is out of a read's reach from the other), and is read on as above. Read from the right flank, a read the rule decides counts although it starts past POS, as for an RNA exon edge below. - A change beside the event. A change that recurs within 25 bases of the event, at a fixed distance, on the fragments carrying either allele's core (at least 3 fragments, overlapping mates once, and 20% of those reading the base) is a separate event, as another indel outside the window is (RJ-8), masked in both windows so a read carrying the given allele counts whatever it shows there: a change the REF fragments carry too (a germline SNP or MNP, background on both haplotypes, a run included; on 3 or more REF fragments, as one REF read's error is no background), and an isolated change of the ALT fragments' own (no own change at either neighbouring base, the base past the window included). A run of changes on the ALT fragments alone is an indel or MNP beside the event: another allele, not masked. The isolated changes on the ALT fragments alone (under 5% of the REF fragments) make their haplotype larger than the given allele (beside a run none is named: the per-base tallies cannot spell an indel; the observed allele names what the reads' alignments carry), a phased neighbour or an event the input under-describes; the counts stay the given allele's, and gbcms_diagnostic names the larger allele under the rule every observed allele follows (OBSERVED_ALLELE or COEXISTING_ALLELE: more carriers than the given allele, 5% of the reads read, not already an input row). - A change inside the grown event (RJ-23, amended 2026-10-14). A change the REF molecules carry inside the repeat the event was grown through, off the given allele's own bases (a germline SNP in the run), is background too: masked in both windows where it sits on each allele (left of the given allele at the same place, right of it shifted by the length change), so a REF read carrying it counts REF and a carrier carrying it counts ALT. It is found from the molecules placed by their bases (below), so a read written with its gap elsewhere is still read. A change on the ALT molecules alone stays unmasked for an event changing the length (it may be the event's own sequence); for an MNP, which has nothing to place, the run is beside it, so an isolated one is its phased neighbour, masked and named. Carriers all clipped short of the windows' reach cannot show such a neighbour; a carrier aligned with a clip or an indel at the block counts once its fellows show it. An MNP's bases it keeps inside its core sit on both alleles as well: background there (a germline SNP, in cis or in trans) is masked too, on both routes; a change there on the ALT molecules alone is another allele. A mask inside the event's span takes the reference's base and the background's only: a third base there is still another allele's (beside the event a mask takes any base). Inside a delins the bases kept depend on the alignment, so a change there is not masked. - Masks that hide the difference (RJ-27, amended 2026-10-20, #294). In a repeat whose units differ at single bases, the base telling a one-unit insertion from REF can be the one the guard masks (a germline SNP turning the next unit into the inserted one, so both haplotypes read alike over a short read). The windows then read on to the next base telling the alleles apart, which prep fetches far enough to hold whatever the guard masks; a read short of it is depth only unless its own bases, held as far as before, fit one allele alone and exactly. Leaving such a mask unapplied was tried and rejected: it counted the SNP's reads as the insertion. - The ALT's own base, shifted. Where the aligner writes an ALT read without its gap, the read shows past the shorter allele's tract the ALT's base there, which looks like a change on the REF fragments; a germline SNP can give the same base. Where it is the REF fragments' most common change there, background or on more than 1% of them, the fragments placed apart from the event tell them apart: those reading the base whose bases where the reference and that shifted ALT differ past it (aligned with it, no length change between) show the reference's more often (within 300 bases of the variant's REF span). Carried by three or more of them, and 3% of them (above the errors at one base, random or most systematic ones), it is a haplotype's change and is masked, though it be under the REF fragments' background share (a minority haplotype, after a tumour's gain or loss, reads as the ALT over a short read). Otherwise, read by 15 or more, a background change there is the ALT's own base: not masked, and spelling no other allele; with fewer it is masked where it is background. Where the reads placed apart come in both pair orientations, carriers too one-sided against that balance for a haplotype's (a binomial chance under 0.1%; twelve all one way, not eight) are an artifact (oxoG, FFPE), not a haplotype's change. A haplotype showing on fewer of the reads placed apart is taken for errors, and its reads count ALT: with 30 to 45 of them, an SNP haplotype on 5–10% of the molecules in about a third of runs; at depth, any under 3% of the molecules. Such a row, carried by two or more of the reads placed apart and over 1% of them, is named in gbcms_diagnostic (POSSIBLE_MINORITY_HAPLOTYPE). A systematic error giving that base on 3% of them in both orientations masks it, and the ALT's reads short of the next telling base are depth only. - An error both readings share decides nothing. A read whose bases place no length change and whose REF and ALT readings miss exactly the same confident bases is read on past the windows with those bases unread, not counted for the allele whose tolerance happens to fit the length the aligner wrote. Past the windows one error there the other allele's reading does not share is forgiven where that reading misses three or more, the read writes that allele's length, and the error recurs on no other fragment (recurring, it is another allele's base); within the windows only a shared error is tolerated. - Reads placed a unit or a copy off are read whole. A read whose alignment is not clean (another event of its own inside the discrimination window, or one its bases place nowhere; or, placing no length change, two or more confident changes from the reference within 300 bases of the variant's REF span) and that is not spliced is placed ungapped on each allele's haplotype, every readable base of it, over that reach of reference, the guard's masks as in the windows (other alleles at the locus are built on the reference before them). An allele fits where its windows read exactly and at most one confident base elsewhere mismatches. The read is that allele where it fits so and the other allele (its event's bases open, weighed at every placement holding as much of it as the shorter allele's windows, from either flank) and every other allele at the locus miss by more. A read fitting the longer allele only from inside its event is depth only (RJ-23); otherwise the read keeps its call. In an imperfect repeat a read that reaches no unique sequence either side is still placed by the aligner (#290). - Other alleles at the locus (RJ-25, since 6.7.0). A read that fits the ALT counts ALT only when, for every other allele of the ALT's length at the locus, it reads a base telling the two apart (the ALT's base, at or above --min-baseq); otherwise it is that allele or the ALT, neither with partial evidence. Where the reads show no other allele nothing changes: a carrier reading only the event's length, or an MNP read with a differing base below --min-baseq, counts ALT. The other alleles are every co-annotated sibling of the ALT's length whose span meets the event's site and that differs from it inside the window (a listed allele is always another allele: the input is taken as correct, so no fragment need show it, since 2026-11-02), and those the variant's own reads spell in full on 3 or more fragments and on 5% of the fragments spelling the allele they are nearest (REF or the ALT; the larger when as near), placed by their bases from the left flank as the windows read them, not where the aligner wrote the gap. An SNV row where the same input lists an indel or delins whose span meets the row's site is counted by the sequence judge also under the rules (the single-base check weighs no listed allele that changes the length, so a read that spells one counted ALT; RJ-25, 2026-11-03). MNV rows stay on the rules. Not another allele: one whose every difference from the ALT is a change from the reference the REF molecules carry too (a separate event on both haplotypes); one differing only where the guard masks (the ALT with a separate event); and, for an MNP, one differing only beside the MNP's own bases (the ALT with a phased neighbour). A molecule is REF's or the ALT's when its bases, the event's included, favour that allele over the other placed the same way, with at most two changes beside the event. The base-by-base MNP check below follows the same rule for every read it calls ALT. A read that fits REF follows it beside the other alleles of REF's length (not another allele: one whose every difference from REF is a change from the ALT's base the ALT molecules carry too, at a base the event keeps), on both routes, as far as the read reaches (since 2026-11-02): an allele whose first to last difference from REF the read's span does not overlap does not contest it, as a read ending inside a long event is judged by the bases it has (RJ-9); one the read reaches and cannot tell (the base unread, below --min-baseq or masked) leaves it neither. The core of a delins that replaces 50 REF bases or more (4 bases in from both ends) is as a long deletion's deleted bases: a substitution wholly inside it is a context of REF's molecules, so a read that carries it keeps its REF call (RJ-24). A read the MNP check calls REF that is ALT for a co-annotated sibling there is partial evidence (RJ-5); on the exact-carrier route such a read does not fit REF, and is neither. The 1-fragment-and-5% test that made a listed component of an event part of it (2026-10-14) is gone for listed alleles (2026-11-02): a carrier that reads only a component's base cannot tell the event from it and is partial evidence. It stays for an allele nobody listed (3 fragments and 5%). - MNP reads with an indel. An MNP read with an insertion or deletion in the block or right beside it, or a clip in it, is judged by this rule; an aligner may write a block shifted by one base as an insertion before it and a deletion after. Its ALT calls with every window base read count toward mnp_confirmed_alt. Other MNP reads, including those with an indel further off (a germline one nearby), are compared base by base. - Carriers in soft-clipped bases (DNA). A read the rule decides from its own bases, with a window that reads its soft-clipped bases, counts even when its aligned span stops short of the variant position, provided every base read lies inside a well-defined fragment (past the fragment end a clip is adapter). Undecided clipped reads stay out of depth. See counting metrics. - Spliced reads (RNA). A read spliced in a window's flank is judged on windows built over the reference spliced at its own junctions (the far exon read from the FASTA): its next bases come from the next exon, so both alleles continue there, and growth and flank are measured on the spliced sequence. A junction starting inside the bases where the alleles differ and running past them (or ending inside them) splices the haplotypes at the event's edge; any other splice through those bases leaves the read depth only. An RNA read's soft clip that reaches an exon edge or junction end is not read (it may hold the next exon's bases); other clips are. - Same on both backends. The rule uses no alignment scoring, so pairhmm and sw give identical counts for complex variants. - When the rule cannot run. prepare_variants fetches reference until it holds the event, grown through any repeat, with its flank and the neighbouring-change reach of 25 bases, and until the windows, with every base the guard could mask masked, reach a base telling the alleles apart (up to 16,384 bases each side, clamped at a contig end: the window holds the bases that exist). For reading reads whole it also holds a whole read's reach, 300 bases each side of the change (Variant.whole_ref, not exposed to Python), which every bin's fetch covers. A variant without that reference (one built without prepare_variants, or an event whose two flank bases would pass a contig end) keeps the previous classifier below. Every variant of a normal run is prepared. A record stored without its bases (SEQ *) is not read.

Previous classifier (variants without a reference; pure-indel Phase 3)¶

The phases below still classify complex variants the exact-carrier rule cannot judge (see above). They are also the Phase-3 fallback that check_insertion and check_deletion use for pure indels. They are a three-phase algorithm with quality-aware matching and Smith-Waterman fallback.

Phase Overview¶

stateDiagram-v2
    direction LR
    [*] --> Phase1 : read enters check_complex

    Phase1 : Phase 1\nHaplotype Reconstruction\n─ CIGAR walk → read_seq
    Phase2 : Phase 2\nMasked Comparison\n─ quality-masked haplotype match
    Phase25 : Phase 2.5\nEdit Distance\n─ Levenshtein fast-path
    Phase3 : Phase 3\nWFA → PairHMM\n─ dual-haplotype alignment

    Phase1 --> Phase2 : always
    Phase2 --> REF : exact REF match
    Phase2 --> ALT : exact ALT match
    Phase2 --> Neither : ambiguous, or partial ALT bases
    Phase2 --> Phase25 : recon matches NEITHER allele length
    Phase2 --> Phase3 : case mismatch (no partial) / skip guards
    Phase25 --> REF : edit dist REF wins (>1 margin)
    Phase25 --> ALT : edit dist ALT wins (>1 margin)
    Phase25 --> Phase3 : tied or guard-skipped
    Phase3 --> REF : confident REF call
    Phase3 --> ALT : confident ALT call
    Phase3 --> Neither : below threshold

    REF : ✅ REF
    ALT : 🔴 ALT
    Neither : ⚪ Neither

Phase 1 + 2 Detail¶

flowchart TD
    Start(["🧬 Complex Check"]):::start --> Walk["Walk CIGAR → rebuild read_seq"]

    Walk --> Ops{"CIGAR op type?"}
    Ops -->|"M / = / X"| Match["Append overlapping bases + quals"]
    Ops -->|"I"| InsOp["Append inserted bases + quals"]
    Ops -->|"D / N"| Del["Advance ref_pos only"]
    Ops -->|"S (soft-clip)"| SC{"Overlaps variant?"}
    SC -->|Yes| SCAdd["Append soft-clipped bases"]
    SC -->|No| SCSkip["Skip"]
    Match --> LoopBack( )
    InsOp --> LoopBack
    Del --> LoopBack
    SCAdd --> LoopBack
    SCSkip --> LoopBack
    LoopBack --> |"more ops?"| Ops
    LoopBack --> |"done"| Guards

    Guards["Post-reconstruction guards"]
    Guards --> LargeGuard{"ref_len > max(50, read_len/3)
AND recon < 10% of ref?"}
    LargeGuard -->|Yes| P3Skip(["→ Phase 3 directly"]):::next
    LargeGuard -->|No| LenRatio{"ref_len > 2 × alt_len AND
recon matches ALT length only?"}
    LenRatio -->|Yes| P3Direct(["→ Phase 3 directly"]):::next
    LenRatio -->|No| LenCheck{"Recon length matches?"}

    LenCheck -->|"ALT and REF"| CaseA["Case A: Dual compare
+ ambiguity detection"]
    LenCheck -->|"ALT only"| CaseB["Case B: ALT-only compare"]
    LenCheck -->|"REF only"| CaseC["Case C: REF-only compare"]
    LenCheck -->|"Neither"| P25(["→ Phase 2.5"]):::next

    CaseA --> Ambig{"Reliable bases
match BOTH?"}
    Ambig -->|Yes| Neither1(["⚪ Neither — ambiguous"]):::neither
    Ambig -->|No| Matches{"Which allele?"}
    Matches -->|ALT| Alt1(["🔴 ALT"]):::alt
    Matches -->|REF| Ref1(["✅ REF"]):::ref
    Matches -->|"Neither, some reliable
ALT-matching bases"| NeitherP(["⚪ Neither + partial"]):::neither
    Matches -->|"Neither, none"| P3F(["→ Phase 3"]):::next

    CaseB --> AltMatch{"0 mismatches on
reliable bases?"}
    AltMatch -->|Yes| Alt2(["🔴 ALT"]):::alt
    AltMatch -->|"No, some reliable
ALT-matching bases"| NeitherP
    AltMatch -->|"No, none"| P3F

    CaseC --> RefMatch{"0 mismatches on
reliable bases?"}
    RefMatch -->|Yes| Ref2(["✅ REF"]):::ref
    RefMatch -->|No| P3F

    classDef start fill:#9b59b6,color:#fff,stroke:#7d3c98,stroke-width:2px;
    classDef next fill:#3498db,color:#fff,stroke:#2471a3,stroke-width:2px;
    classDef ref fill:#27ae60,color:#fff,stroke:#1e8449,stroke-width:2px;
    classDef alt fill:#e74c3c,color:#fff,stroke:#c0392b,stroke-width:2px;
    classDef neither fill:#95a5a6,color:#fff,stroke:#7f8c8d,stroke-width:2px;

Phase 2.5 Detail — Edit Distance Fast-Path¶

flowchart TD
    P25in(["→ from Phase 2
(recon matches NEITHER allele length)"]):::entry
    P25in --> ReconLen{"recon_len ≥ 2?"}
    ReconLen -->|No| P3(["→ Phase 3"]):::next
    ReconLen -->|Yes| SizeGuard{"REF ≤ 50bp AND ALT ≤ 50bp?"}
    SizeGuard -->|No| P3
    SizeGuard -->|Yes| RatioGuard{"ref_len > 2 × alt_len?"}
    RatioGuard -->|Yes| P3
    RatioGuard -->|No| EditDist["Compute Levenshtein distance
to REF and ALT alleles"]
    EditDist --> EDMargin{"> 1 edit margin?"}
    EDMargin -->|"ALT closer"| Alt25(["🔴 ALT"]):::alt
    EDMargin -->|"REF closer"| Ref25(["✅ REF"]):::ref
    EDMargin -->|"Tied"| P3

    classDef entry fill:#9b59b6,color:#fff,stroke:#7d3c98,stroke-width:2px;
    classDef next fill:#3498db,color:#fff,stroke:#2471a3,stroke-width:2px;
    classDef ref fill:#27ae60,color:#fff,stroke:#1e8449,stroke-width:2px;
    classDef alt fill:#e74c3c,color:#fff,stroke:#c0392b,stroke-width:2px;

Phase 3 Detail — WFA → PairHMM (see Alignment Fallback)¶

flowchart LR
    P3in(["→ from Phase 2.5"]):::entry --> Extract["Extract raw read window"]
    Extract --> Backend{"--alignment-backend?"}

    Backend -->|"pairhmm (default)"| WFA["WFA edit-distance triage
vs pangenomic haplotypes"]
    WFA --> WFAClear{"Clear winner?"}
    WFAClear -->|"Yes (~70-80%)"| WFADone{"Which allele?"}
    WFADone -->|ALT| AltW(["🔴 ALT"]):::alt
    WFADone -->|REF| RefW(["✅ REF"]):::ref
    WFAClear -->|"Ambiguous"| HMM["Marginalized PairHMM LLR"]
    HMM --> LLR{"LLR > threshold?"}
    LLR -->|Confident| HMMDone{"Which allele?"}
    LLR -->|"Below threshold"| N4(["⚪ Neither"]):::neither
    HMMDone -->|ALT| AltH(["🔴 ALT"]):::alt
    HMMDone -->|REF| RefH(["✅ REF"]):::ref

    Backend -->|"sw"| SW["Semiglobal SW alignment
vs REF + ALT haplotypes"]
    SW --> Margin{"Score diff ≥ 2?"}
    Margin -->|"Confident"| ConfCheck{"Dual-trigger?
(borderline / poor)"}
    Margin -->|"Tied"| N3(["⚪ Neither"]):::neither
    ConfCheck -->|No| SWDone{"Which allele?"}
    ConfCheck -->|Yes| LocalSW["Local SW alignment
(soft-clips bad flanks)"]
    LocalSW --> LocalDone{"Score diff ≥ 2?"}
    LocalDone -->|ALT| AltL(["🔴 ALT"]):::alt
    LocalDone -->|REF| RefL(["✅ REF"]):::ref
    LocalDone -->|"Tied"| N2(["⚪ Neither"]):::neither
    SWDone -->|ALT| AltS(["🔴 ALT"]):::alt
    SWDone -->|REF| RefS(["✅ REF"]):::ref

    classDef entry fill:#9b59b6,color:#fff,stroke:#7d3c98,stroke-width:2px;
    classDef ref fill:#27ae60,color:#fff,stroke:#1e8449,stroke-width:2px;
    classDef alt fill:#e74c3c,color:#fff,stroke:#c0392b,stroke-width:2px;
    classDef neither fill:#95a5a6,color:#fff,stroke:#7f8c8d,stroke-width:2px;

Phase 1: Haplotype Reconstruction¶

Walks the CIGAR to rebuild what the read shows for the genomic region [pos, pos + ref_len). Each CIGAR operation contributes differently:

Phase 2: Masked Comparison¶

Instead of exact matching, bases with quality below --min-baseq or N bases are masked out — they cannot vote for either allele. Three cases based on reconstructed sequence length:

Phase 2.5: Edit Distance¶

When the reconstruction's length matches neither allele (so no Case A/B/C comparison is possible), gbcms measures the Levenshtein edit distance between the reconstruction and each allele. It runs only for small alleles (REF and ALT both ≤50bp, ref_len ≤ 2×alt_len). A failed Case A/B comparison never reaches this phase — it returns Neither (with partial_alt when reliable ALT-matching bases exist) or falls to Phase 3. This catches cases where the reconstruction is off by 1-2 bases due to an incomplete variant definition:

Phase 3: Alignment Fallback¶

When Phase 2 and Phase 2.5 both fail, the engine expands to the full ref_context window and performs dual-haplotype alignment against both REF and ALT haplotypes.

Default backend: pairhmm (two-stage WFA + PairHMM)

1. Pangenomic haplotype matrix — build REF_hap and ALT_hap extended with flanking context
2. WFA fast-path triage (wfa2lib-rs) — edit-distance alignment resolves ~70-80% of reads instantly at O(s²) cost; if the edit-distance gap between REF and ALT is unambiguous, the read is classified immediately
3. Marginalized PairHMM (escalated only when WFA is ambiguous) — integrates per-base quality probabilities into alignment scoring, producing a Log-Likelihood Ratio (LLR) confidence score; default threshold 2.3 ≈ ln(10) — the read is aligned whole against each haplotype, whose ends are free, as the sw backend's semiglobal alignment does: every read base is explained, and haplotype bases the read does not reach cost nothing

Alternative backend: sw

Runs Smith-Waterman directly on every Phase 3 read (no WFA pre-filter):

1. Build haplotypes: REF_hap = left_ctx + REF + right_ctx, ALT_hap = left_ctx + ALT + right_ctx
2. Mask: Replace low-quality bases with N (scores 0 against any base)
3. Align: Semiglobal alignment — read (query) is fully consumed, haplotype (text) has free overhangs
4. Score: ALT wins if alt_score ≥ ref_score + 2; REF wins if ref_score ≥ alt_score + 2
5. Dual-trigger check: For indel/complex variants, if the semiglobal result is low-confidence, retry with local alignment

Dual-Trigger Local Fallback¶

For complex variants (both alleles > 1bp with different lengths), semiglobal alignment can produce confident but incorrect calls when the MAF/VCF definition is incomplete (e.g., a complex variant missing an adjacent SNV). The ALT haplotype then has a "frameshifted flank" — right-context bases that don't match the biological read. Semiglobal forces gap penalties through this invalid flank.

[!NOTE] The dual-trigger only applies to complex variants (e.g., TCC→CT, ATGA→CATG), not to pure insertions or deletions. Pure indels are well-handled by semiglobal alignment, and applying local fallback would risk false positives in homopolymer regions.

Two conditions detect low-confidence semiglobal results:

When either trigger fires, the engine retries with local alignment (Aligner::local()), which soft-clips the bad flank and finds the best matching substring without penalizing overhangs on either side. The call then rests on the local scores, and so do the checks after it: a REF call's nearby evidence (ALT within twice the margin: partial_alt) and a tie's (the best score reaching half the read) (#195).

Only complex rows the exact-carrier rule cannot judge reach this step (direct API use, a failed reference fetch, a contig end). On every arm measured on 2026-10-30 (about 4.5 million read calls under each backend) none did.

Ambiguous Tie — Unbiased VAF Preservation¶

For small Complex and MNP variants, biological reads heavily affected by surrounding genetic polymorphism can result in 50% partial matches against the Alternate array. Mathematically, this scores an exact numerical tie between the REF and ALT haplotypes (alt_score = 11, ref_score = 11).

These ambiguous reads are routed to neither (is_ref = false, is_alt = false). They still contribute to physical Total Depth (DP) via the anchor overlap gate, but do not inflate RD or AD. This preserves an unbiased VAF = AD / (RD + AD) — critical for low-VAF cfDNA detection where even small RD inflation can push a variant below the clinical Limit of Detection.

Informative Reads for Indels¶

An indel inside a repeat can sit anywhere along its shift-equivalence region: the stretch it slides over without changing the haplotype (the whole tract for a homopolymer, every full and partial copy for an STR, the event itself in unique sequence). A read that starts or ends inside that region reads the same with or without the event, so the aligner places no gap and IGV cannot show which allele it carries. Counting such a read as REF biases VAF down. For a 1bp deletion in a 10-base homopolymer, with 10 REF and 10 ALT reads spanning the tract and 10 of each ending inside it, VAF came out 25% instead of 50%.

A REF call on a pure insertion or deletion therefore needs an informative read. Reading inward from one flank of the region, REF and ALT agree until one base, the first discriminating base. The read must cover that flank base and run one base past the first discriminating base on that side. Either side will do.

The extra base past the discriminating one is an aligner margin. An ALT read ending on the discriminating base keeps a single terminal mismatch (cheaper than a clip) and would look like REF. With one more base, the mismatches cost more than the gap or the clip.

Prep measures each pure indel's region over its own reference fetch, sized to the event (Variant.shift_region). A tandem duplication (an ITD, for example) slides over its whole duplicated segment, which can be far longer than the repeat context kept for alignment.

A read spans a window with one aligned block, between its splices (RJ-17, 6.6.0): a splice is not reference coverage, so an RNA read spliced inside the window shows none of it past the splice, as a read ending there shows none. An exon piece that a junction at a pure deletion's edge begins and the read's end closes inside the deleted bases holds part of them and no flank base of its own: it is depth only, never REF, with BAQ on or off (amended 2026-10-30,

292).¶

Reads that fail the rule count toward DP and DPF but are neither REF nor ALT, the same as GATK's AD, which counts only informative reads. The rule does not apply to substitution-bearing events (SNV, MNP, delins): they have no shift region, and their first base already discriminates.

A worked example: counting by base against by allele (#253, measured 2026-10-08). A 9bp deletion in a 51bp trinucleotide repeat, in an IMPACT tumor and its matched normal. By base at the indel's position, as IGV, bam-readcount, the original GetBaseCounts and the clinical counts do: 159 reads show the reference base, but only 16 of them cross the tract. The deletion count of 24 is every length together; the given 9bp deletion is 15 of them. By allele, the tumor's vaf is 15/31 ≈ 48% and the normal's crossing reads show 0/23 (a working note of the same measurement says 0/24). So it is a high-fraction somatic event, which by-base counting reports at 6.6%. Whether allele-level VAF needs a spanning-bias correction (the shorter allele crosses the tract more often) is measured against germline hets in C41 #264.

The ALT side (6.6.0, #188). An ALT call stands from the CIGAR alone when the read spans C10's windows read on the ALT haplotype, [lo − 1, hi + 2) or [lo − 2, hi + 1). For an insertion these are its REF windows. For a deletion they are an insertion's, because a carrier's reference extent counts its own gap: spanning the REF windows, a −AA carrier in a run of six A's may hold only G AAAA, which fits both alleles. Otherwise the call stands only when the read's own bases tell the alleles apart: read where they sit, rightwards from its aligned base left of the window or leftwards from its aligned base right of it, the read must read, unmasked, a base where the alleles differ, every base up to it fitting the ALT. There is no margin base past it (C10's margin guards CIGAR-only REF calls against a hidden terminal mismatch; here the deciding base itself is read, on a read the CIGAR already calls ALT). A carrier that ends inside the repeat holds only shared bases, so its gap alone is placement: it counts toward depth only. A long insertion's carrier spends its span inside the insert and often cannot reach the far flank, yet its inserted bases discriminate, so its bases decide. A masked base where the alleles first differ is skipped: a later unmasked one where they differ decides. Past the reference fetched around the event a base decides nothing.

One change across the window (6.6.0, #188). An indel at the variant's own junction counts ALT only when the read carries no other insertion or deletion across the discrimination window; a cancelled pair or a split longer allele (+A and +A for a +A row) is a distinct allele, partial evidence. The order the aligner writes an I/D pair in does not matter: a deletion right after an insertion (M I D M), or an insertion right after a deletion (M D I M), is inspected as one right after an aligned block is; a deleted anchor followed by an insertion is judged by the read's bases (a re-inserted anchor with the insert is the ALT). A same-length insertion of other bases beside another indel in the window is another allele, not Phase 3's to arbitrate. At 50bp or more the junction judge reads the read's bases at the POS breakpoint (RJ-24). Splices there remain the RNA rules' business.

Multi-Allelic Behavior¶

When multiple variants have overlapping REF spans at the same locus, reads carrying one variant's ALT allele could be incorrectly counted as REF for another variant. The engine addresses this with a two-phase approach:

Phase 1: Annotation¶

During normalization, assign_multi_allelic_groups() groups co-annotated variants with a fixed-point sweep over sorted (chrom, pos) coordinates, under two criteria: variants whose REF spans intersect (any types — tagged MULTI_ALLELIC), and length-changing variants whose scan windows (max(5, repeat_span+2) each side) overlap (tagged TRACT_CLUSTER when window-only). Groups close transitively and may be non-contiguous in position order; members share a multi_allelic_group ID and the tag is appended to gbcms_status_reason (the verdict stays PASS). Within a group the engine assigns AD exclusively: a read's ALT call is demoted to partial_alt when a sibling explains it strictly better by span-explanation cost, or exactly as well as the row's ALT does (the two alleles differ only at bases the read has masked, so the read is either), and an alignment-phase ALT not confirmed exactly by the read's own span reconstruction is demoted as ambiguous; anchor-exact CIGAR evidence is never contested. An ALT call the exact-carrier rule made is contested by its bases instead: the whole read is placed ungapped anywhere on this row's ALT and on each sibling's, and a sibling placing it with as few confident mismatches or fewer claims it; the aligned REF re-test is not run for it (RJ-5, amended 2026-10-22).

One allele given on two rows is not its own sibling: the same row given verbatim (a cohort MAF listing a recurrent indel once per sample), or one allele written two ways (a repeat unit inserted at either end of the repeat, which normalization makes one variant). Before 6.6.0 each twin was the other's: an indel carrier the aligner wrote anywhere but the junction fit both exactly, so it counted at neither row, while one written at the junction counted at both (and --rescue-mnp skipped an MNP given twice as grouped). Each row's sibling list now holds one entry per other allele and never its own, so each row counts every carrier; the rows keep their group and tag. A row given verbatim is noted at INFO; one allele given two ways warns, naming the rows.

Phase 2: Sibling ALT Exclusion¶

During counting, a read classified as REF for a variant is checked against the siblings whose change lies inside the variant's discrimination window: its shift-equivalence region (see Informative Reads for Indels; the span for other variants) plus one base on each side. For each such sibling, the full check_allele_with_qual() pipeline determines whether the read carries the sibling's ALT. If it does, the read is excluded from REF for the current variant and counted as partial_alt.

A carrier of a sibling whose change lies outside the window shows the reference across every base that could tell this variant's alleles apart, so it stays REF, as IGV shows it and as GATK counts REF at a site unless an event overlaps it.

The exclusion runs before fragment evidence is recorded, so ref_count and ref_count_fragment drop the same molecules. The main and per-transcript read loops apply the same rule.

SW Gap Penalties¶

Phase 3's Smith-Waterman aligners use fixed affine gap penalties: SW_GAP_OPEN = -6, SW_GAP_EXTEND = -1 (match +1, mismatch −1, N scores 0). bio 4 counts a gap's first base in the open penalty (a k-base gap scores SW_GAP_OPEN + SW_GAP_EXTEND·(k − 1), rust-bio#660), so −6 is bio 3's −5 open plus one extension: a k-base gap scores −5 − k, as it did. They are documented constants, not tuned per locus: an earlier dynamic_sw_gap_extend logistic curve rounded to −1 for every repeat_span, and traced real runs (ACCESS duplex, MSI-high) confirmed SW scores nothing under the default backend on well-formed input, so the curve was removed (issue #92). The traced runs were an ACCESS duplex sample and the heaviest MSI-high sample: 42 tract indels and 31 Phase-3 classifications, with zero Smith-Waterman fallback invocations (#92 comment). Under the default PairHMM backend the pangenomic matrix builds on well-formed input (prep supplies ref_context), and on the explicit SW backend the old function was a constant −1 anyway. Making the relaxation real (a curve reaching 0 for deep repeats) was not adopted: it would change SW alignment scores on a path that never fires, at the cost of a validation pass (REJ-20260923-005 in .agents/learnings/REJECTED.md).

Smith-Waterman has two roles:

• The explicit scorer under --alignment-backend sw — the primary Phase-3 engine there, kept for cross-backend concordance. There are no SW gap flags.
• A last-resort fallback under the default pairhmm backend, used only when the pangenomic haplotype matrix cannot be built for a variant (its reference context is missing because prep's fetch failed, does not contain it, or the ALT haplotype cannot be constructed) — i.e. upstream input was malformed. The fallback is kept but never silent: each affected variant logs one WARN naming the reason and outcome, and the row carries SW_FALLBACK(n) in gbcms_diagnostic (n = depth-contributing reads the matrix could not evaluate). Where SW cannot build its haplotypes either — always the case without a reference context — those reads end as NEITHER, and the flag is how that loss becomes visible. (MNPs carry no reference context by design and use no Phase 3, so their reads are never counted here.)

The PairHMM gap flags retune PairHMM probabilities only; they never reach the SW penalties.

The repeat_span is computed during normalization using find_tandem_repeat() at the first changed base of the alleles (not the shared anchor, which sits one base left of a left-aligned tract) and stored on the Variant struct. It drives the PairHMM gap blending, the windowed-scan width, and adaptive context padding (the SW penalties are constants and do not use it).

Limitations¶

1. Windowed scan range — Indels shifted beyond the context padding from their expected position won't be detected by the CIGAR-based check. Phase 3 can catch some of these via ref_context, but only if the read shows evidence (indels/clips) near the variant. Adaptive context padding (enabled by default) utilizes a multi-anchor footprint sweep to dramatically widen the padded detection bounds natively for INDEL clusters.

2. Score margin ≥ 2 — The SW margin is fixed at 2 points to prevent ambiguous calls. Reads failing to achieve definitive spacing are routed to neither — they contribute to DP but not to RD or AD, preserving unbiased VAF.

3. Soft-clip recovery — Phase 1 includes soft-clipped bases that overlap the variant window, but only when ref_pos is within the variant region. Soft clips at the edge of reads far from the variant are not considered.

4. MNP strict matching — MNP evaluation is per-position and terminal: LowQuality and ThirdAllele results return neither (with partial_alt tracking) and are deliberately not rescued by the complex chain. Only structural issues (no coverage, indel inside the block) fall back to check_complex.

5. Incomplete variant definitions — When the MAF/VCF represents a complex event incompletely (e.g., TCC→CT omitting an adjacent SNV), reads may carry a different CIGAR signature than expected. The dual-trigger local fallback in Phase 3 mitigates this by soft-clipping "frameshifted flanks" caused by the definition mismatch.

Related¶

• Variant Normalization — How variants are prepared before classification
• Read Filters — Which reads reach the allele checker
• Counting & Metrics — How classifications become counts
• Glossary — Term definitions