本文へ移動
cccskills
無料GitHub で公開

bio-tumor-fraction-estimation

Estimates tumor fraction (the genome-wide proportion of cfDNA molecules that are tumor-derived, the cfDNA analogue of bulk-tumor purity) from shallow whole-genome sequencing with ichorCNA, an HMM over 1 Mb bins that jointly EM-estimates tumor fraction, ploidy, and subclonal prevalence over a normal/ploidy grid. Encodes the load-bearing reframes: tumor fraction is the quantity that travels across assays and is NOT mutation VAF (clonal-het VAF approximately TF/2), CNA-based estimation has a hard ~3 percent limit-of-detection floor, and near-diploid or copy-neutral-LOH genomes return a falsely low value. Selects the estimator by data type (sWGS to ichorCNA, deep panel to max-VAF, methylation to deconvolution, sub-3 percent to fragmentomics or methylation). Use when quantifying tumor burden from a liquid biopsy, picking a tumor-fraction estimator for a given assay, or reconciling a TF estimate against a panel VAF.

インストール方法を見る

含まれるファイル(3)

  • SKILL.md14.9 KB
  • examples/estimate_tumor_fraction.R3.5 KB
  • usage-guide.md3.6 KB

SKILL.md(原文)

インストールする前に、エージェントに与えられる指示の中身を確認できます。

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

EstimatorClassInputStrengthFails when
ichorCNACNA / depth (HMM)sWGS 0.1-1xMutation-agnostic genome-wide burden; calibrated standardTF < ~3%; near-diploid / copy-neutral-LOH genome
TitanCNACNA + allelic (B-allele)deeper WGS with het-SNP depthResolves CNLOH via allelic imbalanceNeeds informative het-SNP coverage (not 0.1x)
max-VAF / clonal-cluster MAFMutation / paneldeep targeted or WESSensitive to <0.1% VAF with UMI/duplexNeeds callable variants + CHIP filtering; CN-sensitive
Methylation deconvolution (CelFiE, CelFEER)MethylationWGBS/EM-seq or methyl panelDense per-molecule signal reaches below CNA floorNeeds a tumor-type methylation reference atlas
Fragmentomics (Griffin, DELFI)FragmentomicsWGSCN-independent corroboration at low TFQuantifies "tumor signal," not a calibrated molecular fraction

Decision Tree by Data Type

Data availableTF regimeRecommendedWhy
sWGS 0.1-1x, no known variants, aneuploid tumor>= ~3%ichorCNAMutation-agnostic genome-wide burden; the standard
sWGS, tumor type known to be near-diploid / quietanymutation or methylationichorCNA underestimates with no depth signal
sWGS, TF suspected < 3%< 3%deep-panel max-VAF, methylation, or fragmentomicsBelow the CNA floor (see fragment-analysis, methylation-based-detection)
Deep targeted / WES paneldown to <0.1% VAFmax-VAF excl. CHIP, or clonal-cluster MAFPer-locus sensitivity; convert via TF approximately 2*VAF with CN care (see ctdna-mutation-detection)
Methylation (WGBS/EM-seq/panel)very lowmethylation deconvolutionDense per-molecule signal; needs reference atlas
Targeted panel, want CN-based TF>= few %ichorCNA on off-target readsRecovers genome-wide CN from off-target coverage

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.

Unit Confusion -- TF vs VAF vs ctDNA%

These three are routinely conflated; the relation is exact and copy-number-dependent. For a variant at local copy number Cn with mutant-copy multiplicity m:

VAF = (TF * m) / [ TF * Cn + 2 * (1 - TF) ]

For a clonal heterozygous SNV in a diploid region (Cn=2, m=1) this collapses to VAF approximately TF/2, equivalently TF approximately 2*VAF. The common errors:

  • LOH variant (mutant on both copies, normal copy lost): m=Cn, so VAF -> TF, not TF/2 — treating it as TF/2 doubles the estimate.
  • Amplified mutant allele inflates VAF above TF/2; a mutant on a deleted copy deflates it — so max-VAF over-estimates TF for amplified drivers and under-estimates for deleted ones.
  • Subclonal variants carry an extra cancer-cell-fraction factor and understate TF.
  • Germline heterozygous SNPs sit at VAF approximately 0.5 regardless of TF — never feed them into a TF-from-VAF calculation.
  • CHIP (clonal hematopoiesis) variants are blood-derived, not tumor — exclude them from any max-VAF TF proxy.
  • ctDNA% is loosely used for either TF or max-VAF; always pin down which a lab means.

