GPTomics/bio-ctdna-mutation-detection

Version Compatibility

Reference examples tested with: pysam 0.22+, pandas 2.2+, VarDictJava 1.8+, GATK 4.5+, Ensembl VEP 111+

Before using code patterns, verify installed versions match. If versions differ:

• Python: pip show <package> then help(module.function) to check signatures
• CLI: <tool> --version then <tool> --help to confirm flags

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

Notes specific to this skill: VarDict's -c -S -E -g are 1-based BED COLUMN INDICES, not genomic coordinates; var2vcf_valid.pl's -E suppresses the END tag (opposite meaning to VarDict's -E). VEP gnomAD flags are --af_gnomade (exomes)/--af_gnomadg (genomes); the bare --af_gnomad is a legacy alias that returns only exome AF, so prefer the explicit forms.

ctDNA Mutation Detection

"Detect mutations in my cfDNA sample" -> Either scan a panel for unknown low-VAF somatic variants (de novo calling) or quantify a pre-specified mutation set across samples (tumor-informed tracking) — two different statistical problems.

• CLI: vardict-java | teststrandbias.R | var2vcf_valid.pl for de novo low-VAF calling on a consensus BAM
• CLI: gatk Mutect2 with the read-orientation model for de novo calling with artifact filtering
• Python: pysam pileup of ref/alt counts at fixed loci for known-variant tracking (MRD)

The Single Most Important Modern Insight -- detection is a signal-vs-noise problem, and tracking is not the same problem as calling

Low-VAF ctDNA detection is set by two limits that no caller can overcome: the per-base error floor (raw Illumina ~1e-3 caps naive VAF detection near 0.5-1%) and the number of tumor molecules physically present in the tube (1 ng cfDNA ~= 303 haploid genome-equivalents; at 0.01% VAF in 10 ng the expected mutant count is ~0.3 copies — there is nothing to detect at any depth). The achieved limit of detection is the worse of the two. Error suppression (UMI consensus -> ~1e-5, duplex -> <1e-7) and input mass move the floor; swapping VarDict for Mutect2 does not.

Critically, de novo CALLING and known-variant DETECTION are different statistical problems. De novo calling scans every covered position for an unknown alt and pays a multiple-testing tax across 1e5-1e6 loci, so per-locus thresholds must be stringent (practical LoD ~0.1-0.5% on UMI consensus). Tumor-informed detection tests ONE hypothesis — "is tumor present?" — by integrating signal across a pre-specified set of N patient-specific loci; the multiple-testing penalty collapses and per-locus signal that is individually indistinguishable from noise sums into a confident panel-level call. This is why per-locus LoD is poor while panel-integrated LoD reaches single-ppm. Conflating the two is the most common conceptual error in the field.

Methods Landscape

| Method | Class | Citation | Role | When | |--------|-------|----------|------|------| | VarDict / vardict-java | de novo caller | AstraZeneca-NGS | sensitive low-VAF amplicon/capture calling with explicit strand-bias test | de novo panel calling on a UMI-consensus BAM | | Mutect2 (tumor-only) | de novo caller | GATK | local-assembly somatic caller + learned orientation-bias artifact model | de novo calling needing FFPE/OxoG artifact filtering, PoN, germline resource | | umi-varcal | UMI-aware caller | Sater 2020 Bioinformatics 36(9):2718 | own UMI-aware pileup + per-position Poisson test against local background | UMI-tagged BAM where a consensus-aware caller is wanted (floor ~0.3%) | | CAPP-Seq / iDES | tumor-informed integration | Newman 2014 Nat Med 20:548; 2016 Nat Biotechnol 34:547 | hybrid-capture deep panel + molecular barcoding + in-silico background polishing | de novo ctDNA to ~0.02%; iDES stacks ~15x error suppression | | INVAR | tumor-informed integration | Wan 2020 Sci Transl Med 12:eaaz8084 | integrate variant reads across 100s-1000s patient loci, background-weighted | MRD/monitoring with tumor WES; quantifies to ~1e-5, best ~2.5 ppm | | MRDetect | tumor-informed integration | Zviran 2020 Nat Med 26:1114 | shallow WGS vs patient SNV compendium, read-level SVM noise model | MRD trading depth for breadth (~35x WGS, thousands of SNVs); ~1e-5 | | PhasED-Seq | tumor-informed integration | Kurtz 2021 Nat Biotechnol 39:1537 | enrich phased (co-occurring) variants to suppress single-molecule error | sub-ppm MRD where phased variants are available |

Per-locus LoD for any of these is error- and sampling-limited (~0.1-0.5%); the tumor-informed methods reach ppm only by integrating across a known, large, patient-specific variant set. Methodology evolves — verify current best practice against each tool's live docs before committing to one.

Decision Tree by Scenario

| Scenario | Recommended | Why | |----------|-------------|-----| | Tumor tissue available, MRD/monitoring of a known cancer | Tumor-informed tracking (INVAR/MRDetect/Signatera-class), or the pysam tracker below for a fixed list | Integrating across N pre-specified loci is the only route to ppm; CHIP excluded by construction (CHIP variants are not on the tumor list) | | No tumor tissue, screening/discovery | de novo panel calling (VarDict or Mutect2) + matched WBC | Must scan for unknown variants; CHIP subtraction is mandatory or most calls are not tumor | | VAF regime > 1% | Any standard caller on a deduplicated BAM | Above the raw error floor; consensus not strictly required | | VAF regime 0.1-1% | UMI single-strand consensus + VarDict/umi-varcal/Mutect2 | Below the raw 1e-3 floor; consensus needed to recover real signal from error | | VAF regime < 0.1% | Duplex consensus + tumor-informed integration | Single-strand consensus cannot remove one-strand deamination/oxidation; only duplex + panel integration reaches this regime | | No matched WBC available | Do NOT report de novo calls as somatic-tumor | Without WBC subtraction, CHIP (the majority of non-germline cfDNA variants) is indistinguishable from tumor |

CHIP -- the dominant false positive, not background

Clonal hematopoiesis of indeterminate potential (CHIP) is the single largest source of false-positive somatic calls in plasma, and it is the null hypothesis for any low-VAF cfDNA variant. Razavi 2019 sequenced cfDNA with matched white-blood-cell DNA (508 genes, >60,000x) and found that 53.2% of non-germline cfDNA variants in cancer patients and 81.6% in non-cancer controls had features consistent with clonal hematopoiesis; only ~24.4% of cfDNA somatic variants in patients were also in the matched tumor (the remainder split between white-cell CHIP and variants of uncertain origin). These are bona fide somatic mutations — in cancer genes — that come from lysed leukocytes, not tumor. No error-suppression tier removes them because they are not errors.

The biology: CHIP arises in hematopoietic stem cells and rises steeply with age (Jaiswal 2014: ~10% prevalence over age 70), enriched for PPM1D/TP53/CHEK2 clones after prior chemo/radiation — exactly the monitored population. The canonical genes are DNMT3A, TET2, ASXL1 (the big three), then PPM1D, TP53, JAK2, SF3B1, SRSF2, GNB1, GNAS, CBL, ATM, CHEK2. TP53 and ATM are both CHIP genes and bona fide tumor suppressors, so a low-VAF TP53 cfDNA call is the ambiguous case par excellence.

The only reliable filter is matched buffy-coat/WBC subtraction: sequence the WBC fraction of the same draw at comparable depth and remove any cfDNA variant also present in WBC. gnomAD filtering removes germline only — CHIP variants are somatic and absent from germline databases, so they sail straight through. A canonical-CHIP-gene list (the example's CHIP_GENES) is a heuristic flag for extra scrutiny, NOT a substitute for WBC subtraction. See analytical-validation for the LoB/LoD statistics that quantify how confidently a subtracted call clears background.

De Novo Calling with VarDict

Goal: Scan a target panel for unknown low-VAF somatic variants on a UMI-consensus BAM.

Approach: Run vardict-java with a lowered -f, pipe through the strand-bias test, then convert to VCF — matching -f across both stages so the threshold is not silently re-applied.

AF_THR=0.005   # 0.5% — practical UMI-consensus de novo floor; below this approaches the per-base error floor
vardict-java -G ref.fa -f $AF_THR -N sample -b consensus.bam \
  -c 1 -S 2 -E 3 -g 4 targets.bed | \
  teststrandbias.R | \
  var2vcf_valid.pl -N sample -E -f $AF_THR > sample.vcf

Key flags: -G indexed reference; -f min VAF (VarDict default 0.01); -N sample name; -b BAM. -c 1 -S 2 -E 3 -g 4 are the 1-based BED COLUMN INDICES for chrom/start/end/gene in a standard 4-column BED — they are column positions, not genomic values. On var2vcf_valid.pl, -E means "do NOT print the END tag" (unrelated to VarDict's -E); its -f default is 0.02, so set it to match. For PCR/amplicon data add -P 0 (positional std is expected to be ~0). For paired tumor/normal use the testsomatic.R | var2vcf_paired.pl path instead.

De Novo Calling with Mutect2 and the Orientation-Bias Model

Goal: Call de novo somatic variants while filtering FFPE-deamination (C>T) and OxoG (G>T) strand-biased artifacts that dominate low-VAF false positives.

Approach: Collect F1R2/F2R1 counts during calling, learn the orientation-bias prior, then apply it during filtering alongside a panel of normals and germline resource.

gatk Mutect2 -R ref.fa -I consensus.bam --f1r2-tar-gz f1r2.tar.gz \
  --germline-resource af-only-gnomad.vcf.gz --panel-of-normals pon.vcf.gz \
  -O unfiltered.vcf.gz
gatk LearnReadOrientationModel -I f1r2.tar.gz -O read-orientation-model.tar.gz
gatk FilterMutectCalls -R ref.fa -V unfiltered.vcf.gz \
  --ob-priors read-orientation-model.tar.gz -O filtered.vcf.gz

Mutect2 is run tumor-only here (no normal sample arg); the orientation model is the load-bearing low-VAF filter. At true ctDNA VAFs Mutect2 is underpowered relative to a dedicated UMI/duplex + background-polishing pipeline, and local assembly can miss extremely low-AF alt support — it is a reasonable de novo caller on consensus reads with the orientation model + PoN (and ideally a matched normal, not shown in this tumor-only command), not a substitute for tumor-informed integration at ppm.

Track Known Mutations Across Serial Samples

Goal: Quantify the VAF of a pre-specified mutation set at fixed loci for MRD monitoring — the detection (not calling) problem.

Approach: For each target mutation, pileup reads at the position, count ref/alt/other alleles, and compute VAF with depth; aggregate across loci as the panel-level detection signal. The single-base pileup below tracks SNVs only — indel reporters (e.g. EGFR exon-19 deletions) need read.indel/CIGAR-aware counting; a single-base comparison silently scores every indel read as other and reports the locus as cleared.

import pysam

def track_known_variants(bam_file, variants):
    '''Pileup ref/alt counts at fixed (chrom, pos, ref, alt) SNV loci; pos is 1-based.
    SNVs only - indel reporters need read.indel/CIGAR handling, not a single-base compare.'''
    bam = pysam.AlignmentFile(bam_file, 'rb')
    rows = []
    for chrom, pos, ref, alt in variants:
        counts = {'ref': 0, 'alt': 0, 'other': 0}
        for col in bam.pileup(chrom, pos - 1, pos, truncate=True):
            for read in col.pileups:
                if read.is_del or read.is_refskip:
                    continue
                base = read.alignment.query_sequence[read.query_position]
                counts['alt' if base == alt else 'ref' if base == ref else 'other'] += 1
        depth = sum(counts.values())
        rows.append({'chrom': chrom, 'pos': pos, 'ref': ref, 'alt': alt,
                     'depth': depth, 'alt_count': counts['alt'],
                     'vaf': counts['alt'] / depth if depth else 0.0})
    bam.close()
    return rows

Annotate calls for interpretation with Ensembl VEP (--cache --offline --fasta --vcf --everything); the gnomAD allele-frequency flags are --af_gnomade (exomes) and --af_gnomadg (genomes) — the bare --af_gnomad is a legacy alias returning only exome AF. gnomAD presence separates germline; only WBC presence separates CHIP.

Per-Method Failure Modes

CHIP misclassified as tumor

Trigger: de novo calling without matched WBC. Mechanism: leukocyte-derived clonal somatic variants in cancer genes look identical to tumor signal and pass gnomAD filtering. Symptom: low-VAF calls in DNMT3A/TET2/TP53; "tumor" mutations not in the matched tissue. Fix: subtract matched buffy-coat/WBC genotype; never report somatic-tumor without it.

Strand-biased / deamination artifacts at low VAF

Trigger: FFPE-style C>T or oxidative G>T at VAF near the floor. Mechanism: damage on one template strand is inherited by every PCR copy, so single-strand UMI consensus votes unanimously for the artifact. Symptom: alt support concentrated on one strand. Fix: VarDict strand-bias test or Mutect2 orientation model; for sub-0.1% require duplex consensus.

Calling below the error floor

Trigger: lowering -f to e.g. 0.001 on a non-consensus BAM. Mechanism: the raw ~1e-3 error rate manufactures alt reads at that frequency. Symptom: a flood of low-VAF calls scaling with depth. Fix: do consensus upstream; do not set a VAF threshold below the demonstrated error floor of the input.

Germline-vs-somatic confusion at low coverage

Trigger: classifying by VAF alone when depth is low. Mechanism: a true 50% het reads 3/12 = 0.25 by chance. Symptom: germline hets mislabeled subclonal somatic. Fix: gnomAD + matched-WBC presence (germline ~0.5 in WBC; CHIP at clone VAF; tumor-only absent from WBC).

Per-locus LoD quoted as the assay LoD

Trigger: reporting a single-variant sensitivity for a multi-locus tracking assay (or vice versa). Mechanism: panel-integrated LoD is orders of magnitude below per-locus LoD. Symptom: a "0.1%" claim that does not match observed ppm-level tracking. Fix: state per-locus vs panel-integrated explicitly; see analytical-validation.

Quantitative Thresholds

| Threshold | Source | Rationale | |-----------|--------|-----------| | Raw Illumina error ~1e-3 caps naive VAF near 0.5-1% | Schmitt 2012 PNAS 109:14508 | Per-base miscall rate sets the per-locus VAF floor; alt support below it is mostly error | | UMI single-strand consensus -> ~1e-5; duplex -> <1e-7 | Schmitt 2012; Newman 2016 Nat Biotechnol 34:547 | Family consensus erases PCR/sequencing error; duplex strand concordance also catches one-strand damage | | CAPP-Seq de novo LoD ~0.02% at 96% specificity | Newman 2014 Nat Med 20:548 | Deep hybrid-capture + reporter set; demonstrates the de novo panel floor | | iDES ~15x error suppression (UMI ~3x x polishing ~3x) | Newman 2016 Nat Biotechnol 34:547 | Molecular consensus and in-silico background polishing are orthogonal and stack | | INVAR quantifies to ~1e-5, detects to ~2.5 ppm | Wan 2020 Sci Transl Med 12:eaaz8084 | Integrating variant reads across 100s-1000s of patient loci collapses multiple testing | | MRDetect ~1e-5 tumor fraction at ~35x WGS, 95% spec | Zviran 2020 Nat Med 26:1114 | Breadth (thousands of SNVs) + read-level SVM (~14.4x error reduction) substitutes for depth | | CHIP = 53.2% (cancer pts) / 81.6% (controls) of cfDNA variants | Razavi 2019 Nat Med 25:1928 | Most non-germline cfDNA variants are not tumor; matched WBC is mandatory | | Depth >= 1000-5000x unique consensus for panels | community / Phallen 2017 Sci Transl Med 9:eaan2415 (~30,000x) | Detecting <1% VAF needs enough unique molecules sampled at each locus | | ~303 genome-equivalents per ng cfDNA (3.3 pg/haploid) | standard constant | Input mass sets a hard Poisson ceiling on detectable VAF independent of sequencing | | LoB / LoD / LoD95 per CLSI EP17 | CLSI EP17-A2 | A bare VAF without input mass + replicate detection rate is not a sensitivity spec |

Common Errors

| Error / symptom | Cause | Solution | |-----------------|-------|----------| | Flood of low-VAF calls scaling with depth | -f set below the input's error floor on non-consensus reads | Do UMI/duplex consensus first; keep -f >= demonstrated floor | | VarDict emits nothing or wrong regions | -c -S -E -g read as genomic values | They are 1-based BED column indices; use -c 1 -S 2 -E 3 -g 4 for a 4-column BED | | var2vcf re-filters away VarDict calls | var2vcf_valid.pl -f default 0.02 mismatched | Set var2vcf -f to match VarDict's -f | | Real amplicon calls dropped as positional artifacts | var2vcf -P (filter pstd=0) on by default | Add -P 0 for PCR/amplicon data | | "Tumor" variants absent from matched tissue | CHIP not subtracted | Sequence and subtract matched WBC; flag CHIP-gene hits | | --af_gnomad returns only exome AF | bare flag is a legacy exome-only alias | Use --af_gnomade (exomes) / --af_gnomadg (genomes) | | ppm "LoD" not reproducible | per-locus LoD quoted for a tracking assay | Report panel-integrated LoD with input mass and LoD95 |

References

• Schmitt MW, Kennedy SR, Salk JJ, Fox EJ, Hiatt JB, Loeb LA. 2012. Detection of ultra-rare mutations by next-generation sequencing. Proc Natl Acad Sci USA 109(36):14508-14513. — Duplex Sequencing; DCS error <1e-7.
• Newman AM, Bratman SV, To J, et al. 2014. An ultrasensitive method for quantitating circulating tumor DNA with broad patient coverage. Nat Med 20(5):548-554. — CAPP-Seq; ~0.02% LoD.
• Newman AM, Lovejoy AF, Klass DM, et al. 2016. Integrated digital error suppression for improved detection of circulating tumor DNA. Nat Biotechnol 34(5):547-555. — iDES; UMI + background polishing.
• Phallen J, Sausen M, Adleff V, et al. 2017. Direct detection of early-stage cancers using circulating tumor DNA. Sci Transl Med 9(403):eaan2415. — TEC-Seq deep panel.
• Razavi P, Li BT, Brown DN, et al. 2019. High-intensity sequencing reveals the sources of plasma circulating cell-free DNA variants. Nat Med 25(12):1928-1937. — CHIP is the majority of cfDNA variants; matched WBC.
• Wan JCM, Heider K, Gale D, et al. 2020. ctDNA monitoring using patient-specific sequencing and integration of variant reads. Sci Transl Med 12(548):eaaz8084. — INVAR; integration to ~2.5 ppm.
• Zviran A, Schulman RC, Shah M, et al. 2020. Genome-wide cell-free DNA mutational integration enables ultra-sensitive cancer monitoring. Nat Med 26(7):1114-1124. — MRDetect; shallow WGS + read SVM.
• Kurtz DM, Soo J, Co Ting Keh L, et al. 2021. Enhanced detection of minimal residual disease by targeted sequencing of phased variants in circulating tumor DNA. Nat Biotechnol 39(12):1537-1547. — PhasED-Seq; phased-variant enrichment.
• Sater V, Viailly P-J, Lecroq T, et al. 2020. UMI-VarCal: a new UMI-based variant caller that efficiently improves low-frequency variant detection in paired-end sequencing NGS libraries. Bioinformatics 36(9):2718-2724. — UMI-aware pileup + per-position Poisson test.

Related Skills

• cfdna-preprocessing - UMI/duplex consensus input that sets the error floor
• analytical-validation - LoD/LoB and the panel-integration math behind detection
• longitudinal-monitoring - track detected variants across serial samples
• tumor-fraction-estimation - orthogonal burden estimate to cross-check
• variant-calling/variant-calling - general somatic calling principles
• clinical-databases/variant-prioritization - clinical annotation and interpretation