Counting & Metrics¶

How allele classifications become counts — read-level metrics, fragment counting, strand bias, and output columns.

Overview¶

After each read is classified as REF, ALT, or neither, gbcms accumulates counts at two levels: individual reads and collapsed fragments. It then computes strand bias statistics using Fisher's exact test.

flowchart LR
    Class(["🧬 Classification"]):::start

    subgraph ReadLevel ["📖 Read-Level Path"]
        direction LR
        RL1["Read Counts\nDP / RD / AD"] --> SB["Strand Bias\n(Fisher's test)"]
    end

    subgraph FragLevel ["🔗 Fragment-Level Path"]
        direction LR
        FL1["Fragment Tracking\n(QNAME hash)"] --> FL2["Quality Consensus\n(R1 vs R2 resolve)"] --> FL3["Fragment Counts\nDPF / RDF / ADF"] --> FSB["Fragment Strand Bias\n(Fisher's test)"]
    end

    Class --> ReadLevel
    Class --> FragLevel

    classDef start fill:#9b59b6,color:#fff,stroke:#7d3c98,stroke-width:2px;
    classDef metrics fill:#3498db,color:#fff,stroke:#2471a3,stroke-width:2px;
    class SB,FSB metrics;

Read-Level Counting¶

Each read passing filters is counted independently. No deduplication is applied at this level.

flowchart TD
    Fetch(["📥 Fetch reads ±5bp window"]) --> Classify["🧬 Allele classification"]
    Classify --> Anchor{"Read start ≤ variant POS?"}
    Anchor -->|"Yes — normal"| DP["✅ Count in DP"]
    Anchor -->|"No"| ClassCheck{"Classified as REF or ALT?"}
    ClassCheck -->|"Yes — shifted indel"| DPShifted["✅ Count in DP"]:::shifted
    ClassCheck -->|"No"| Skip(["⏭️ Skip — outside footprint"]):::skip
    DP --> Allele{"Allele?"}
    DPShifted --> Allele
    Allele -->|REF| RD["RD += 1"]:::ref
    Allele -->|ALT| AD["AD += 1"]:::alt
    Allele -->|Neither| Neither["DP only (no RD/AD)"]:::neither

    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 skip fill:#bdc3c7,color:#000,stroke:#95a5a6;
    classDef shifted fill:#e67e22,color:#fff,stroke:#bf6516,stroke-width:2px;

Read Metrics¶

Fragment Counting¶

Fragment counting collapses read pairs into a single observation per fragment. This is critical for cfDNA sequencing (like MSK-ACCESS) where the same DNA fragment is sequenced from both ends — counting reads would double-count each fragment.

Fragment Tracking¶

Fragments are tracked using QNAME hashing — each read's QNAME is hashed to a u64 key for memory-efficient lookup. The FragmentEvidence struct records the best base quality seen for REF and ALT across both reads in the pair, plus a sticky has_structural_alt flag that captures whether any read had a direct CIGAR I/D op matching the target INDEL.

Phase A — Fragment Tracking (per-read pass):

flowchart LR
    Read(["📖 Read passes filters"]):::start --> Hash["Hash QNAME → u64 key"]
    Hash --> Seen{"Fragment seen before?"}
    Seen -->|No| Create["Create FragmentEvidence\n{ref_qual=0, alt_qual=0}"]
    Seen -->|Yes| Update["Update best quality score"]
    Create --> Store["Store allele + base quality"]
    Update --> Store
    Store --> Next(["→ next read …"]):::next

    classDef start fill:#3498db,color:#fff,stroke:#2471a3,stroke-width:2px;
    classDef next fill:#3498db,color:#fff,stroke:#2471a3,stroke-width:2px;

Phase B — Quality-Weighted Consensus (after all reads loaded):

