Skip to content

gbcms dna

Count alleles at variant positions across one or more DNA/cfDNA BAM/CRAM files.

Migrated from gbcms run

The deprecated gbcms run command was removed in v4.1.0. Use gbcms dna instead — all arguments are identical.

Synopsis

gbcms dna [OPTIONS] --variants <FILE> --bam <NAME:PATH>... --fasta <FILE> --output-dir <DIR>

--variants, --fasta and --output-dir are required, and so is one of --bam or --bam-list; a command without them exits with a usage error.

Required Arguments

Option Description
--variants, -v VCF or MAF file with variant positions (.vcf, .vcf.gz, .vcf.bgz, or .maf). Unsupported extensions are rejected immediately.
--bam, -b BAM or CRAM file path (can repeat). Optionally prefix with name: for sample naming, e.g. --bam tumor:tumor.bam. If no name given, the filename stem is used.
--bam-list, -L File containing BAM/CRAM paths (one per line, optionally sample_name path). Alternative to repeated --bam.
--fasta, -f Reference FASTA file (with .fai index). Required for all runs; also used for CRAM decoding.
--lenient-bam Skip missing --bam paths and missing --bam-list entries, continuing with the remaining samples (default: all missing paths are reported and the run exits non-zero). A missing or unreadable --bam-list file, or a duplicate sample name, always fails regardless.

CRAM Support (v5.3.0+)

--bam accepts both BAM (.bam) and CRAM (.cram) files transparently. The --fasta reference is automatically used for CRAM decoding at both the Rust engine and Python header-validation layers. No additional flags are needed.

Index discovery is automatic: gbcms searches for .bam.bai / .bai (BAM) or .cram.crai / .crai (CRAM) next to the alignment file. You can also provide the index path explicitly via the Nextflow samplesheet bai column.

Output Options

Option Default Description
--output-dir, -o required Output directory
--format vcf Output format (vcf or maf)
--suffix '' Suffix for output filenames
--column-prefix '' Prefix for gbcms count columns in MAF output. Only letters, digits, and underscores ([A-Za-z0-9_]) are allowed; invalid characters exit immediately.
--preserve-barcode false Keep original Tumor_Sample_Barcode from input MAF. No-op (with warning) when input is not MAF.
--show-normalization false Append norm_* columns showing left-aligned coordinates
--context-padding 5 Minimum flanking bases for haplotype construction. Range 1–50, enforced at parse time. Auto-increased in repeat regions when --adaptive-context is enabled.
--adaptive-context true Dynamically increase context padding in tandem repeat regions
--threads 1 Total worker-thread budget for this run (minimum 1; 0 is rejected). All parallelism — read classification and any decode threads — stays within this budget, so gbcms never oversubscribes a small allocation. Multi-sample parallelism is the orchestrator's job (e.g. N concurrent Nextflow processes, each pinned to task.cpus), not this flag's.

mFSD Options

Mutant Fragment Size Distribution (mFSD) analysis compares insert-size distributions for REF- vs ALT-classified fragments at each variant position, enabling detection of short-fragment enrichment associated with tumor-derived cfDNA (see mFSD Metrics).

Option Default Description
--mfsd false Enable mFSD analysis. Adds 41 mFSD columns (KS test, LLR, mean sizes, pairwise comparisons, derived metrics) to MAF output and 13 MFSD_* INFO fields to VCF.
--mfsd-parquet false Write a companion <sample>.fsd.parquet with per-variant raw fragment size arrays (ref_sizes, alt_sizes). Enables downstream visualizations. Requires --mfsd.
--mfsd-report false Generate an interactive HTML report (<sample>.mfsd_report.html) with per-variant fragment size distributions, dual-axis histograms, and graded fragment-size evidence (leans somatic, no size evidence, insufficient). Implies --mfsd and --mfsd-parquet (both are auto-enabled). Requires --format maf: the report reads its statistics from the MAF, so --format vcf is refused. See mFSD Report.
--mfsd-report-min-alt 3 Minimum ALT fragment count to include a variant in the HTML report.
--mfsd-report-max-variants 20 Maximum variants in the HTML report (selected by highest ALT count). Use -1 for no limit.

Tip

To generate summary statistics, raw Parquet data, and an interactive HTML report: 

gbcms dna --mfsd-report --format maf \
    --variants variants.maf --bam tumor:tumor.bam --fasta hg19.fa -o ./results
