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

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

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.