flowchart TD
    ForEach(["For each unique fragment"]):::start --> HasBoth{"Both REF and ALT evidence?"}
    HasBoth -->|No| Direct["Assign to whichever allele was seen"]
    HasBoth -->|Yes| Structural{"Structural INDEL evidence?\n(CIGAR I/D op on ALT read)"}
    Structural -->|Yes| StructAlt(["🔴 Count as ALT\n(structural priority)"]):::structural
    Structural -->|No| QualCheck{"Quality difference > threshold?"}
    QualCheck -->|"REF qual >> ALT"| Ref(["✅ Count as REF"]):::ref
    QualCheck -->|"ALT qual >> REF"| Alt(["🔴 Count as ALT"]):::alt
    QualCheck -->|"Within threshold"| Discard(["⚪ Discard — ambiguous"]):::discard
    Direct --> Count
    StructAlt --> Count
    Ref --> Count
    Alt --> Count
    Discard --> DPF["Counted in DPF only"]:::info
    Count["Update RDF / ADF / DPF + strand counts"]

    classDef start 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 structural fill:#9b59b6,color:#fff,stroke:#7d3c98,stroke-width:2px;
    classDef discard fill:#95a5a6,color:#fff,stroke:#7f8c8d,stroke-width:2px;
    classDef info fill:#3498db25,stroke:#3498db;

Quality-Weighted Consensus¶

When R1 and R2 of a fragment disagree (one supports REF, the other ALT), the engine resolves the conflict:

The threshold is configurable via --fragment-qual-threshold (default: 10).

Fragment Orientation¶

Fragment strand is determined by read 1 orientation (preferred) or read 2 if read 1 is not available. This is consistent with standard library preparation conventions.

Fragment Metrics¶

VAF Calculation¶

VAF = AD / (RD + AD)

Where AD and RD are the read-level alternate and reference counts. Fragment-level VAF can be similarly computed as ADF / (RDF + ADF).

Strand Bias (Fisher's Exact Test)¶

Strand bias is detected by building a 2×2 contingency table from strand-specific counts and running Fisher's exact test:

flowchart LR
    Counts(["2×2 Strand Counts"]):::start --> Fisher["Fisher's Exact Test"]
    Fisher --> PVal["SB_pval"]
    Fisher --> OR["SB_OR"]
    PVal --> Interp{"p < 0.05?"}
    Interp -->|Yes| Bias(["⚠️ Possible Strand Bias"]):::warn
    Interp -->|No| NoBias(["✅ No Strand Bias"]):::pass

    classDef start fill:#3498db,color:#fff,stroke:#2471a3,stroke-width:2px;
    classDef warn fill:#e74c3c,color:#fff,stroke:#c0392b,stroke-width:2px;
    classDef pass fill:#27ae60,color:#fff,stroke:#1e8449,stroke-width:2px;

Computed at both levels:

The p-value is the exact two-sided Fisher test as R's fisher.test defines it: the sum of the probabilities of every table with the observed margins that is no more likely than the observed table (with R's tie tolerance, 1 + 10⁻⁷). It is computed in log space, so it is exact at any depth. Before 6.6.0 it was 0 once both strands were deep (from about 1,030 reads on a strand-balanced table), and strongly biased tables were floored near 10⁻¹⁰.

Complete Output Column Reference¶

The core counting fields of the BaseCounts struct returned by count_bam_binned() — the binned counting entry point that groups variants into ~10kb genomic bins for a single bam.fetch() per bin before classifying reads. See Architecture → Genomic Binning for how bins are built. (The gated mFSD, RNA, and ASJD fields — plus QC fields like mq0_count, alt_dist_end_median/ref_dist_end_median, and non_discriminating_locus — are documented in their own sections above and in the output-formats reference.)

Diagnostic Column Invariants¶

The diagnostic columns maintain strict structural invariants:

VCF INFO/FORMAT Tags¶

N-Base Handling¶

N bases arise from:

• Duplex masking (fgbio): disagreeing R1/R2 bases → masked to N (BQ ≈ 2)
• Sequencer failure: uncalled bases

N bases are strictly uninformative — they do NOT contribute to RD, AD, any_alt, or partial_alt. They are tracked separately via n_count to enable:

• QC: n_count / DP ratio flags duplex masking hotspots
• mFSD: N-classified fragments form a separate size class for distribution analysis

mFSD — Mutant Fragment Size Distribution¶

mFSD compares insert-size distributions for REF-classified vs ALT-classified fragments at each variant position. Short-fragment enrichment in the ALT class is graded evidence toward a tumor origin, never an origin call (mFSD report classes).