This produces <sample>.maf (with 41 mFSD columns), <sample>.fsd.parquet (raw fragment sizes), and <sample>.mfsd_report.html (interactive visualization). Using --mfsd-report auto-enables --mfsd and --mfsd-parquet.

Observation Options

Flag Default Description
--observations-parquet false Write a companion <sample>.observations.parquet with the per-molecule allele call at each variant (which molecule carried REF/ALT/N/OTHER). Enables read-backed phasing and allelic-imbalance analyses downstream. Counts output is unchanged. See Per-Molecule Observations.

Rows are written from Rust, so memory stays flat at panel scale. Nothing is written unless the flag is passed.

Filtering Options

Option Default Description
--min-mapq 20 Minimum MAPQ
--min-baseq 20 Minimum base quality. A base below it is read as unknown, like N: it supports no allele and counts against none. Reads are not dropped: a read still counts toward depth and is judged on its other bases. 0 masks only N bases. See Read Filters → Quality Thresholds
--filter-duplicates true Filter duplicate reads
--filter-secondary true Filter secondary alignments
--filter-supplementary true Filter supplementary alignments
--filter-qc-failed true Filter QC failed reads
--filter-improper-pair false Filter improperly paired reads
--filter-indel false Filter reads with indels
--fragment-qual-threshold 10 Quality difference threshold for fragment consensus (see Fragment Counting)

BAQ Options

Base Alignment Quality (BAQ) heuristically downgrades base qualities near a read's insertions and deletions to prevent systematic errors from realignment artifacts.

Option Default Description
--apply-baq/--no-baq off Enable BAQ quality downgrade near indels

Internal constants (not configurable — derived from Li 2011):

Constant Value Description
BAQ_RADIUS 5 bp Bases within 5 bp of an indel boundary are penalized
BAQ_PENALTY 20 Subtracted from BQ (clamped to 0)

When to Enable BAQ for DNA

Enable --apply-baq when your upstream pipeline does not include BQSR, fgbio consensus, or any other base quality recalibration. In these "raw BQ" pipelines, bases near indels may retain inflated sequencer-assigned quality scores. BAQ compensates by penalizing BQ within 5 bp of indel boundaries.

Do Not Enable for Pre-Calibrated BAMs

If your DNA pipeline includes BQSR or consensus calling (e.g., fgbio, Marianas), leave BAQ off (the default). Applying BAQ on already-recalibrated BAMs may over-penalize, causing undercounting of legitimate variant evidence.

RNA Mode Default

In gbcms rna, BAQ is on by default because RNA pipelines typically lack BQSR. A splice junction (CIGAR N) lowers no base's quality. See RNA Splice-Junction Handling.

UMI Options

Unique Molecular Identifier (UMI) support for molecule-level deduplication.

Option Default Description
--umi-tag (none) BAM tag for UMI barcode (e.g. RX). Enables UMI-aware fragment grouping.

UMI-Aware Fragment Counting

When --umi-tag is set, two reads are considered the same fragment only if they share both QNAME and UMI barcode. This prevents UMI-collapsed reads from different original molecules being incorrectly merged into a single fragment, which would deflate fragment-level allele counts.

gbcms dna --umi-tag RX \
    --variants variants.vcf --bam sample.bam --fasta ref.fa -o results/

MNP Rescue Options

The rescue pass is for annotated MNPs (multi-nucleotide polymorphisms) whose carriers hold only a component of the haplotype — e.g. two SNVs on different molecules annotated as one MNP. gbcms correctly reports such a haplotype as absent: the carriers land in partial_alt, not alt_count. Rescue re-counts each discriminating position as an SNV and reports the best-supported component instead. See Architecture → MNP Rescue Pass for the full design.

Option Default Description
--rescue-mnp false Enable the rescue pass. Candidates are PASS MNPs flagged MNP_RESCUE_ELIGIBLE with partial_alt > alt_count, outside co-annotated groups, and whose full haplotype no read shows. When the best component beats the MNP's alt_count, that component's full counts (every count, fragment, strand and mFSD column) replace the row's, the row is flagged RESCUED_COMPONENT(chrom:pos:REF>ALT) in gbcms_diagnostic, and a warning is logged; the MNP's own counts and the per-position split go to gbcms_rescue. MAF and VCF output carry the same counts, flags and audit.
--rescue-mnp-threshold 1.0 Maximum discriminating/length ratio for MNP rescue eligibility (0.0–1.0). 1.0 = all MNPs are eligible (C++ gbcms compatible, default). 0.5 = conservative sparse-only mode (≤50% discriminating positions). 0.0 = disable rescue eligibility (MNP_DISC_RATIO diagnostics are still emitted). Only used when --rescue-mnp is enabled.

