Skip to content

Troubleshooting

Common issues and solutions for gbcms. Issues are grouped by phase — work top-to-bottom for systematic diagnosis.

Installation Issues

Rust compilation fails: error: linker 'cc' not found

Install system build tools:

xcode-select --install

apt-get install build-essential

yum groupinstall "Development Tools"
htslib linking error: cannot find -lhts

The Rust backend compiles htslib statically. Ensure compression library headers are present:

apt-get install zlib1g-dev libbz2-dev liblzma-dev
pip builds from source, or GLIBC_2.28 not found

The one prebuilt wheel (any CPython from 3.10, Linux x86_64) needs glibc 2.28 or newer (RHEL 8 and newer) and pip 20.3 or newer. With an older glibc or an older pip, pip builds the source distribution instead (Platform support). On RHEL 8, install a Python from 3.10 on and use a virtual environment (python3.11 -m venv, then pip install --upgrade pip); see the Legacy Linux install path. For glibc older than 2.28, use that page's micromamba + clangdev source install.

LIBCLANG_PATH not set error during build

Required for bindgen (htslib headers). Set: 

export LIBCLANG_PATH=$CONDA_PREFIX/lib
Make sure you have clangdev (not clang) from conda-forge in your environment.

BAM / Input Issues

BAI index not found for sample.bam

Create or verify the index: 

samtools index sample.bam
# Expected: creates sample.bam.bai (or sample.bai)
samtools quickcheck sample.bam  # Sanity check the BAM
Chromosome X not found in reference / chromosome mismatch

gbcms looks up reference sequence by the chromosome name in the variant file. Ensure chromosome naming is consistent across VCF/MAF, BAM, and FASTA:

# Check names in FASTA
grep "^>" reference.fa | head -5
# Check names in BAM
samtools view -H sample.bam | grep "^@SQ" | head -5
# Check names in MAF/VCF
head -20 variants.maf | cut -f5   # Chromosome column

chr1 vs 1 and M/MT differences between the BAM and the variant file are auto-reconciled — you do not need to pre-harmonize them. A contig that has no match in the BAM is skipped with a one-time WARN (not a silent zero-count). The reference FASTA, however, is fetched by the variant's contig name, so the FASTA must contain that contig — otherwise the variant is rejected with verdict FAIL, reason FETCH_FAILED.

FETCH_FAILED reason (verdict FAIL) for all variants on a chromosome

This means the reference FASTA anchor could not be fetched for that contig. Because BAM↔variant contig naming is auto-reconciled, this is now almost always a FASTA problem, not a BAM one. Check: - The contig exists in the FASTA under the variant's name (grep "^>" reference.fa) - FASTA .fai index is present (samtools faidx reference.fa) - No truncated / corrupted FASTA

Allele Count Issues

Start with the QC summary at the end of the run log: it counts the rows by what they need (fix the input, check the call, or nothing), names the first rows to look at and says what to do about each flag (End-of-run summary).

alt_count = 0 for a variant I can see reads for in IGV

Check in this order:

  1. gbcms_status column — if FAIL, the variant was excluded from counting; gbcms_status_reason says why (each reason). See Normalization Issues.

  2. Complex variants count exact carriers only — a delins, or a deletion whose anchor also changes (e.g., GC→T, AG→T), counts a read as ALT only when its bases carry the whole given allele with two flank bases each side (exact-carrier rule). Reads carrying a different allele there are neither (partial_alt when closer to ALT), and gbcms_diagnostic names the allele they carry (OBSERVED_ALLELE / COEXISTING_ALLELE). Use --trace to see each read's call: 

    gbcms dna --trace --variants variants.maf --bam sample.bam \
    --fasta ref.fa --output-dir /tmp/debug/ 2>&1 | grep "read call"
    See Complex Indels for detailed case studies.

  3. Left-alignment shift — for large deletions (≥5bp), BWA left-aligns the anchor; the CIGAR D position may be several bases right of the anchor. Run normalization to see the left-aligned coordinates: 

    gbcms normalize --variants variants.maf --fasta ref.fa --output /tmp/norm/normalized.tsv