Cross-check: for a clonal heterozygous driver in a diploid region, ichorCNA TF and 2*(panel VAF) should agree. TF >> 2VAF implies a subclonal/deleted variant or a ploidy mis-call; TF << 2VAF implies a near-diploid/CNLOH tumor or an amplified/LOH driver. Never average the two blindly (see ctdna-mutation-detection).

Per-Method Failure Modes

~3 percent CNA floor

Trigger: TF below 0.03 at 0.1x sWGS. Mechanism: the log2 deflection from a single-copy event is proportional to TF (±0.02 at TF=0.03), smaller than per-bin sampling noise; only averaging over hundreds of bins recovers it. Symptom: TF collapses toward 0; replicate variability (MNSD) rises sharply. Fix: the floor scales with aneuploidy magnitude and coverage — it needs roughly one >100 Mb gain AND one >100 Mb loss; sequence deeper (>1-5x) or switch estimator class (fragment-analysis, methylation-based-detection).

Near-diploid / copy-neutral-LOH

Trigger: quiet tumor type or CNLOH-rich genome. Mechanism: CNLOH has identical total coverage to diploid, indistinguishable on depth alone; ichorCNA is also tuned conservative and "may underestimate." Symptom: falsely low TF with a flat genome-wide plot. Fix: treat a flat low call as uninformative, not negative; escalate to a mutation/methylation assay; TitanCNA can use allelic imbalance if het-SNP depth exists.

Mismatched panel of normals / references

Trigger: PoN, GC/map/centromere WIG, or build does not match the library prep, bin size, or genome build. Mechanism: the PoN models protocol-specific coverage bias; a mismatched PoN injects its own bias as spurious CN waviness. Symptom: wavy log-ratio, implausible TF. Fix: build/obtain a PoN from healthy-donor cfDNA on the exact protocol at the same bin size and build; keep hg19 vs hg38 and 1 vs chr1 consistent end-to-end.

Ploidy aliasing

Trigger: ploidy/TF degeneracy. Mechanism: the max-loglik solution is occasionally a ploidy-3 alias of a ploidy-2 truth. Symptom: doubled ploidy with halved TF. Fix: read all .params.txt solutions, inspect the plot; for low-TF samples force --ploidy "c(2)".

Quantitative Thresholds

ThresholdSourceRationale
Coverage 0.1-1x sWGS; 1 Mb binsAdalsteinsson 2017; ichorCNA wikiFiner bins add noise at 0.1x; ~0.1x is the calibrated ULP-WGS operating point
~3% TF limit of detection at ~0.1xAdalsteinsson 2017 (95% sens, 1125/1288 mixtures; 91% spec, 20/22 donors at 0.03 TF cutoff)Below 0.03 the depth deflection falls under per-bin noise
97.2-100% sensitivity to detect 3% TF (1x and 0.1x)J Mol Diagn 2024 assay validationIndependent dilution/replicate validation; MNSD rises sharply below 3%, establishing 3% as the LOD
GC-Map Correction MAD < 0.15 good; > 0.3 distrustichorCNA FAQResidual post-correction noise; high MAD means the depth signal is unreliable
Manual-curation band 0.03-0.10 TFichorCNA wikiModel can pick the wrong solution and tends to underestimate near the floor; inspect the plot
Low-TF recipe: --normal "c(0.95,0.99,0.995,0.999)" --ploidy "c(2)" --maxCN 3 --estimateScPrevalence FALSE --scStates "c()"ichorCNA wikiSeeds EM near TF 5/1/0.5/0.1%; ploidy and subclonality are unidentifiable when CN signal is weak