Diagnostic Flags

Every MNP row carries MNP_DISC_RATIO(n/m), with or without --rescue-mnp, and MNP_RESCUE_ELIGIBLE when n/m is at most --rescue-mnp-threshold. Both describe the annotated MNP's shape, so a rescued row keeps them next to RESCUED_COMPONENT. Definitions: QC Flags → Diagnostics.

A rescued row reports a component, not the annotated MNP

When reads show the whole MNP (e.g. a somatic change on top of a germline SNP), rescue keeps the MNP's counts (outcome=haplotype_confirmed). But where the MNP itself is absent — typically a fillout of other timepoints or normals — a germline SNP component can still be adopted, and the rescued VAF is then the germline VAF. Read gbcms_rescue (outcome, original counts, per-position split) before interpreting a rescued row.

Choosing the threshold

# Permissive (default) — rescue all MNPs, matching C++ gbcms behavior
gbcms dna --rescue-mnp --variants input.maf --bam sample:sample.bam --fasta ref.fa -o out/

# Conservative — rescue only sparse MNPs (≤50% discriminating positions)
gbcms dna --rescue-mnp --rescue-mnp-threshold 0.5 --variants input.maf --bam sample:sample.bam --fasta ref.fa -o out/

Homopolymer Twin Option

Option Default Description
--rescue-homopolymer false Deprecated. Deprecated in 6.7.0; removed in 6.8.0; it did nothing useful the diagnostics do not (OBSERVED_ALLELE and COEXISTING_ALLELE). Accepted and ignored, with a warning; the row counts the given allele.

Off by default: the row counts the given allele

Both forms accept near-matches, so the winner can report another allele's reads under the row's label. By default gbcms counts the allele it is given, and when the reads carry a different one, gbcms_diagnostic names it: OBSERVED_ALLELE when no read carries the given allele, COEXISTING_ALLELE when it is present beside a more frequent one (see output formats).

Debugging Options

Option Default Description
--verbose, -V false Enable verbose debug logging
--trace, -T false Enable per-read Rust trace logging (slow). Implies --verbose. Shows detailed per-read classification diagnostics.

Alignment Backend

Phase 3 (haplotype-based) classification uses a two-stage pipeline:

  1. WFA fast-path (Wavefront Alignment, wfa2lib-rs) — edit-distance triage against the pangenomic haplotype matrix. Resolves ~70-80% of reads instantly at O(s²) cost where s is the edit distance. If REF and ALT scores differ clearly, the read is classified immediately.

  2. Marginalized PairHMM (escalated only when WFA is ambiguous) — integrates per-base quality probabilities into alignment scoring, producing a log-likelihood ratio (LLR) confidence score. More sensitive in noisy, low-quality, or repeat-dense regions.

Option Default Description
--alignment-backend pairhmm Phase 3 backend: pairhmm (WFA + PairHMM, default) or sw (Smith-Waterman only, no WFA triage). Invalid values are rejected at parse time.
--llr-threshold 2.3 PairHMM LLR threshold for confident calls (≈ ln(10))
--gap-open-prob 1e-4 PairHMM gap-open probability for non-repeat regions
--gap-extend-prob 0.1 PairHMM gap-extend probability for non-repeat regions
--repeat-gap-open-prob 1e-2 PairHMM gap-open probability for tandem repeat regions
--repeat-gap-extend-prob 0.5 PairHMM gap-extend probability for tandem repeat regions

pairhmm vs sw

pairhmm (default) uses WFA edit-distance triage first, then escalates to PairHMM only for ambiguous reads. This is faster than running PairHMM on every read and more accurate in low-quality or repeat-dense regions.

sw runs Smith-Waterman on every Phase 3 read (no WFA pre-filter). Use only if you need exact reproducibility with versions <3.0.0.

Read Judgment

