Skip to content

gbcms merge

Merge per-BAM-type genotyped MAFs (e.g., duplex, simplex) into a single output MAF with type-prefixed count columns.

Polars Batch Engine

This command uses Polars for lazy-evaluated, memory-efficient joins and aggregation. It operates on the complete file (batch mode), unlike the streaming per-variant counting in gbcms dna/rna.


Quick Start

# Merge duplex + simplex MAFs with combined columns
gbcms merge \
    --input duplex:sample1-duplex.maf \
    --input simplex:sample1-simplex.maf \
    --output sample1.merged.maf

# Merge 3 BAM types, no combined columns
gbcms merge \
    --input duplex:sample1-duplex.maf \
    --input simplex:sample1-simplex.maf \
    --input standard:sample1-unfiltered.maf \
    --output sample1.merged.maf \
    --no-combined

Arguments

Required

Flag Type Description
--input TYPE:PATH BAM type label and path to its genotyped MAF. Repeat for each input (minimum 2). Labels become column prefixes.
--output PATH Output path for the merged MAF.

Optional

Flag Default Description
--no-combined false Skip computing simplex_duplex_* combined columns.
--legacy-naming false Use t_{metric}_{type} column naming for genotype_variants compatibility.

How It Works

flowchart LR
    subgraph Inputs
        D["duplex.maf"]
        S["simplex.maf"]
    end

    subgraph Engine ["Merge Engine (Polars)"]
        Scan["Lazy Scan"] --> Prefix["Auto-prefix\ncount columns"]
        Prefix --> Join["Outer Join\non variant key"]
        Join --> Fill["Null → 0\n(missing counts)"]
        Fill --> Combine["Additive\ncombined columns"]
        Combine --> SB["Fisher Strand Bias\n(Rust)"]
    end

    subgraph Output
        M["merged.maf"]
    end

    D --> Scan
    S --> Scan
    SB --> M
Use mouse to pan and zoom