# Compare original_pos with norm_pos in the output TSV

  4. BAM coverage — verify the BAM actually has reads at this position: 

    samtools view -c sample.bam "CHROM:START-END"  # the event's span

  5. Verbose per-read trace — most diagnostic tool: 

    # Per-read trace lines print 0-BASED positions; variant-level lines are 1-based.
gbcms dna --trace ... 2>&1 | grep -E "<pos0>|<pos1>"  # the anchor, 0-based and 1-based
ref_count = 0 for a large deletion

A pure large deletion counts a clean read spanning its anchor as REF (check_deletion). A large delins (a deletion with a replacement) is judged by the exact-carrier rule: one pair of equal-length windows read from the POS flank, so a REF read counts when it carries the reference from that flank as far as the ALT read must carry the ALT.

If ref = 0, check: - The variant was normalized correctly (gbcms_status is PASS) - The BAM slice has coverage at the event's ends (samtools depth -a sample.bam -r CHROM:START-END) - Reads actually reach a junction: reads mapping entirely inside a large deleted span, or ending inside the event, carry no information about the variant and count as depth only

See NF2 Case Study.

alt_count is lower than expected for a deletion ≥5bp

BWA left-alignment can shift the anchor further left than where the CIGAR D appears. S3 checks the haplotype a placement gives, so a deletion written elsewhere in its shift region counts ALT; one written elsewhere with compensating mismatches counts ALT when the read's bases spell the ALT across the window.

For deletions <5bp failing S3 sequence validation, CIGAR-definitive REF is used — this is intentional. A 1-4bp deletion in the wrong reference context is almost certainly spurious noise, not a left-alignment artifact.

To diagnose deletion ≥5bp still returning low alt, read the per-read calls: each read call line names the rule that decided the read (rule=), and --verbose adds a Phase stats line per variant: 

gbcms dna --trace ... 2>&1 | grep "read call"
Large Δalt vs sign-out t_alt for MNP / ONP / DNP variants

Several expected causes — not gbcms bugs:

Cause Pattern Diagnosis
TRANS configuration Two MNP alleles on opposite chromosomes; sign-out counts 2× reads VAF ~2× expected
MAF annotation artifact ONP defined at wrong position or missing adjacent SNV gbcms normalize shows position shift
Sign-out duplication Variant captured by two MAF entries (DEL+INS pair both counted) Two nearby rows in MAF with same gene
Component carriers Carriers hold only part of the annotated haplotype (e.g. one of two G>A changes): alt_count ≈ 0, partial_alt ≈ sign-out t_alt, PARTIAL_DOMINANT Correct for the annotated allele. --rescue-mnp reports the best-supported component; gbcms_rescue shows the per-position split

Use --verbose to log per-variant read counts, or --trace for per-read detail.

dp = 0 for a variant
  • BAM slice doesn't cover this sample: samtools view -c sample.bam "chrom:pos-pos" returns 0
  • BAM index missing: see BAI issue above
  • Chromosome mismatch: see above
  • Variant outside BAM region: the BAM file only covers certain genomic regions
DPF is much lower than DP (large gap)

The gap DPF − (RDF + ADF) does not mean data loss — it counts fragments discarded due to R1/R2 disagreement within the quality threshold.

DPF - (RDF + ADF) = ambiguous fragments (both R1 and R2 covered but base quality within --fragment-qual-threshold)

A large gap at a specific locus is a quality signal indicating the site has many discordant fragments — typical of noisy or error-prone positions. Adjust with --fragment-qual-threshold (default 10).

ref_count + alt_count > total_count (INVARIANT violation)

This should never happen. If observed:

  1. Record the exact MAF row (chrom, pos, ref, alt, Tumor_Sample_Barcode)
  2. Run with --trace to capture per-read classification
  3. File a bug report at github.com/msk-access/gbcms/issues with the MAF row, BAM details, and trace log

The invariant ref_count + alt_count ≤ total_count is enforced throughout the counting engine.

Normalization Issues

REF_MISMATCH in gbcms_status — what does it mean?