How each read is judged against a row. The default is the size-specific rules, unchanged. The sequence judge is opt-in while it is measured toward becoming the default (#304; the default flip is #331).

Option Default Description
--read-judge rules rules: the size-specific read rules of the Read Judgment register. sequence: the sequence judge. Each read is scored whole against the row's candidate sequences (REF, the ALT and the other alleles at its site, each with every context of the other listed rows nearby and the recurring background), and counts REF or ALT only where its own bases favour one by at least one confident base; another allele winning by as much is partial evidence, and anything closer is depth only. haplotype, its former name, is taken as sequence with a warning. Invalid values are rejected at parse time. 
gbcms dna \
    --variants mutations.maf \
    --bam tumor:tumor.bam \
    --fasta reference.fa \
    --format maf \
    --read-judge sequence \
    --output-dir results/

Under --read-judge sequence:

  • Contexts explain, never credit. The other listed rows and the recurring background only explain a read's other differences: the judge does no linkage and no phasing, so a variant linked to the row elsewhere on the molecule never credits a read.
  • Four diagnostics can appear in gbcms_diagnostic: MNV_SHARED(chrom:pos:REF>ALT:n/m) at an SNV whose change a co-listed MNV holds (the MNV's carriers count at both rows, RJ-28), its indel twin INDEL_SHARED(chrom:pos:REF>ALT:n/m) at a row a listed indel or delins holds whole with listed neighbours (a delins that is two listed deletions together: its carriers count at each), UNEXPLAINED_READS(n/m), the row's depth reads that fit none of its candidate sequences (RJ-29), and NO_TELLING_BASE at an insertion of 50 bp or more whose equal windows no base within reach tells apart, every read depth only (RJ-23). See QC Flags. No column is added.
  • The Phase 3 options do not apply to the row's counts (--alignment-backend and the PairHMM probabilities above): the judge has its own scorer.
  • Spliced reads are not judged yet. A read whose CIGAR splices (an N op) is depth only, and gbcms rna refuses --read-judge sequence.
  • It is slower than the rules.
  • GBCMS_JUDGE_* environment variables are measurement switches for the judge's development, not options: leave them unset.

The design, the rule and every measurement round are in the sequence judge proposal.

Examples

gbcms dna \
    --variants mutations.vcf \
    --bam sample:sample.bam \ # (1)!
    --fasta reference.fa \
    --output-dir results/
  1. The sample: prefix sets the output filename. Without it, the BAM filename stem is used.
gbcms dna \
    --variants mutations.maf \
    --bam tumor:tumor.bam \
    --bam normal:normal.bam \
    --fasta reference.fa \
    --format maf \ # (1)!
    --output-dir results/
  1. MAF output preserves all input MAF columns and appends gbcms count columns.
gbcms dna \
    --variants mutations.vcf \
    --bam sample:sample.bam \
    --fasta reference.fa \
    --filter-duplicates \
    --min-mapq 30 \
    --output-dir results/

gbcms dna \
    --variants mutations.maf \
    --bam sample:sample.bam \
    --fasta reference.fa \
    --format maf \
    --show-normalization \ # (1)!
    --output-dir results/
  1. MAF output appends norm_Start_Position, norm_End_Position, norm_Reference_Allele and norm_Tumor_Seq_Allele2 (after any --column-prefix) showing the left-aligned variant; VCF output carries NORM_POS, NORM_REF and NORM_ALT INFO fields.
gbcms dna \
    --variants mutations.vcf \
    --bam sample:umi_deduped.bam \
    --fasta reference.fa \
    --umi-tag RX \ # (1)!
    --output-dir results/
  1. The RX tag is the standard SAM tag for UMI barcodes (fgbio, gencore).

Output

See Output Formats for a complete column-level schema reference covering:

  • VCF output: ##INFO fields, ##FORMAT fields, and annotated examples
  • MAF output: VCF→MAF vs MAF→MAF column sets, Tumor_Sample_Barcode behaviour, and column prefix options
  • mFSD columns (with --mfsd): all 41 mFSD fields
  • Normalization columns (with --show-normalization)
  • Companion files: .fsd.parquet (with --mfsd-parquet) and .observations.parquet (with --observations-parquet)

The run log ends with a QC summary: which QC flags fired, on how many rows, which rows need a look and what to do about each, then a WARNING when any row needs a look or an INFO line when none does. It reads the output after writing it and changes nothing. See QC Flags → End-of-run summary.