Bio Tumor Fraction Estimation

Version Compatibility

Reference examples tested with: ichorCNA 0.6.0+ (GavinHaLab fork), HMMcopy 1.40+, R 4.2+

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

• R: packageVersion('<pkg>') then ?function_name to verify parameters
• 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: ichorCNA is NOT an importable R function runIchorCNA() — it is a command-line script invoked as Rscript scripts/runIchorCNA.R with optparse flags, preceded by HMMcopy readCounter to build the WIG. Use the GavinHaLab fork (v0.6.0, 22 Nov 2024) for new work; the original broadinstitute/ichorCNA holds the wiki. The flag --repTimeWig does not exist — do not invent it.

Tumor Fraction Estimation

"Estimate tumor fraction from my cfDNA sample" -> Estimate the genome-wide proportion of cfDNA molecules that are tumor-derived, mutation-agnostic, from copy-number signal.

• CLI: readCounter (HMMcopy) to bin the BAM, then Rscript scripts/runIchorCNA.R for the HMM
• R: parse .params.txt (tumor fraction = 1 - n); a Python subprocess wrapper is a thin alternative

The Single Most Important Modern Insight -- tumor fraction is the quantity that travels across assays; it is NOT VAF and it is blind below ~3 percent

Tumor fraction (TF) is the fraction of cfDNA molecules that are tumor-derived — the cfDNA analogue of bulk-tumor purity. It is the burden metric that is comparable across assays and over time, which is exactly why it is the right unit to report. The two errors that dominate cfDNA work are unit confusion and floor confusion. Unit confusion: TF is NOT mutation VAF — a clonal heterozygous SNV in a diploid region sits at VAF approximately TF/2, so reporting a max-VAF as "tumor fraction" halves the true burden (and the factor changes entirely under LOH, amplification, or subclonality). Floor confusion: ichorCNA derives TF from copy-number deflection averaged over hundreds of 1 Mb bins, and that signal has a hard ~3 percent limit of detection. Below it the depth shift is smaller than per-bin sampling noise; a near-diploid or copy-neutral-LOH tumor returns a falsely low TF even at high true burden because it carries no depth signal. A low ichorCNA value is "low burden" only if the genome-wide plot is genuinely flat; otherwise it is uninformative, not negative.

Estimator Landscape

Decision Tree by Data Type

Methodology evolves; verify current best practice against the live ichorCNA wiki and the relevant tool docs before committing to an estimator.

ichorCNA Mechanics

ichorCNA is a hidden Markov model over copy-number states across 1 Mb bins. The emission per bin is the GC- and mappability-corrected log2 read-depth ratio (tumor vs a panel of normals). The HMM simultaneously segments the genome, calls large-scale CNAs (HOMD/DLOH/NEUT/GAIN/AMP/HLAMP up to maxCN), and by EM jointly estimates three global latent parameters: tumor fraction (via n), tumor ploidy (phi), and subclonal prevalence. The observed copy at a bin is a mixture: copy approximately 2*(1-TF) + TF*(tumor copy), and the sample ploidy identity is 2*(1-TF) + TF*tumor.ploidy. Because TF, ploidy, and per-bin tumor copy are all unknown, the same log-ratio can be explained by (low TF, large CN swing) or (high TF, small CN swing) — this ploidy/TF degeneracy is why ichorCNA fits over a grid of (normal, ploidy) start points and selects the maximum-likelihood solution.

Bin the BAM and Run the HMM

Goal: Produce a calibrated tumor-fraction estimate plus genome-wide CN segments from a single sWGS BAM.

Approach: Bin coverage into 1 Mb WIG with HMMcopy readCounter (chromosome naming must match the BAM @SQ style), then run runIchorCNA.R with build-matched GC/map/centromere references and a protocol-matched panel of normals; read .params.txt.

# Step 1: 1 Mb bins. --chromosome style ('1' vs 'chr1') MUST match the BAM @SQ names.
readCounter --window 1000000 --quality 20 \
  --chromosome "1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,X,Y" \
  tumor.bam > tumor.wig

# Step 2: the HMM. NOT an R function call -- it is a script with optparse flags.
Rscript scripts/runIchorCNA.R \
  --id tumor --WIG tumor.wig \
  --gcWig gc_hg38_1000kb.wig --mapWig map_hg38_1000kb.wig \
  --centromere GRCh38.centromere.txt \
  --normalPanel HD_ULP_PoN_1Mb_median.rds \
  --normal "c(0.5,0.6,0.7,0.8,0.9)" --ploidy "c(2,3)" --maxCN 7 \
  --estimateNormal TRUE --estimatePloidy TRUE --estimateScPrevalence TRUE \
  --scStates "c(1,3)" --txnE 0.9999999 --txnStrength 1e7 \
  --minMapScore 0.9 --genomeBuild hg38 --genomeStyle UCSC \
  --outDir ichor_out/

Key flags (verified defaults from runIchorCNA.R): --maxCN 7 (lower to 3 for low-TF); --normal "0.5" and --ploidy "2" are grid start points, not fixed values (--estimateNormal/--estimatePloidy still estimate them; these are optparse type=logical flags so they need an explicit TRUE/FALSE, not a bare flag); --txnE 0.9999999 and --txnStrength 1e7 set the segment-length prior; --minMapScore 0.9 drops low-mappability bins; --gcWig/--mapWig/--centromere/--normalPanel must all match the BAM's build and the 1 Mb bin size.

Parse the Optimal Solution

Goal: Extract the calibrated tumor fraction, ploidy, and QC from ichorCNA output.

Approach: Read .params.txt; TF = 1 - n_est for the selected (max-loglik) solution; gate on the GC-Map MAD; inspect subclonal fractions and the genome-wide plot before trusting a borderline call.

parse_ichor <- function(params_file) {
    p <- read.table(params_file, header = TRUE, sep = '\t', stringsAsFactors = FALSE)
    list(
        tumor_fraction = 1 - p$n_est[1],   # TF = 1 - normal fraction; selected solution is row 1
        ploidy = p$phi_est[1],
        loglik = p$loglik[1]
    )
}

The .params.txt also carries Tumor Fraction (= 1 - n), Tumor Ploidy (phi), Fraction Genome Subclonal, Fraction CNA Subclonal, and GC-Map Correction MAD (the data-noise QC). Companion outputs: .cna.seg (per-bin CN and log-ratio), .seg (IGV-compatible Viterbi segments), .RData (all grid solutions), and the genome-wide plot PDF — always inspect it for borderline calls because the ploidy/TF degeneracy can select a ploidy-3 alias of a ploidy-2 truth.