Step-by-Step

  1. Scan — Each input MAF is lazily scanned with all columns as strings
  2. Per-input names — Every gbcms column (counts, status, diagnostics, strand bias, mFSD, RNA; the set is the writer's) is renamed to {type}_{name} (e.g., ref_count → duplex_ref_count). Counts written with --column-prefix (duplex_ as the pipeline runs it, t_) are found under that prefix; the other columns are written unprefixed and found as they are.
  3. Outer Join — Progressive full outer join on the variant: Chromosome, Start_Position, Reference_Allele, Tumor_Seq_Allele2. End_Position follows from Start and REF and is not joined on, so inputs that write it differently for one variant still join (the output keeps the first input's, and the log names each input that differs); a variant only a later input has takes its End_Position from that input. When an input lists one variant twice with different End_Position, it joins on End_Position too, so each row pairs with its own counterpart. The output's columns are the first input's: without End_Position there, the output has none. When every input is VCF-derived (carries vcf_pos / vcf_ref / vcf_alt), the VCF record is added to the key: two VCF records can describe one MAF record (TCT>TCG and T>G at the changed base), and each should stay its own row. The MAF columns stay in the key, so a variant only a later input has keeps its coordinates and alleles. An input that lists one join key more than once gets a WARN, since each such row pairs with every matching row of the other inputs. Merge MAFs genotyped by the same gbcms version: merge warns when the inputs' version lines differ, and stops when one input is a VCF-input MAF from before 6.5.0 (vcf_pos without vcf_ref/vcf_alt) and another is not, since 6.5.0 changed VCF-input coordinates and their rows would not join. A row only a later input has takes the first input's annotation columns from the earliest later input that has it. Chromosome is compared by the same rule counting uses to reconcile contigs (chr1 ~ 1, chrM ~ MT), so inputs counted from differently named variant files still join. Each contig is written one way. That is the first input's name for it, or, for a contig the first input lacks, the name used by the earliest input that has it. When a later input names contigs differently, the log says so (one INFO line per input). An input column named _contig_key, or starting with _row_, _chrom_, _end_ or _ann_, is rejected: those names are reserved for the join. If one input names a contig two ways (e.g. chrM and MT), a variant listed under both joins once per name, and merge logs a WARN.
  4. Absent rows — A row an input lacks gets "0" for that input's read and fragment counts, "" for its status columns, and empty mFSD/RNA columns; the log counts the rows each input lacks. A row an input has keeps its cells as written: a missing count stays missing.
  5. Combined Columns — If both simplex and duplex are present and --no-combined is not set:
    • Phase 1: Additive sums (12 columns: read + fragment + strand counts); NA where either flavor's cell is missing or non-numeric, warned once per column
    • Phase 2a: Derived totals (total_count, total_count_fragment)
    • Phase 2b: Derived VAFs (vaf, vaf_fragment)
    • Phase 3: Fisher's exact test for strand bias (read + fragment level, via Rust)
  6. Write — Materialized DataFrame written as tab-separated MAF, after provenance lines (its own #gbcms/#command, and one #input line per input with that input's version line), atomically (a temp file renamed into place). Rows follow the inputs: the first input's rows in its order, then rows only a later input has, in that input's order. The output is identical on every run.

Provenance Comment Lines (v5.3.0+)

Starting in v5.3.0, gbcms MAF output includes #gbcms and #command comment lines before the TSV header. The merge engine handles these automatically — Polars scan_csv is configured with comment_prefix="#", which skips all #-prefixed lines. No manual pre-processing is needed.


Combined Columns

When both simplex and duplex inputs are present, 20 combined simplex_duplex_* columns are computed (assuming all strand-level counts are present in the input):

Phase Column Formula
Sum simplex_duplex_ref_count simplex_ref_count + duplex_ref_count
Sum simplex_duplex_alt_count simplex_alt_count + duplex_alt_count
Sum simplex_duplex_ref_count_forward Additive
Sum simplex_duplex_ref_count_reverse Additive
Sum simplex_duplex_alt_count_forward Additive
Sum simplex_duplex_alt_count_reverse Additive
Sum simplex_duplex_ref_count_fragment Additive
Sum simplex_duplex_alt_count_fragment Additive
Sum simplex_duplex_ref_count_fragment_forward Additive
Sum simplex_duplex_ref_count_fragment_reverse Additive
Sum simplex_duplex_alt_count_fragment_forward Additive
Sum simplex_duplex_alt_count_fragment_reverse Additive
Derived simplex_duplex_total_count ref + alt
Derived simplex_duplex_total_count_fragment ref_fragment + alt_fragment
Derived simplex_duplex_vaf alt / total (0 when total=0)
Derived simplex_duplex_vaf_fragment alt_frag / total_frag
Fisher simplex_duplex_strand_bias_p_value Fisher 2×2 (read strand)
Fisher simplex_duplex_strand_bias_odds_ratio Fisher 2×2 (read strand)
Fisher simplex_duplex_fragment_strand_bias_p_value Fisher 2×2 (fragment strand)
Fisher simplex_duplex_fragment_strand_bias_odds_ratio Fisher 2×2 (fragment strand)

Schema-Aware

If the input MAFs do not contain strand-level columns (e.g., from an older gbcms version), only the available columns are summed. Missing metrics are logged at INFO level and skipped — the pipeline does not fail.

Mixed MNP rescue between duplex and simplex

With inputs genotyped using --rescue-mnp, a row can be rescued in one flavor (its counts are then a component SNV's) but not the other (the MNP's own counts), or rescued to different components. The two flavors' columns then describe different alleles in one row, so the row's combined simplex_duplex_* cells are NA (summing them would add two alleles). gbcms merge names every such row in a WARNING (with or without --add-combined) and logs a per-run count; the per-flavor counts are left unchanged. A row one flavor lacks is not mixed. If only one flavor was genotyped with --rescue-mnp, merge logs that once and treats every row rescued in that flavor as mixed, since the other flavor reports the MNP there. Check duplex_gbcms_rescue / simplex_gbcms_rescue for those rows.

Validated on real ACCESS data (2026-09-23)

The additive combination rests on duplex and simplex being disjoint molecule sets. Checked on real clinical data first: zero QNAME overlap between the duplex and simplex BAMs at three loci. Then a flavor-pair truth harness (#92 comment): 4 samples and 256 signed-out variants, each sample genotyped twice (duplex BAM and simplex BAM separately, default filters, one streamed BAM at a time), combined with gbcms merge --add-combined, and scored as simplex_duplex_*_fragment against sign-out.

  • With t_alt of 5 or more (n=252): median relative error +0.0% (IQR +0.0% to +0.6%); by type SNV +0.0% (n=191), DEL +2.9% (n=44), INS +7.3% (n=17). No missed variants at any t_alt.
  • The only deviations over 25% were five deletions in one hypermutation cluster of about 150bp in BRCA2, each over-counted 2-5× before the cluster contest of 6.5.0 (three 2bp deletions, a 33bp deletion, a 14bp deletion and four 1-3bp inserts interleaved in an AG-rich tract); see Read Judgment, evidence behind RJ-5.
  • This was also the first real-data exercise of the combined-sums path.

Nextflow Integration

When running via the Nextflow pipeline, enable merge with:

nextflow run main.nf \
    --input samplesheet.csv \
    --variants mutations.maf \
    --fasta reference.fa \
    --format maf \
    --merge_counts \
    -profile slurm

The samplesheet should include a bam_type column:

sample,bam,bai,suffix,bam_type
sample1,sample1.duplex.bam,,-duplex,duplex
sample1,sample1.simplex.bam,,-simplex,simplex

When bam_type is set, the DNA module automatically derives --column-prefix from the type label, so counts are pre-prefixed in the per-BAM MAF output.