Decides how to preprocess plasma cfDNA sequencing data so the recoverable signal survives - library-prep-aware fragment expectations (dsDNA vs ssDNA/adaptase prep), UMI/duplex consensus with fgbio (ExtractUmisFromBam, GroupReadsByUmi --strategy paired for duplex, CallMolecularConsensusReads vs CallDuplexConsensusReads, FilterConsensusReads min-reads "total s1 s2"), the align->group->consensus->RE-align ordering, and the cfDNA dedup trap where naive coordinate dedup collapses nucleosome-coincident independent molecules. Covers when single-strand consensus suffices vs when duplex is mandatory, the singleton/sensitivity tax at low input, and reading the insert-size histogram as a pre-analytical QC instrument. Use when processing plasma cfDNA reads before fragmentomics, ctDNA mutation calling, or tumor-fraction estimation.
Decides how to preprocess plasma cfDNA sequencing data so the recoverable signal survives - library-prep-aware fragment expectations (dsDNA vs ssDNA/adaptase prep), UMI/duplex consensus with fgbio (ExtractUmisFromBam, GroupReadsByUmi --strategy paired for duplex, CallMolecularConsensusReads vs CallDuplexConsensusReads, FilterConsensusReads min-reads "total s1 s2"), the align->group->consensus->RE-align ordering, and the cfDNA dedup trap where naive coordinate dedup collapses nucleosome-coincident independent molecules. Covers when single-strand consensus suffices vs when duplex is mandatory, the singleton/sensitivity tax at low input, and reading the insert-size histogram as a pre-analytical QC instrument. Use when processing plasma cfDNA reads before fragmentomics, ctDNA mutation calling, or tumor-fraction estimation.
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: fgbio flag semantics drift across major versions - in the is a permissive PRE-filter (fgbio issue #1009), and in the () takes up to three values "total strand1 strand2" where if values two and three differ the more stringent value must come first. Confirm both against the installed before scripting.
CallDuplexConsensusReads
--min-reads
FilterConsensusReads
--min-reads
-M
fgbio --help
cfDNA Preprocessing
"Preprocess my plasma cfDNA reads" -> Convert reads to consensus molecules and an analysis-ready BAM without destroying the fragment-size signal or collapsing independent molecules.
Python: shell the fgbio/bwa chain with subprocess; read insert sizes with pysam for QC
The Single Most Important Modern Insight -- Pre-Analytics and Library Prep Set the Ceiling; Consensus Is Error Suppression, Not Just Dedup
What a cfDNA assay can ever see is fixed before any bioinformatics runs. cfDNA is not randomly sheared - it is a nucleosome footprint with a mononucleosome mode at ~167 bp and a ctDNA-enriched short tail (134-144 bp), so fragment length is structured biological signal, not noise to be normalized away. Two upstream choices determine what survives: (1) the blood draw and plasma prep (a delayed draw dumps leukocyte gDNA into the denominator and irreversibly dilutes tumor fraction - no algorithm recovers it; defer that QC to liquid-biopsy/analytical-validation), and (2) the library chemistry (dsDNA ligation polishes native ends and discards the sub-100 bp population; ssDNA/adaptase prep recovers short and damaged molecules). Reverse engineering a fragment the prep threw away, or a molecule a bad draw never delivered, is impossible.
The second reframe: UMI/duplex consensus is error suppression, not merely duplicate removal. Grouping reads by UMI and majority-voting a consensus erases PCR/sequencing errors that are not shared across a family - but a single-strand consensus votes unanimously for damage (C->T deamination, G->T 8-oxoG) that was on the template before amplification. Only DUPLEX consensus, requiring the same base on both independently-copied strands, removes those lesions. Single-strand consensus cannot. Choosing simplex vs duplex is choosing an error floor, and it trades against molecular recovery at low input (the singleton tax below).
requiring duplicate observation discards the only evidence for a low-VAF variant
No UMIs at all, quantitative readout
do NOT coordinate-dedup; document the bias (see Failure Modes)
nucleosome-positioned ends make naive dedup delete real molecules
Methodology evolves: confirm current fgbio flag semantics and prep-vendor size behavior against live docs before committing a pipeline.
The fgbio Consensus Pipeline
Goal: Turn UMI-tagged raw reads into error-suppressed consensus molecules with correct coordinates.
Approach: Extract UMIs into the RX tag, align (consensus needs coordinates to group), group by UMI + approximate position, call consensus (which emits UNMAPPED reads because the consensus sequence differs from any input read), RE-align the consensus, then apply the real quality gate with FilterConsensusReads. The two-pass alignment is mandatory.
# 1. Extract inline UMIs into RX. Read structure tokens: M=UMI, S=skip/stem, T=template, +=all remaining.# 6M11S+T per end = 6 bp UMI, 11 bp stem (the S that bleeds into T if omitted), rest = insert.
fgbio ExtractUmisFromBam --input raw.unmapped.bam --output with_umis.bam \
--read-structure 6M11S+T 6M11S+T --single-tag RX
# 2. Align. -Y soft-clips supplementaries so tag-bearing short-fragment sequence is not dropped.
bwa mem -t 8 -Y reference.fa with_umis.bam | samtools sort -o aligned.bam -
samtools index aligned.bam
# 3. Group by UMI. adjacency = simplex default; paired = MANDATORY for duplex (reconstructs strand pairing).
fgbio GroupReadsByUmi --input aligned.bam --output grouped.bam \
--strategy paired --edits 1 # use --strategy adjacency for single-strand UMIs# 4a. SIMPLEX: single-strand consensus, --min-reads takes ONE value.
fgbio CallMolecularConsensusReads --input grouped.bam --output consensus.unmapped.bam --min-reads 1
# 4b. DUPLEX: --min-reads here is a permissive PRE-filter (fgbio #1009) - call low, filter later.
fgbio CallDuplexConsensusReads --input grouped.bam --output consensus.unmapped.bam --min-reads 1
# 5. RE-align: consensus reads are emitted UNMAPPED by design. Re-map, then ZipperBams# transfers the consensus/UMI tags from the unmapped BAM onto the new alignments.
samtools fastq consensus.unmapped.bam | bwa mem -t 8 -Y -p reference.fa - \
| fgbio ZipperBams --unmapped consensus.unmapped.bam --ref reference.fa \
| samtools sort -o consensus.bam -
# 6. The REAL quality gate. --min-reads "total strand1 strand2"; "2 1 1" = true duplex (both strands seen).# If values two and three differ, the more stringent must come first (e.g. "6 3 0", not "0 3").
fgbio FilterConsensusReads --input consensus.bam --output filtered.bam --ref reference.fa \
--min-reads 2 1 1 --max-read-error-rate 0.025 --max-base-error-rate 0.1 \
--min-base-quality 40 --reverse-per-base-tags
Key flags: --strategy paired is non-negotiable for duplex (adjacency cannot reconstruct A/B strand pairing). --min-reads 1 1 0 on the filter accepts single-strand consensus too (one strand may be absent); 2 1 1 requires both strands = true duplex. --reverse-per-base-tags makes per-base depth/error tags read in genomic orientation after alignment.
Reading the Insert-Size Histogram (QC)
Goal: Use the fragment-length distribution as a pre-analytical instrument before trusting any downstream number.
Approach: Tabulate proper-pair template lengths from the BAM, locate the mode, and compare the shape against the expected nucleosome footprint - the mode, the sawtooth, and the long-fragment fraction each diagnose a specific failure.
import pysam
import numpy as np
definsert_size_qc(bam_path, max_size=600):
'''Summarize cfDNA fragment lengths as a QC readout.'''
bam = pysam.AlignmentFile(bam_path, 'rb')
sizes = [r.template_length for r in bam.fetch()
if r.is_proper_pair andnot r.is_secondary and0 < r.template_length <= max_size]
bam.close()
sizes = np.array(sizes)
long_frac = np.mean(sizes > 250) # excess >250 bp signals gDNA/leukocyte-lysis contaminationreturn {'n': len(sizes), 'mode_bp': int(np.bincount(sizes).argmax()),
'median_bp': float(np.median(sizes)), 'frac_over_250bp': float(long_frac)}
A mode at ~167 bp with a ~10.4 bp sawtooth below it is healthy. A mode drifting up plus excess mass >250 bp is gDNA contamination. A mode ~10 bp low is the adaptase chemistry signature, not contamination. A ~120-130 bp spike is adapter dimer (trim adapters before reading the histogram).
Per-Method Failure Modes
The dedup trap (naive coordinate dedup on no-UMI cfDNA)
Trigger: running Picard MarkDuplicates / samtools markdup on cfDNA without UMIs. Mechanism: cfDNA ends pile up non-randomly at nucleosome/linker boundaries, so many independent molecules genuinely share the same ~167 bp start+end; coordinate dedup declares them PCR duplicates and keeps one. Symptom: deflated unique-molecule counts and depressed VAF for true low-frequency variants - worst exactly where coverage and nucleosome positioning are strongest. Fix: use UMIs and group by family; if none, do not dedup by position alone for quantitative readouts and document the bias. This is the strongest single argument for putting UMIs on a cfDNA assay.
Strict filtering on the consensus caller instead of FilterConsensusReads
Trigger: setting CallDuplexConsensusReads --min-reads high. Mechanism: it is a PRE-filter (fgbio #1009) that discards data before the real filter sees it. Symptom: lower molecular recovery than expected with no specificity gain. Fix: call permissively (--min-reads 1), apply strictness in FilterConsensusReads.
adjacency grouping on duplex data
Trigger:--strategy adjacency for a duplex library. Mechanism: adjacency cannot link the two strands of one duplex (which carry the UMI in opposite order). Symptom: duplex consensus finds no two-strand molecules; output collapses to simplex. Fix:--strategy paired.
Blanket MAPQ filtering kills short ctDNA fragments
Trigger: a global MAPQ >= 30/60 filter. Mechanism: a 40-80 bp insert has fewer anchoring bases and lands in the low-MAPQ tail even when correctly placed. Symptom: the short, ctDNA-enriched fragments are preferentially deleted - the molecules size-selection was meant to keep. Fix: filter MAPQ with the size distribution in mind; use a gentler threshold for fragmentomics.
Reflexive --min-reads >= 2 at low input (the singleton tax)
Trigger: requiring duplicate observation per family on a picogram-input library. Mechanism: a large fraction of unique cfDNA molecules are sequenced once (singletons); requiring two reads discards genuine, often the only, evidence for a variant. Symptom: sensitivity loss disguised as quality. Fix: for detection favor recovery (call permissively); reserve 2 1 1 for genotyping/MRD where error suppression dominates.
Fragmentomics on a size-selected library
Trigger: computing end-motif/ratio/VAF features after in-vitro (Pippin) or in-silico size selection. Mechanism: selection conditions on length, so length-derived features are biased by construction. Symptom: distorted end-motif spectra and fragment-ratio features that do not reproduce. Fix: size-select for detection sensitivity only; never report length-derived features from a length-selected library; keep the selection step in metadata.
Quantitative Thresholds
Threshold
Source
Rationale
Mononucleosome mode 167 bp; ~10.4 bp sub-mode periodicity
Snyder 2016 Cell 164:57
core particle (~147 bp) + linker; periodicity = helical pitch of nucleosome-bound DNA, a QC sanity check
TF/CTCF-footprint short fragments 35-80 bp
Snyder 2016 Cell 164:57
sub-nucleosomal protection; only surfaces with ssDNA prep
ctDNA principal length 134-144 bp vs 167 bp germline
Underhill 2016 PLoS Genet 12:e1006162
tumor fragments shorter; BRAF V600E mutant 132-145 bp vs WT 165 bp
ctDNA-enrichment selection windows 90-150 bp (and 250-320 bp)
Mouliere 2018 Sci Transl Med 10:eaat4921
selecting the short window enriches mutant allele fraction at no extra sequencing cost
per-read vs per-base vs fraction/count switch (<1.0 = fraction, >=1.0 = count)
gDNA-contamination signature: excess mass >180-250 bp / ladder distortion
Snyder 2016 (biology); pre-analytical QC
leukocyte-lysis gDNA is long; a "too clean, too long" library is contaminated, not pristine
Common Errors
Error / symptom
Cause
Solution
Consensus BAM has garbage coordinates
skipped re-alignment (consensus emitted unmapped)
RE-align consensus reads before FilterConsensusReads
UMI bases corrupt mapping / consensus
omitted the S skip in the read structure
include the stem: e.g. 6M11S+T not 6M+T
Duplex run finds no two-strand molecules
grouped with --strategy adjacency
regroup with --strategy paired
FilterConsensusReads rejects valid strand spec
--min-reads 0 3 (less stringent first)
put the more stringent value first: 3 0 (or 6 3 0)
Deflated VAF / low unique-molecule count
coordinate dedup on no-UMI cfDNA
use UMI families; never naive-dedup quantitative cfDNA
Short ctDNA fragments missing after filtering
blanket high MAPQ filter
lower MAPQ threshold; account for short-fragment mapping
Apparent ~10 bp size shift "needs correcting"
adaptase (Swift/Accel-1S) chemistry signature
record the prep; do not correct a real chemistry effect
References
Snyder MW, Kircher M, Hill AJ, Daza RM, Shendure J. 2016. Cell-free DNA comprises an in vivo nucleosome footprint that informs its tissues-of-origin. Cell 164:57-68. -- 167 bp mode, ~10.4 bp periodicity, 35-80 bp TF/CTCF footprints; the nucleosome-positioning basis of the dedup trap.
Underhill HR, Kitzman JO, Hellwig S, Welker NC, Daza R, Baker DN, Gligorich KM, Rostomily RC, Bronner MP, Shendure J. 2016. Fragment length of circulating tumor DNA. PLoS Genet 12:e1006162. -- ctDNA 134-144 bp vs 167 bp germline.
Mouliere F, Chandrananda D, Piskorz AM, Moore EK, Morris J, Ahlborn LB, Mair R, Goranova T, Marass F, Heider K, et al. 2018. Enhanced detection of circulating tumor DNA by fragment size analysis. Sci Transl Med 10:eaat4921. -- 90-150 bp and 250-320 bp ctDNA-enrichment windows.
Burnham P, Kim MS, Agbor-Enoh S, Luikart H, Valantine HA, Khush KK, De Vlaminck I. 2016. Single-stranded DNA library preparation uncovers the origin and diversity of ultrashort cell-free DNA in plasma. Sci Rep 6:27859. -- ssDNA prep recovers the ultrashort mitochondrial/microbial fraction dsDNA prep discards.
Newman AM, Lovejoy AF, Klass DM, Kurtz DM, Chabon JJ, Scherer F, Stehr H, Liu CL, Bratman SV, Say C, et al. 2016. Integrated digital error suppression for improved detection of circulating tumor DNA. Nat Biotechnol 34:547-555. -- barcode + in-silico error suppression are complementary (~3x each).
Kennedy SR, Schmitt MW, Fox EJ, Kohrn BF, Salk JJ, Ahn EH, Prindle MJ, Kuong KJ, Shen JC, Risques RA, Loeb LA. 2014. Detecting ultralow-frequency mutations by Duplex Sequencing. Nat Protoc 9:2586-2606. -- both-strand concordance and the <1e-7 error floor that single-strand consensus cannot reach.
Related Skills
analytical-validation - the LoD/molecule-counting framework input quality feeds
fragment-analysis - fragmentomics consumes the preprocessed fragment ends