References

  • Adalsteinsson VA, Ha G, Freeman SS, et al. 2017. Scalable whole-exome sequencing of cell-free DNA reveals high concordance with metastatic tumors. Nat Commun 8(1):1324. — ichorCNA primary method; the ~3% LOD benchmark (95% sensitivity, 91% specificity at a 0.03 TF cutoff, ~0.1x).
  • Assay Validation of Cell-Free DNA Shallow Whole-Genome Sequencing to Determine Tumor Fraction in Advanced Cancers. 2024. J Mol Diagn 26(5):413-422 (PMC11090203). — Independent validation: 97.2-100% sensitivity at 3% TF (1x and 0.1x); MNSD rising below 3% establishes 3% as the LOD.
  • broadinstitute/ichorCNA and GavinHaLab/ichorCNA GitHub repositories and wiki (Usage, Output, Parameter-tuning, Create-Panel-of-Normals, FAQ). — readCounter command, runIchorCNA.R flag defaults, .params.txt fields, MAD QC thresholds, CNLOH/near-diploid underestimation, PoN construction.

Related Skills

  • cfdna-preprocessing - sWGS BAM input and minimal-processing path
  • fragment-analysis - the estimator to use below the ~3% CNA floor
  • ctdna-mutation-detection - max-VAF cross-check and the TF-vs-VAF reconciliation
  • analytical-validation - the ~3% floor framed as a limit of detection
  • copy-number/cnvkit-analysis - copy-number calling concepts
  • copy-number/copy-ratio-segmentation - segmentation concepts

レビュー

まだレビューはありません。使ってみた感想をお寄せください。

同じリポジトリのスキル

概要と使いどころ

Predicts ADMET properties using ADMETlab 3.0 (119 platform features, including 77 prediction models with modeled-endpoint uncertainty), ADMET-AI, DeepChem MolNet, and chemprop D-MPNN with explicit handling of OECD QSAR principles, applicability domain assessment, calibration, hERG/CYP/AMES endpoints, and PAINS / Lipinski / Ro5 / Veber / BBB druglikeness filters. Use when filtering compounds for drug-likeness, prioritizing leads by predicted safety, or building an in-house ADMET QSAR model.

日本語の概要は準備中です。原文の説明を表示しています。

GPTomics/bioSkills1,2192026年8月15日 更新

Trim PCR primers from aligned reads in amplicon-panel BAMs using samtools ampliconclip. Use when processing SARS-CoV-2 ARTIC, hereditary cancer panels, ctDNA hot-spot panels, or any amplicon assay where primer-derived bases would falsely confirm reference at primer footprints.

日本語の概要は準備中です。原文の説明を表示しています。

GPTomics/bioSkills1,2192026年8月15日 更新

Filter alignments by flags, mapping quality, and regions using samtools view and pysam. Use when extracting specific reads, removing low-quality alignments, or subsetting to target regions.

日本語の概要は準備中です。原文の説明を表示しています。

GPTomics/bioSkills1,2192026年8月15日 更新

Create and use BAI/CSI indices for BAM/CRAM files using samtools and pysam. Use when enabling random access to alignment files or fetching specific genomic regions.

日本語の概要は準備中です。原文の説明を表示しています。

GPTomics/bioSkills1,2192026年8月15日 更新

Read, write, and convert multiple sequence alignment files using Biopython Bio.AlignIO. Supports Clustal, PHYLIP, Stockholm, FASTA, Nexus, and other alignment formats for phylogenetics and conservation analysis. Use when reading, writing, or converting alignment file formats.

日本語の概要は準備中です。原文の説明を表示しています。

GPTomics/bioSkills1,2192026年8月15日 更新

Parse and analyze multiple sequence alignments using Biopython. Extract sequences, identify conserved regions, analyze gaps, work with annotations, and manipulate alignment data for downstream analysis. Use when parsing or manipulating multiple sequence alignments.

日本語の概要は準備中です。原文の説明を表示しています。

GPTomics/bioSkills1,2192026年8月15日 更新

GPTomics のスキルをすべて見る

このスキルの問題を報告する