The REF allele in the variant file does not match the reference FASTA at the stated position. gbcms requires ≥90% base match to auto-correct; below that, the variant is excluded.

Common causes: - Wrong reference genome build (GRCh37 data with GRCh38 reference or vice versa) - Chromosome naming mismatch (chr1 vs 1) - Trailing-base error in MAF annotation (fixed by tolerant validation if ≥90% match) - The REF written a few bases off (e.g. VCF alleles that keep the anchor base, with Start advanced as if it were trimmed): gbcms_diagnostic then says REF_AT_OFFSET(k), where the given REF matches exactly. The row is not moved or counted; correct the input's coordinates.

Diagnose with: 

gbcms normalize --variants variants.maf --fasta ref.fa --output /tmp/norm/normalized.tsv
grep "REF_MISMATCH" /tmp/norm/normalized.tsv | head -5   # gbcms_status_reason; REF_AT_OFFSET in gbcms_diagnostic
WARN_REF_CORRECTED reason (verdict PASS) — is this a problem?

The REF allele was ≥90% match — corrected to FASTA sequence and counting proceeded normally. The original MAF REF is preserved in the original_ref column (when --show-normalization is set).

This is not a problem in most cases. Common benign cause: 27bp EGFR exon 19 deletion with 26/27 bases correct and 1 trailing base that's a MAF annotation artifact.

If you want to audit all corrected variants: 

grep "WARN_REF_CORRECTED" output.maf | wc -l   # matches the gbcms_status_reason column
WARN_HOMOPOLYMER_DECOMP reason (verdict PASS) — what changed?

Only files from before 6.7.0 hold it: --rescue-homopolymer is deprecated in 6.7.0 (removed in 6.8.0) and ignored. The run used --rescue-homopolymer, and the variant is a homopolymer run called as a larger deletion (e.g. CCCCCC→T). In this sample, the corrected allele (the run with its last base replaced, CCCCCC→CCCCCT) got more ALT support, so its counts were used. used_decomposed = True in the output, and the flag is set per sample. Without the flag, the row counts the given allele, and OBSERVED_ALLELE in gbcms_diagnostic names the allele the reads carry.

It suggests the caller collapsed a smaller change at the run's end into one complex variant. The comparison is a heuristic: both alleles can claim reads carrying other forms of the run. Inspect the reads at the locus before relying on either allele's counts (see Variant Normalization).

After upgrading to 6.5.0, MAF rows from a VCF have different Start_Position / alleles

Since 6.5.0, VCF input is written to MAF exactly as vcf2maf writes it: the leading bases REF and ALT share are trimmed, so e.g. TCT>TCG becomes the SNP T>G at the changed base, and TTAC>A stays a 4bp deletion. The VCF record itself is kept in vcf_pos / vcf_ref / vcf_alt; key on those to compare runs. Counts are unchanged. Re-genotype every flavor with the same version before gbcms merge: MAFs from different versions do not join. gbcms convert shows the conversion for a variant file without counting. See Output Formats.

WARNING: Skipped VCF ALT ... — not countable

The VCF reader skips alleles that name no sequence to count: *, symbolic <...> alleles, breakends, a missing ., other non-sequence alleles, and a REF that is not a base sequence. The first few are logged individually, then one WARNING gives the totals by reason. They get no output row. See Input Formats.

RNA Mode Issues

All counts are 0 in RNA mode

  1. Strandedness filter with the wrong protocol — under the default --strandedness reverse (dUTP), a forward-stranded or unstranded library's sense reads read as antisense and are excluded. Set --strandedness forward or unstranded to match the library.

  2. MAPQ filter too strict — default MAPQ=1 with NH:i:1 rescue is correct for STAR. If your aligner assigns MAPQ differently, check --min-mapq.

  3. BAM not STAR-aligned — RNA mode assumes STAR BAMs with NH tag. Other aligners may lack the NH:i:1 tag needed for MAPQ rescue. Check: 

    samtools view sample.bam | head -1 | tr '\t' '\n' | grep "NH"
rna_antisense_depth is 0 and every read counts as sense