Enabled with --mfsd. All 41 columns (and their VCF INFO equivalents) are absent from output without this flag.

Fragment Classes¶

A fragment whose reads all start or end inside an indel's repeat tract carries no readable allele (informative reads). It counts in fragment depth but in none of the four classes.

Fragment size¶

A fragment's size is its molecule's length: TLEN (the reference span of the pair) less the deletions and splices, plus the insertions, its reads' CIGARs write, the smallest over its reads. An ALT fragment is the molecule carrying the ALT: when its span covers the variant but no read's CIGAR writes the allele's length change there (a carrier written as a soft clip, or judged by its bases), the change is applied too, so a deletion carrier is not sized longer by the deleted bases (since 6.7.0; before, a fragment of a deletion over about 830 bases fell outside the 1000-bp limit). A carrier aligned from past a deletion never spanned the deleted bases, and keeps its TLEN.

MAF Columns (41 total)¶

Raw Counts¶

Log-Likelihood Ratios¶

Mean Fragment Sizes¶

Pairwise KS Statistics (6 pairs × 3 values = 18 columns, + 1 FDR q-value)¶

Each pair yields: delta (mean difference in bp), ks (KS D-statistic), pval (KS p-value). When either class has fewer than 5 fragments the test does not run (mfsd_ks_valid = False): ks is NA and pval is 1.0000, a placeholder; delta is still the mean difference when both classes have fragments. Read mfsd_ks_valid before the KS columns.

The pval is the exact two-sample KS p-value whenever the two classes span at most 10⁷ lattice cells (n·m) — every class pair a targeted panel produces — and the asymptotic Kolmogorov series with Stephens' finite-sample correction above. Exact matters where cfDNA decides: a few ALT fragments against thousands of REF fragments, where the uncorrected series overstates p (about 1.7x at 5 ALT fragments near p = 0.05, and up to 45x for a strong shift). Fragment sizes are integers, so ties are common; the exact p-value treats sizes as continuous, which is conservative with ties. Below a handful of fragments the test has little power (about 8% at 5 ALT fragments on real cfDNA), so a non-significant p is not evidence of "no shift".

An additional column, mfsd_qval_alt_ref, carries the Benjamini-Hochberg FDR q-value for the ALT-vs-REF KS p-value, corrected across all variants with a valid ALT-vs-REF test in the sample. The mFSD report's LEANS-SOMATIC class uses this q-value, not the raw p-value. When the KS test did not run it stays the p-value placeholder (1.0000) and is left out of the correction; a variant on a contig absent from the BAM is likewise left out.

Derived Metrics¶

Nucleosomal Fractions¶

CH Gene Flag¶

VCF INFO Fields (13 total)¶

Added to ##INFO header and per-variant INFO column when --mfsd is set.

Parquet Output (--mfsd-parquet)¶

When --mfsd-parquet is also set, a companion <sample>.fsd.parquet is written with raw arrays:

Written natively by the Rust engine (no pyarrow dependency).

RNA-Specific Output Columns¶

When running in RNA mode (gbcms rna), additional columns capture transcriptome-specific metrics. These columns are absent from DNA mode output.

Biological Context¶

RNA-seq reads exhibit orientation biases (dUTP strandedness), splice junctions (CIGAR N operations representing intron skips), and post-transcriptional modifications (A-to-I editing by ADAR enzymes). These metrics enable downstream filtering of RNA-specific artifacts and biological characterization of variants.

MAF Columns (5 additional)¶

VCF Fields¶

UMI-Aware Fragment Counting¶

When --umi-tag is specified (e.g., --umi-tag RX), fragment identity incorporates the UMI barcode in addition to the QNAME:

Comparison with Original GBCMS¶

For the measured count differences from GetBaseCountsMultiSample 1.2.5, one line each with its cause, see Known differences from 1.2.5.

Related¶

• Output Formats — Column schemas for VCF and MAF output (DNA and RNA)
• Allele Classification — How reads are classified
• Read Filters — Which reads reach counting
• Architecture — System design
• Glossary — Term definitions