Version Compatibility
Reference examples tested with: BWA 0.7.17+, VarDict 1.8+, fgbio 2.1+, ichorCNA 0.6.0+, FinaleToolkit 0.9+ (the delfi param was autosomes before 0.9, renamed chrom_sizes), MethylDackel 0.6+, numpy 1.26+, pandas 2.2+, pysam 0.22+, samtools 1.19+
Before using code patterns, verify installed versions match. If versions differ:
- Python:
pip show <package> then help(module.function) to check signatures
- 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.
Liquid Biopsy Analysis Pipeline
"Analyze my liquid biopsy cfDNA data end-to-end" -> Orchestrate UMI-aware preprocessing (fgbio), ctDNA mutation detection (VarDict), tumor fraction estimation (ichorCNA), fragmentomics analysis, and longitudinal monitoring for treatment response.
This is a workflow skill: it owns the chaining decisions and hand-offs, not the internals of any one step. Every step below cross-references the component skill that teaches its mechanism.
The governing principle
A liquid-biopsy result is decided at four seams, two of them set before the sequencer ever runs.
- Pre-analytics is irreversible and is the sensitivity ceiling — set at the blood draw, not the pipeline. WBC lysis dumps high-quality germline/CHIP DNA into the denominator; a 0.5% tumor fraction diluted to 0.1% by a processing delay is a false negative with NO bioinformatic recovery. Tube chemistry (Streck ~7d RT vs EDTA <6h), double-spin, and hemolysis are the FIRST gate the whole pipeline is conditional on. Mixing tube types within a longitudinal series is a batch confound.
- The error-suppression / UMI scheme is committed at library prep and cannot be upgraded after sequencing. Single-strand UMI consensus floors error ~1e-4 to 1e-5 but does NOT remove deamination (C>T) or oxidation (G>T) — those lesions are on the template, so every PCR copy inherits them and the family votes unanimously for the artifact. Only DUPLEX (both-strand concordance) reaches <1e-7. Call on the consensus BAM, never the raw BAM.
- A VAF is undefined without input genome-equivalents, and TF != VAF. 0.1% on 100 GE is noise; on 30,000 GE it is solid. The LoD is set by input GE and assay design (tumor-informed bespoke integrates 16-50 loci to ppm; tumor-naive fixed panel is error/sampling-limited ~0.1-0.5%; ichorCNA sWGS floors ~3% TF). TF ~ 2x VAF only for a clonal, heterozygous, diploid locus — copy number breaks the factor of 2.
- Matched buffy-coat/WBC is the CHIP-subtraction commitment. ~81.6% of cfDNA variants in controls and ~53.2% in cancer patients are CHIP, not tumor (Razavi 2019). Without matched WBC, a tumor-naive cfDNA call is presumptively CHIP-contaminated; a gene-list filter is a weak fallback. Decide at study design.
Made-once commitments
| Commitment | Consequence inherited downstream |
|---|
| Pre-analytics (tube/time-to-plasma/double-spin) | Irreversible TF dilution; the whole pipeline's sensitivity ceiling; consistent across a longitudinal series |
| Error-suppression scheme (single-strand vs DUPLEX) | The achievable error floor; single-strand cannot remove deamination/oxidation; cannot upgrade post-hoc |
| Assay design (tumor-naive vs tumor-informed; panel vs sWGS) + input GE | The LoD and whether it is per-locus or panel-integrated; a VAF is undefined without input GE |
| Matched WBC available? | Whether CHIP is definitively subtracted (WBC) or only gene-list-filtered (weak) |
Pipeline Overview
Pre-analytical QC -> cfDNA Preprocessing -> Fragment QC
↓
┌─────────────────┴─────────────────┐
↓ ↓
sWGS Branch Panel Branch
↓ ↓
ichorCNA VarDict/smCounter2
(Tumor Fraction) (Mutation Detection)
↓ ↓
└─────────────────┬─────────────────┘
↓
Longitudinal Tracking
Step 0: Pre-Analytical QC
def check_preanalytical_quality(sample_metadata):
'''
Pre-analytical factors critical for cfDNA quality.
Requirements:
- Streck tube: up to 7 days at room temperature
- EDTA tube: process within 6 hours
- Avoid hemolysis
- Store extracted DNA at -80C
'''
issues = []
if sample_metadata['tube_type'] == 'EDTA':
if sample_metadata['processing_delay_hours'] > 6:
issues.append('EDTA tube processed > 6 hours - risk of gDNA contamination')
if sample_metadata['hemolysis_score'] > 1:
issues.append('Hemolysis detected - expect cellular DNA contamination')
return issues
Step 1: cfDNA Preprocessing with UMI Consensus
fgbio ExtractUmisFromBam \
--input raw.bam \
--output with_umis.bam \
--read-structure 3M2S+T 3M2S+T \
--single-tag RX
samtools fastq -T RX with_umis.bam | \
bwa mem -C -p -t 8 -Y reference.fa - | \
fgbio ZipperBams --unmapped with_umis.bam --ref reference.fa --output aligned.bam
fgbio GroupReadsByUmi \
--input aligned.bam \
--output grouped.bam \
--strategy adjacency \
--edits 1 \
--family-size-histogram qc/family_sizes.txt
fgbio CallMolecularConsensusReads \
--input grouped.bam \
--output consensus.bam \
--min-reads 1
fgbio FilterConsensusReads \
--input consensus.bam \
--output final.bam \
--ref reference.fa \
--min-reads 2
Step 2: Fragment QC Checkpoint
import pysam
import numpy as np
def verify_cfdna_quality(bam_path):
'''
QC Checkpoint: Verify cfDNA fragment profile.
Expected: peak at ~167bp (mononucleosome)
'''
bam = pysam.AlignmentFile(bam_path, 'rb')
sizes = []
for read in bam.fetch():
if read.is_proper_pair and not read.is_secondary and 0 < read.template_length <= 400:
sizes.append(read.template_length)
bam.close()
sizes = np.array(sizes)
modal_size = np.bincount(sizes).argmax()
mono_frac = np.sum((sizes >= 150) & (sizes <= 180)) / len(sizes)
qc_pass = 150 <= modal_size <= 180 and mono_frac > 0.3
return {
'modal_size': modal_size,
'mononucleosome_fraction': mono_frac,
'qc_pass': qc_pass,
'message': 'Good cfDNA profile' if qc_pass else 'Atypical fragment distribution'
}
Step 3a: Tumor Fraction Estimation (sWGS)
ichorCNA is a command-line script (Rscript scripts/runIchorCNA.R), NOT an importable runIchorCNA() function, and it is preceded by HMMcopy readCounter to bin the BAM. The ~3% tumor-fraction floor is an analytical limit of detection; below it, route to fragmentomics or methylation rather than trusting a low value (see liquid-biopsy/analytical-validation and liquid-biopsy/tumor-fraction-estimation).
readCounter --window 1000000 --quality 20 \
--chromosome "chr1,chr2,chr3,chr4,chr5,chr6,chr7,chr8,chr9,chr10,chr11,chr12,chr13,chr14,chr15,chr16,chr17,chr18,chr19,chr20,chr21,chr22,chrX" \
sample.bam > sample.wig
Rscript scripts/runIchorCNA.R \
--id sample_id --WIG sample.wig \
--gcWig gc_hg38_1000kb.wig --mapWig map_hg38_1000kb.wig \
--centromere GRCh38.GCA_000001405.2_centromere_acen.txt \
--normalPanel HD_ULP_PoN_hg38_1Mb_normAutosomes_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 \
--outDir ichor_results/
Step 3b: Mutation Detection (Targeted Panel)
vardict-java \
-G reference.fa \
-f 0.005 \
-N sample_id \
-b consensus.bam \
-c 1 -S 2 -E 3 -g 4 \
panel.bed | \
teststrandbias.R | \
var2vcf_valid.pl \
-N sample_id \
-E \
-f 0.005 \
> sample.vcf
Step 4: CHIP Filtering
Clonal hematopoiesis (CHIP) is the dominant false-positive source in plasma: ~81.6% of cfDNA variants in controls and ~53.2% in cancer patients trace to white blood cells (Razavi 2019 Nat Med 25:1928). A gene-list filter is a weak fallback; the definitive control is sequencing matched buffy-coat/WBC DNA and subtracting any variant present there. See liquid-biopsy/ctdna-mutation-detection.
CHIP_GENES = ['DNMT3A', 'TET2', 'ASXL1', 'PPM1D', 'JAK2', 'SF3B1', 'SRSF2', 'TP53']
def filter_chip(variants_df, wbc_variants=None, chip_genes=CHIP_GENES):
'''Subtract WBC-matched variants when available; else fall back to a CHIP gene list.'''
if wbc_variants is not None:
key = ['chrom', 'pos', 'ref', 'alt']
wbc_keys = set(map(tuple, wbc_variants[key].itertuples(index=False, name=None)))
in_wbc = variants_df[key].apply(lambda r: tuple(r) in wbc_keys, axis=1)
return variants_df[~in_wbc], variants_df[in_wbc]
chip = variants_df[variants_df['gene'].isin(chip_genes)]
somatic = variants_df[~variants_df['gene'].isin(chip_genes)]
return somatic, chip
Step 5: Fragmentomics Analysis (Optional)
FinaleToolkit (MIT license, not DELFI software) exposes real hyphenated CLI subcommands and an underscored finaletoolkit.frag Python API; delfi GC-corrects the short/long ratio (raw ratios are dominated by GC and sequencing batch). DELFI is a methodology and a company, not a pip install-able tool.
finaletoolkit delfi consensus.bam hg38.chrom.sizes hg38.2bit bins_100kb.bed -g gaps.bed --no-remove-nocov -o sample.delfi.bed
finaletoolkit end-motifs consensus.bam hg38.2bit -o sample.end_motifs.tsv
finaletoolkit mds sample.end_motifs.tsv
from finaletoolkit.frag import delfi
def run_fragmentomics(bam_path, chrom_sizes, reference, bins_bed, gap_bed):
'''GC-corrected DELFI short/long profile (MDS comes from end_motifs().motif_diversity_score()).
Python positional order is (input, chrom_sizes, bins_file, reference_file) - note this differs
from the CLI order (input, chrom_sizes, reference, bins), so pass by keyword to be safe.'''
return delfi(bam_path, chrom_sizes=chrom_sizes, bins_file=bins_bed,
reference_file=reference, gap_file=gap_bed)
Step 6: Longitudinal Tracking
import pandas as pd
import numpy as np
def track_longitudinal(samples_df):
'''
Track ctDNA over treatment.
samples_df columns: [sample_id, timepoint, tumor_fraction, mutations...]
'''
samples_df = samples_df.sort_values('timepoint')
baseline = samples_df.iloc[0]['tumor_fraction']
samples_df['log2_fc'] = np.log2(samples_df['tumor_fraction'] / baseline)
nadir = samples_df['tumor_fraction'].min()
response = 'unknown'
if nadir < 0.001:
response = 'Complete molecular response'
elif nadir < baseline * 0.01:
response = 'Major molecular response (>2 log)'
elif nadir < baseline * 0.5:
response = 'Partial molecular response'
return samples_df, response
Complete Pipeline Script
def run_liquid_biopsy_pipeline(sample_config):
'''
Complete liquid biopsy analysis pipeline.
sample_config: dict with keys:
- bam_file: Input BAM
- data_type: 'swgs' or 'panel'
- reference: Reference FASTA
- bed_file: Panel BED (for panel data)
- output_dir: Output directory
'''
results = {}
if sample_config.get('has_umis'):
preprocessed_bam = preprocess_with_fgbio(sample_config['bam_file'])
else:
preprocessed_bam = sample_config['bam_file']
frag_qc = verify_cfdna_quality(preprocessed_bam)
if not frag_qc['qc_pass']:
print(f"WARNING: {frag_qc['message']}")
results['fragment_qc'] = frag_qc
if sample_config['data_type'] == 'swgs':
results['tumor_fraction'] = run_ichorcna(preprocessed_bam)
elif sample_config['data_type'] == 'panel':
variants = call_variants(preprocessed_bam, sample_config['bed_file'])
wbc = call_variants(sample_config['wbc_bam'], sample_config['bed_file']) if sample_config.get('wbc_bam') else None
somatic, chip = filter_chip(variants, wbc_variants=wbc)
results[] = somatic
results[] = chip
sample_config.get():
results[] = run_fragmentomics(preprocessed_bam)
results
Choosing the branch (tumor-naive vs tumor-informed)
Pipeline-level selection only; mechanism lives in the component skills.
| Situation | Branch | Hand off to |
|---|
| Screening / no known tumor / MCED | tumor-NAIVE (fixed panel, sWGS, or methylation) | liquid-biopsy/methylation-based-detection |
| MRD/recurrence of a KNOWN tumor | tumor-INFORMED bespoke (16-50 clonal variants from tissue WES) | liquid-biopsy/longitudinal-monitoring |
| Tumor fraction from sWGS (0.1-1x) | ichorCNA (copy-number-based, floor ~3% TF) | liquid-biopsy/tumor-fraction-estimation |
| Low-VAF mutations from a deep panel | VarDict / smCounter2 on the UMI-consensus BAM | liquid-biopsy/ctdna-mutation-detection |
| Below mutation LoD, still need signal | fragmentomics (DELFI/end-motifs) or methylation | liquid-biopsy/fragment-analysis |
| Stating/trusting a sensitivity claim | LoB/LoD/LoD95/LoQ, per-locus vs panel-integrated | liquid-biopsy/analytical-validation |
Tumor-informed reaches ppm by integrating across 16-50 loci (beats the single-locus Poisson floor) but is BLIND by design to new/resistance variants; tumor-naive sees any panel variant but is CHIP-dominated and error-limited.
Common Errors
| Symptom | Cause | Fix |
|---|
| False negative on a real low-TF sample | Bad blood draw irreversibly diluted TF | Pre-analytic gate; cannot be rescued downstream |
| Residual C>T / G>T false positives | Single-strand consensus assumed to remove damage | Duplex, or UDG for deamination; flag the substitution spectrum |
| Money spent, no sensitivity gain | Sequenced past molecular saturation | Report unique-molecule coverage; invest in plasma volume / conversion efficiency |
| Blood clones reported as tumor | CHIP not subtracted (esp. TP53/PPM1D, both CHIP and driver) | Matched buffy-coat/WBC sequencing (gene list is a weak fallback) |
| Apparent burden halved or doubled | VAF reported as TF (or vice versa) | State VAF vs TF; TF ~ 2x VAF only for clonal-het-diploid; prefer CNA-based TF when CN is non-neutral |
| A "low" TF trusted below the assay floor | ~3% ichorCNA LoD ignored | Below the floor = below detection; route to fragmentomics/methylation |
| Longitudinal "response" that is a batch shift | Tube/panel/pipeline changed between draws | Same tube, same panel, same pipeline across timepoints |
References
- Razavi P, Li BT, Brown DN, et al (2019) High-intensity sequencing reveals the sources of plasma circulating cell-free DNA variants. Nature Medicine 25:1928-1937. DOI 10.1038/s41591-019-0652-7. (CHIP dominance of cfDNA variants.)
- Adalsteinsson VA, Ha G, Freeman SS, et al (2017) Scalable whole-exome sequencing of cell-free DNA reveals high concordance with metastatic tumors. Nature Communications 8:1324. DOI 10.1038/s41467-017-00965-y. (ichorCNA; ~3% TF LoD.)
- Schmitt MW, Kennedy SR, Salk JJ, et al (2012) Detection of ultra-rare mutations by next-generation sequencing. PNAS 109:14508-14513. DOI 10.1073/pnas.1208715109. (duplex error floor.)
- Reinert T, Henriksen TV, Christensen E, et al (2019) Analysis of plasma cell-free DNA by ultradeep sequencing in patients with stages I to III colorectal cancer. JAMA Oncology 5:1124-1131. DOI 10.1001/jamaoncol.2019.0528. (tumor-informed bespoke MRD.)
Related Skills
- liquid-biopsy/cfdna-preprocessing - UMI/duplex consensus error suppression
- liquid-biopsy/analytical-validation - molecule-counting limits of detection and honest LoD reporting
- liquid-biopsy/ctdna-mutation-detection - low-VAF calling and CHIP subtraction
- liquid-biopsy/tumor-fraction-estimation - ichorCNA tumor fraction from sWGS
- liquid-biopsy/fragment-analysis - fragmentomics features
- liquid-biopsy/methylation-based-detection - methylation detection and tissue-of-origin
- liquid-biopsy/longitudinal-monitoring - serial MRD tracking