The variant has no gene strand. gbcms resolves it from the --gtf exons at the variant's position, or at an intronic position from the transcripts spanning it. Without one, every read passes as sense, strandedness is not enforced, and a run-level WARNING names such variants. Causes:

  • no --gtf, or a GTF that does not annotate the variant's contig;
  • an intergenic position (no transcript spans it);
  • genes of both strands over the position (overlapping antisense genes): both genes' transcripts carry the allele, so every read counts by design.

Solutions: - Pass --gtf with the annotation the BAM was aligned against. - For an unstranded library, use --strandedness unstranded.

RNA editing flag (rna_editing_site) is always False

Ensure --rna-editing-db points to a valid REDIportal TABLE1 file: 

# Verify file exists and is readable
head -3 TABLE1_hg38.txt.gz | zcat | head -3

# Download if missing
wget http://srv00.recas.ba.infn.it/atlas/download/TABLE1_hg38.txt.gz

Also verify chromosome naming matches (REDIportal uses chr-prefixed names).

Nextflow Issues

Container not found: Unable to find image 'ghcr.io/msk-access/gbcms:X.Y.Z'

Pull manually: 

docker pull ghcr.io/msk-access/gbcms:X.Y.Z
# or for Singularity
singularity pull docker://ghcr.io/msk-access/gbcms:X.Y.Z
Task metrics error: Command 'ps' required by nextflow

Use a newer container image with procps installed. This occurs with minimal base images that lack the ps command needed for Nextflow resource tracking.

Docker permission denied 
sudo usermod -aG docker $USER && newgrp docker
When does gbcms exit non-zero? (orchestration / Nextflow)

gbcms exits 1 if any sample fails (a BAM raises during counting — e.g. a missing index or a Rust panic surfaced as an error), whether all samples failed or only some. Successful runs exit 0.

An empty or fully-rejected variant set is not a failure: a sample that legitimately has no variants called (or whose variants are all rejected at prep) still exits 0, and the rejected variants are written to the output with gbcms_status FAIL and the reason in gbcms_status_reason. This lets a per-sample workflow keep going instead of failing the whole cohort on a sample that genuinely has nothing to count.

Debugging Tools

gbcms has two levels of diagnostic output:

Flag Output Use When
--verbose Per-variant summary: read counts, classification decisions First-pass diagnosis
--trace Per-read detail: CIGAR walk, S1/S2/S3 checks, Phase 2/3 scores Deep debugging a specific variant 
# Verbose: see per-variant decisions
gbcms dna --verbose --variants variants.maf --bam sample.bam \
    --fasta ref.fa --output-dir /tmp/debug/ 2>&1 | tee debug.log

# Trace: per-read detail (slow — run on small BAMs or specific regions).
# Per-read trace lines print 0-BASED positions; variant-level lines are 1-based,
# so grep both spellings of the coordinate.
gbcms dna --trace --variants variants.maf --bam sample.bam \
    --fasta ref.fa --output-dir /tmp/debug/ 2>&1 | grep -E "<pos0>|<pos1>" > trace.log

Inspecting Genomic Bin Construction

gbcms groups variants into ~10kb genomic bins before BAM traversal to reduce bam.fetch() calls (see Architecture → Genomic Binning). Bin construction is logged at INFO by default:

gbcms dna --variants variants.maf --bam sample.bam \
    --fasta ref.fa --output-dir /tmp/out/ 2>&1 | grep -i "bin"
# Example (INFO):
#   Built 42 genomic bins from 1247 variants (window=10000bp)
# Per-bin detail (read counts, filter breakdown) is DEBUG-level — add --verbose.

Expected bin counts

  • Targeted panels (MSK-ACCESS, IMPACT): many variants per gene → very few bins per chromosome — optimal I/O
  • Sparse WGS variants: ~1 variant per bin ≈ per-variant I/O (expected, not a bug)
  • Unexpectedly many bins for a clustered MAF? Check for chromosome name mismatches (chr1 vs 1) — the binner treats different chromosome names as separate bins regardless of coordinate proximity.