DupCaller

DupCaller is a tool for calling somatic mutations and calculating somatic mutational burden from barcoded error-corrected next generation sequencing (ecNGS/duplex sequencing) data with matched normal (e.x. NanoSeq, UDSeq).

Prerequisites

DupCaller requires python>=3.11 to run (setup.py pins numpy==2.3.4/scipy==1.16.2, both of which require >=3.11). Earlier versions are not supported. The complete DupCaller pipeline also requires the following tools for data preprocessing. The Nextflow pipeline (nextflow/) runs exactly these versions and commands. Other versions may work but have not been validated.

bwa-mem2 2.3 (https://github.com/bwa-mem2/bwa-mem2), run as bwa-mem2 mem -C -T 0 (Step 3). (The v2.3 release still prints 2.2.1 from bwa-mem2 version and in the BAM @PG header.)
GATK 4.3.0.0 (https://github.com/broadinstitute/gatk/releases), for MarkDuplicates (Step 4)
samtools, for sorting and indexing BAMs
htslib, which provides tabix and bgzip for compressed, indexed genomic files (recommended installation: conda install bioconda::htslib)
PERF (Pattern-based Exhaustive Repeat Finder) v0.4.6 (MIT license) is bundled with DupCaller, runs on the same Python, and is installed as the PERF command, so it needs no separate installation.

Installation

pip Installation

The tool uses pip for installing scripts and prerequisites. We recommend creating a new environment to install DupCaller:

conda create -n DupCaller python=3.12 bioconda::htslib

To install DupCaller, simply clone this repository and install via pip:

conda activate DupCaller
git clone https://github.com/AlexandrovLab/DupCaller.git
cd DupCaller
pip install .

Docker / Singularity

A pre-built Docker image is available on Docker Hub at yuhecheng62/dupcaller:1.2.11.

Pull and run with Singularity:

Pull the image from Docker Hub (only needed once):

singularity pull dupcaller-1.2.11.sif docker://yuhecheng62/dupcaller:1.2.11

For quick verification:

singularity exec dupcaller-1.2.11.sif DupCaller.py --help

For installation-free execution of DupCaller commands, run all DupCaller.py commands with singularity exec and binding of current directories:

singularity exec --bind $(pwd):$(pwd) dupcaller-1.2.11.sif DupCaller.py {your commands}

Pipeline

Step 1: Index Reference Genome

DupCaller uses a numpyrized reference genome to perform memory-efficient reference fetching, trinucleotide context lookup, and repeat (homopolymer/short tandem repeat, STR) annotation used to improve indel calling near repetitive regions.

STR annotation comes from a single tsv produced by PERF. First, run PERF against the reference FASTA:

PERF -m 2 -M 10 -u 2 -i reference.fa -o repeats.tsv

-m/-M — min/max repeat unit (motif) size to search for. -m 2 -M 10 covers 2–10bp-unit STRs; going much higher increases PERF's runtime and memory non-trivially. Homopolymers (unit length 1) are not taken from this tsv: index derives them directly from the reference sequence and ignores any unit-length-1 rows, so -m 1 only adds runtime.
-u 2 (--min-units) — minimum number of repeat copies to report.

Then index the reference, passing PERF's tsv directly:

DupCaller.py index -f reference.fa -rt repeats.tsv

The command will generate five h5 files in the same folder as the reference: {reference}.ref.h5, {reference}.tn.h5, {reference}.hp.h5, {reference}.str.h5, and {reference}.dbs.h5 — numpyrized reference sequences, trinucleotide contexts, homopolymer annotations, STR (unit length >=2) annotations, and per-position dinucleotide (DBS) classes, respectively.

For human reference genome hg38 and mouse reference genome mm39, we provided pre-built indexes and resource files in the Resources.

The reference also needs a samtools FASTA index and a bwa-mem2 index (for Step 3), built once:

samtools faidx reference.fa
bwa-mem2 index reference.fa   # writes reference.fa.{0123,amb,ann,bwt.2bit.64,pac}; ~90 GB RAM for GRCh38 (28 bytes per base)

Parameters

Short	Long	Description
-f	--reference	Reference genome fasta file (required)
-rt	--repeatTsv	PERF-format repeat tsv (chrom, start, end, motif, length, strand, num_units, motif_repeat) for the reference (required). Repeat unit length and repeat count are read directly from the motif and num_units columns; rows with unit length 1 are ignored.

Step 2: Trim Barcodes

DupCaller.py trim extracts 5-prime barcodes from paired-end FASTQs:

DupCaller.py trim -i read1.fq -i2 read2.fq -p barcode_pattern -o sample_name

read1.fq / read2.fq — FASTQ files from read 1 and read 2. Both unzipped and gzip-compressed files are supported.
barcode_pattern — pattern of barcodes starting from the 5-prime end, with N representing a barcode base and X representing a skipped base (similar notation to UMI-tools). For example, NanoSeq uses a 3-base barcode followed by 4 constant bases, so the pattern should be NNNXXXX.
Barcodes from a fixed list, of varying length (e.g. Tn5-based protocols such as META-CS, where each read starts with one of a set of 11–13 nt transposon barcodes followed by the 19 nt ME): put B in the pattern for "one barcode from the list" and pass the list with -bl/--barcode-list (one barcode per line). Positions after B count from the end of the matched barcode, e.g. -p BXXXXXXXXXXXXXXXXXXX takes the barcode and skips the 19 nt ME. Each read's barcode may have up to -mm/--max-mismatch mismatches (default 1); the barcode with the fewest mismatches wins, then the shortest. Up to 1 mismatch is matched with lookup tables (fast, any list size); larger limits walk all listed barcodes at once, dropping each as soon as it exceeds the limit. The listed (corrected) sequence is written to DB. Read pairs whose read 1 or read 2 matches no barcode, or matches two equally well, are dropped and counted in the summary trim prints.
sample_name — prefix of output paired FASTQs. After the run completes, {sample_name}_1.fastq and {sample_name}_2.fastq will be generated. Barcodes are appended to each read name ({original_read_name}_{read1_barcode}+{read2_barcode}, for optical-duplicate read-name matching only) and also written as a DB:Z:{read1_barcode}-{read2_barcode} FASTQ comment; the comment is promoted to a BAM tag by bwa-mem2 mem -C (Step 3) and is what DupCaller.py call actually reads to assign duplex family membership — it never parses the read name for barcodes.

If the matched normal is prepared in the same way as the sample, apply trimming with the same scheme to the matched normal FASTQs. For traditional bulk normal, trimming is not needed.

Step 3: Align Reads

Align the trimmed FASTQs of both sample and matched normal with bwa-mem2 (bwa-mem2 index from Step 1 must already exist). GATK requires read group fields ID, SM, and PL, so adding those tags during alignment is recommended. FASTQ tags must be kept — for bwa-mem2 mem (as for bwa mem) this requires the -C option. The Nextflow pipeline and the mock-pipeline regression test also use -T 0, which outputs every alignment regardless of its score, so low-scoring alignments are filtered by DupCaller (MAPQ/AS-XS) rather than dropped by the aligner.

bwa-mem2 mem -C -T 0 -t {threads} -R "@RG\tID:{sample_name}\tSM:{sample_name}\tPL:ILLUMINA" \
    reference.fa {sample_name}_1.fastq {sample_name}_2.fastq \
    | samtools sort -@ {threads} > {sample_name}.bam
samtools index -@ {threads} {sample_name}.bam

threads — number of cores used for aligning
reference.fa — reference genome FASTA file, with its bwa-mem2 index next to it
{sample_name}_1.fastq / {sample_name}_2.fastq — trimmed FASTQ files from Step 2

Step 4: Mark Duplicates

Run GATK MarkDuplicates on sample and matched-normal BAMs. Optical and PCR duplicates must be treated differently in ecNGS variant calling:

Set --TAGGING_POLICY OpticalOnly to differentiate optical from PCR duplicates.
Set --DUPLEX_UMI true for duplex UMI handling.
Set --READ_NAME_REGEX as shown below, because the read names were modified in Step 2.

Note: The published analyses used GATK 4.3.0.0. Older GATK releases do not have the --DUPLEX_UMI flag.

gatk MarkDuplicates \
    -I sample.bam -O sample.mkdped.bam -M sample.mkdp_metrics.txt \
    --READ_NAME_REGEX "(?:.*:)?([0-9]+)[^:]*:([0-9]+)[^:]*:([0-9]+)[^:]*$" \
    --DUPLEX_UMI --TAGGING_POLICY OpticalOnly --BARCODE_TAG DB

Step 5: Call Variants

After preprocessing, run DupCaller.py call to call somatic mutations. Usage depends on your experimental design:

Whole genome / whole exome / reduced genome (e.g. NanoSeq) with a matched normal:

DupCaller.py call -b ${sample}.mkdped.bam -f reference.fa -o {output_prefix} \
    -p {threads} -n {normal.bam} -g germline.vcf.gz -m noise_mask.bed.gz

Mutagenesis panel without a matched normal:

DupCaller.py call -b ${sample}.mkdped.bam -f reference.fa -o {output_prefix} \
    -p {threads} -g germline.vcf.gz -m noise_mask.bed.gz -maf 0.1

Note: DupCaller partitions jobs by genomic region; multithreading is less effective for small targeted panels. In this case, use at most one thread per distinct targeted region.

Non-default barcode tags: BAMs produced by DupCaller's own trim (Step 2) carry a shared DB:Z:{bc1}-{bc2} tag per read pair, and this is what call reads by default. If instead you are calling directly from a BAM that stores the duplex barcode differently, pass -bc/--barcode TAG,NORMALIZE,SEP: TAG is the bam tag holding the barcode, and NORMALIZE says whether it still needs splitting/reordering. For example, the Sanger NanoSeq pipeline tags each read with a single RB tag that already identifies the duplex family (e.g. chrom,fragment_start,fragment_end,bc1,bc2) independent of read orientation — pass -bc RB,0,- to read it directly (SEP is ignored in this case). DupCaller's own DB:Z:{bc1}-{bc2} tag is the NORMALIZE=1 case (the default, -bc DB,1,-): it's an unnormalized bc1<SEP>bc2 pair that still needs read-orientation-based reordering to group both strand-copies of a molecule together.

Input validation: DupCaller checks for the existence of all required files (BAM, reference, h5 index files, and optional files) before starting analysis, providing clear error messages for missing files.

Multi-threading: Coverage files from different threads are automatically merged post-processing with region-aware boundary detection and tabix indexing.

See the Results section for descriptions of all output files.

Parameters

Required

Short	Long	Description
-b	--bam	BAM file of ecNGS data
-f	--reference	Reference genome FASTA file
-o	--output	Output directory; created if needed. Files inside are named`{sample}_*`, where `{sample}` is the directory's basename

Short	Long	Description	Default
-r	--regions	Contigs to consider for variant calling. For non-human species, set accordingly (e.g. for mouse:`-r chr{1..19} chrX`)	chr{1..22} chrX
-g	--germline	Indexed germline VCF with AF field	None
-p	--threads	Number of threads	1
	--minChunkLength	Minimum chunk length in bases; short contigs stay unsplit. May use fewer workers than requested. Also accepts`--min-chunk-length`.	10000
-n	--normalBams	BAM file(s) of matched normals. When unavailable, set`-maf` to an appropriate value (e.g. 0.1)	None
-m	--noise	BED interval file(s) masking noisy positions	None
-R	--regionfile	Inclusive BED file specifying target regions	None
-maf	--maxAF	Maximum allele fraction to call a somatic mutation. Must be set when matched normal (`-n`) is unavailable	1
-tt	--trimF	Ignore mutations less than n bp from template ends	7
-tr	--trimR	Ignore mutations less than n bp from read ends	7
-bc	--barcode	Molecular/duplex barcode read tag, as`TAG,NORMALIZE,SEP` (see note above)	DB,1,-

Optional

The effect of changing these parameters should be evaluated before implementation.

Short	Long	Description	Default
--naf		Maximum VAF in matched normal for a mutation to be called	0.01
--rescue	-res	Output discarded variants with reason in the FILTER field	False
-nm	--nmflt	Drop a read family when at least half of its reads have an edit distance (NM) above this value	5
-ax	--minMeanASXS	Minimum mean AS-XS alignment score difference for a read group to be considered	50
-gaf	--germlineAfCutoff	Skip positions with germline AF above this threshold	0.001
-d	--minNdepth	Minimum coverage in normal for called variants	10
-sc	--skipCoveragePass	Skip round 2 (the coverage-only pass): no duplex depth / per-locus`coverage.bed.gz`, no duplex-family composition or by-duplex-group output files, and masked candidates that only clear the final lFDR-refined threshold get no real depth. Only round 1 calling and lFDR-threshold determination run; VCFs and a trimmed `_stats.txt` are still written. Roughly halves total runtime	False

Advanced

These are variant calling model parameters; adjustment is unnecessary for general use.

Short	Long	Description	Default
-E	--errprefix	Prefix for all six error files ({prefix}.amp.tn.srd.txt, {prefix}.amp.hp.txt, {prefix}.amp.str.txt, {prefix}.dmg.tn.txt, {prefix}.dmg.hp.txt, {prefix}.dmg.str.txt); overrides the default (output prefix)	None
-lo	--learnOnly	Stop after estimating/writing the error-rate files (ERROR/ dir); skip variant calling entirely	False
-lfdr	--lfdrThreshold	Target per-channel lFDR (max of`(1-mutation_rate)/(LR*mutation_rate+1-mutation_rate)` over that channel's PASS calls, i.e. its weakest surviving call); channels above this get their LR threshold raised via a single closed-form solve using the round-1 mutation-rate estimate (non-iterative -- no re-simulation/re-estimation loop)	0.05
-a	--pseudocount	Pseudocount (> 0) for each channel's mutation-rate solve: a pseudo-mutations and a pseudo-reference sites, `mu*(Eeff+2a) = sum of site posteriors + a`, which always has exactly one root between 0 and 1. A channel with no effective coverage and no sites gets mutation rate 0. Also the pseudocount of the SBS amplification-error (SRD) fit	0.5
-mlr	--minLR	Minimum log10 LR for a PASS call: floor on every channel's FDR-refined LR threshold, used for calling and for the detection-power simulation behind the burden's sensitivity correction	5
-mq	--mapq	Minimum MAPQ for an alignment to be considered	40
-w	--windowSize	Genomic window size for coverage calculation and BAM partitioning	100000
-bq	--minBq	Bases with quality below this value are zeroed out and excluded from variant calling/learning (usable iff BQ >= minBq)	18
-aq	--minAltQual	Minimum summed per-strand consensus base quality at a position for SBS damage-rate learning to use it	90
--minRef		Minimum per-position read depth (ref+alt combined) for damage-rate learning to use it; the stricter of`--minRef`/`--minAlt` is applied, since the underlying check doesn't separate ref vs alt counts	3
--minAlt		Minimum per-position read depth (ref+alt combined) for damage-rate learning to use it; the stricter of`--minRef`/`--minAlt` is applied, since the underlying check doesn't separate ref vs alt counts	3
--srdMinRead		Minimum reads on a single strand (F1R2 or F2R1) for that strand to be included in SBS single-read-damage (SRD) rate learning; evaluated independently per strand, and also gates the minimum BQ-qualifying bases required for a site to be counted	3
--ssmMinRead		Minimum reads required on EACH strand (F1R2 and F2R1 both) for a duplex family to be considered for single-strand-mutation (SSM)/damage-rate calculation	3
-pt	--p_threshold	Strand independence filter threshold for SBS and indel calls: a PASS call fails with FILTER `strand_independence` if either strand's p-value (INFO `MSP`) is <= this. Per strand, not corrected for the number of calls. Must be in [0, 1); 0 disables the filter (MSP is still reported). See Strand Independence Filter	0.05
-id	--indelbed	Indel enhanced Panel of Normals (ePoN) for indel calling	None
-rt	--regionst	Contigs to consider for error-profile training, if different from`-r`/`--regions`	same as`--regions`
-pd	--maxPileupDepth	Maximum depth for samtools mpileup	1000000
-mr	--muterateprefix	Reuse per-channel mutation-rate tables from an earlier run instead of estimating them from this run's own candidates (useful when calling a small region). Pass`{out}/tmp/{sample}` from that run, which holds `{sample}_sbs96_rate_n1.txt` and `{sample}_indel_rate_by_hp_str.txt` (tables written before the per-run indel Eeff change are refused). LR thresholds are still solved against this run's `-lfdr`	None
	--seed	RNG seed for the Monte Carlo detection-power simulation, for reproducible results across runs and`-p` values. If unset, a random seed is generated and recorded in `{sample}_call_params.log`	random

Germline and Noise Masks

The -m option accepts BED files and will ignore any mutations overlapping an excluded locus. Any custom BED file can be used as input.

Low-Coverage Error-Profile Fallback

DupCaller ships a set of reference error profiles (src/ERROR/fallback_latest.*.txt, installed with the package). When a context in the sample's own learned error profile has fewer than 1000 observations (a trinucleotide row, a homopolymer length × base, or an STR length bin), any error type in that context with zero observed counts takes its rate from the bundled profile instead of the flat pseudocount prior. Error types that were observed, and well-sampled contexts, use the sample's own rates, with one exception: a context whose error rates sum past 1 (after this fallback and the homopolymer monotonicity step) takes the bundled profile's whole row instead. This applies to the SBS amplification (SRD) and damage matrices and to the homopolymer/STR indel matrices.

Strand Independence Filter

The SBS and indel likelihoods compare alt reads against non-alt reads on each strand, so a read family where some reads carry the alt and others do not (for example, 3 alt and 2 reference reads) can score almost as well as a clean one. Such mixed families usually come from an early-cycle PCR error or from two molecules sharing a barcode.

For every PASS SBS and indel call, DupCaller tests each strand (F1R2, F2R1) separately. Under the null hypothesis the strand is a clean mutant and every non-alt read is an error:

SBS: each covering read shows a non-alt base with probability 10^(-BQ/10) plus the amplification error rate for that trinucleotide context. Every read covering the position is used at its actual base quality (no --minBq cut), and the p-value is Poisson-binomial over the reads.
Indel: the reads are the informative ones (showing the indel or spanning the reference allele, with base quality at or above --minBq). Each shows the reference with the learned amplification error rate that undoes the indel in the mutant molecule's own context: for example, a 1 bp deletion from a run of 7 reverts by a +1 error in a run of 6, and a base inserted next to an unrelated run reverts by losing an isolated base. The learned indel rates already include sequencing error, so there is no separate base-quality term; the p-value is binomial.

The p-value is the probability of seeing at least the observed number of non-alt reads. Both strand p-values are written to INFO MSP (F1R2,F2R1) on every SBS and indel call that was PASS going into the filter (so PASS and strand_independence records); every other record gets .. If either strand's p-value is <= --p_threshold (default 0.05), the call gets FILTER strand_independence, is written to the fail VCF, and is not counted by estimate. A DBS whose own read family has a failed SBS gets the same FILTER, and so does that DBS's other SBS, so neither base is counted as a standalone SBS.

A strand whose reads all show the alt has p = 1 and never fails the test, so the filter only acts on families that already contain discordant reads (and on the other SBS of a DBS that failed). For SBS with binned base qualities, a single low-quality (BQ 11) discordant read gives p of about 0.08 and passes, while two such reads, or one high-quality discordant read, fail at 0.05. For indels outside long homopolymers, learned reversion rates are small enough that a single reference read usually fails; in long homopolymers, where slippage is common, one reference read can pass. The threshold is not corrected for the number of calls.

Step 6: Estimate Mutational Burden

After mutation calling, run burden estimation on the call output directory (-i is the same path that was passed to call -o):

DupCaller.py estimate -i sample -f reference.fa -r chr{1..22} chrX

Adjust the -r regions according to the reference genome used.

Per-Gene Coverage

For dNdScv coverage correction, DupCaller can output mean duplex depth per gene using the -gb option:

DupCaller.py estimate -i sample -f reference.fa -r chr{1..22} chrX -gb {target}.bed

The gene BED should have the fourth column formatted as {gene_name}_{exon_number} (e.g. tp53_1). It must be bgzip-compressed and tabix-indexed (bgzip {target}.bed && tabix -p bed {target}.bed.gz).

Re-estimation for Specific Regions

To re-estimate trinucleotide-corrected mutational burden in specific regions without re-running variant calling, use -rb:

DupCaller.py estimate -i sample -f reference.fa -r chr{1..22} chrX -rb {re_estimate}.bed.gz

As with -gb, the BED file must be bgzip-compressed and tabix-indexed.

Parameters

Required

Short	Long	Description
-i	--prefix	Output directory of the`call` command (the `-o` value)
-f	--reference	FASTA file of reference genome (its h5 index files must sit next to it)

Optional

Short	Long	Description	Default
-r	--regions	Contigs to consider for trinucleotide calculation	chr{1..22} chrX
-ot	--outTrinuc	Write the reference trinucleotide composition to this file, for reuse with`-ft`	None
-ft	--refTrinuc	Load a trinucleotide composition written by`-ot` instead of rescanning the reference (`-f` is still required)	None
-d	--dilute	Set when sample and matched normal come from the same starting DNA material: SNVs with a tumor alt allele count (`AC`) above 1 are dropped when their tumor and normal allele counts differ significantly (Barnard's exact test, p <= 0.05), and kept calls are written to `SBS/{sample}_sbs_flt.vcf`	False
-gb	--genebed	bgzipped, tabix-indexed gene BED file for per-gene coverage calculation	None
-rb	--reestimatebed	bgzipped, tabix-indexed BED file for burden re-estimation in specific regions	None

Step 7: Summarize Across Samples

After running estimate on all samples, use DupCaller.py summarize to collate burden metrics and SBS96 profiles into multi-sample tables:

DupCaller.py summarize -i sample1 sample2 sample3 -o cohort_summary.txt

This reads {sample}/{sample}_stats.txt, {sample}/SBS/{sample}_sbs_burden.txt, {sample}/INDEL/{sample}_indel_burden.txt, {sample}/DBS/{sample}_dbs_burden.txt, and {sample}/SBS/{sample}_sbs_96_corrected.txt from each sample folder (all five are required) and writes four output files:

File	Description
`cohort_summary.txt`	One row per sample with SBS, indel, and DBS burden metrics and library statistics
`cohort_summary_SBS96_uncorrected.txt`	96-context raw mutation counts across all samples (SigProfiler-compatible format)
`cohort_summary_SBS96_corrected.txt`	96-context trinucleotide-corrected counts across all samples
`cohort_summary_SBS96_genome.txt`	96-context estimated mutations per genome across all samples

All SBS96 file can be directly input into SigProfilerPlotting

Parameters

Short	Long	Description
-i	--input	One or more sample folders (space-separated)
-o	--output	Output filename for the summary table (e.g.`cohort_summary.txt`)

Optional: Aggregate Error Profiles Across Samples

DupCaller.py aggregate pools the learned error profiles of several call runs into one set, which can then be passed to call -E (for example, for a low-coverage sample from the same library prep):

DupCaller.py aggregate -i sample1 sample2 sample3 -o pooled
DupCaller.py call ... -E pooled

It sums each sample's ERROR/{sample}.*.txt count tables and re-fits the SRD amplification matrix from the summed base-quality histograms in each sample's tmp/{sample}.amp.tn.bqhist.npz, so the tmp/ folder of each input run must still exist. It writes pooled.amp.tn.txt, pooled.amp.tn.srd.txt, pooled.dmg.tn.txt, and the four pooled.{amp,dmg}.{hp,str}.txt files.

Short	Long	Description	Default
-i	--input	One or more`call` output directories	None
-f	--input-file	File with one error-file prefix per line (e.g.`sample1/ERROR/sample1`), as an alternative or addition to `-i`	None
-o	--output	Output prefix for the aggregated error files (required)	None
-a	--pseudocount	Pseudocount for the SRD matrix re-fit (same as`call -a`)	0.5

Results

For detailed column-by-column and field-by-field descriptions of every output file, see the docs/ folder:

docs/call_outputs.md — all files from DupCaller.py call
docs/estimate_outputs.md — all files from DupCaller.py estimate
docs/summarize_outputs.md — the cohort summary table from DupCaller.py summarize

Since the DBS-by-duplex-group and burden overhaul, SBS/indel/DBS-specific outputs (VCFs, corrected-context tables, burden files, and signature plots) are written into per-type SBS/, INDEL/, and DBS/ subfolders under the sample's output directory, rather than at the sample's top level.

Core Output Files

File	Description
`SBS/{sample}_sbs.vcf`	VCF of detected SNVs and MNVs
`INDEL/{sample}_indel.vcf`	VCF of detected short indel mutations
`DBS/{sample}_dbs.vcf`	VCF of detected dinucleotide substitution (DBS) mutations
`{sample}_coverage.bed.gz`	Duplex coverage depths across genomic positions. For multi-threaded runs, files from different threads are automatically merged with tabix indexing
`{sample}_coverage.bed.gz.tbi`	Tabix index for the coverage BED file
`SBS/{sample}_trinuc_by_duplex_group.txt`	Trinucleotide context counts grouped by duplex read number, used for SBS burden estimation
`INDEL/{sample}_indel_by_duplex_group.txt`	Indel context counts grouped by duplex read number, used for indel burden estimation
`DBS/{sample}_dbs_by_duplex_group.txt`	144-class dinucleotide context counts grouped by duplex read number, used for DBS burden estimation
`{sample}_duplex_family_strand_composition.txt`	Strand composition statistics for duplex read families
`{sample}_duplex_family_strand_composition_heatmap.pdf`	Heatmap visualization of duplex family strand composition
`{sample}_call_params.log`	Full record of all resolved parameters used for the`call` run
`{sample}_stats.txt`	Overall sequencing and analysis metrics

Error Profile Files

Written to the ERROR/ subfolder of the call output directory. If a complete set of the six files used for calling already exists there, learning is skipped and they are reused.

File	Description
`ERROR/{sample}.amp.tn.txt`	Raw amplification SBS mismatch profile by trinucleotide context (diagnostic only)
`ERROR/{sample}.amp.tn.srd.txt`	SRD-EM-fitted amplification SBS error rate matrix actually used for calling
`ERROR/{sample}.dmg.tn.txt`	Damage SBS error profile by trinucleotide context
`ERROR/{sample}.amp.hp.txt`	Amplification indel error rates for homopolymer contexts
`ERROR/{sample}.amp.str.txt`	Amplification indel error rates for short-tandem-repeat contexts
`ERROR/{sample}.dmg.hp.txt`	Damage indel error rates for homopolymer contexts
`ERROR/{sample}.dmg.str.txt`	Damage indel error rates for short-tandem-repeat contexts

Burden Estimation Files

File	Description
`SBS/{sample}_sbs_burden.txt`	SBS burden with uncorrected and corrected estimates and 95% confidence intervals
`INDEL/{sample}_indel_burden.txt`	Indel burden with 95% confidence intervals, including masked and unmasked calculations
`DBS/{sample}_dbs_burden.txt`	DBS burden with uncorrected and corrected estimates and 95% confidence intervals
`SBS/{sample}_sbs_96_corrected.txt`	Corrected SBS counts across 96 trinucleotide contexts for signature analysis, including a`mutations_per_opportunity` column
`INDEL/{sample}_indel_83_corrected.txt`	Corrected indel counts across the 83 ID contexts, including a`mutations_per_opportunity` column
`DBS/{sample}_dbs_78_corrected.txt`	Corrected DBS counts across the 78 DBS contexts, including a`mutations_per_opportunity` column
`SBS/{sample}_sbs_burden_by_group_size.txt`	SBS burden stratified by duplex group size, both cumulative ("min group size >= N") and exact ("group size == N"), corrected and uncorrected, with 95% CI
`INDEL/{sample}_indel_burden_by_group_size.txt`	Same stratification as above, for indels
`DBS/{sample}_dbs_burden_by_group_size.txt`	Same stratification as above, for DBS
`{sample}_duplex_allele_counts.txt`	Duplex depths and allele counts for each unique mutation
`{sample}_estimate_params.log`	Full record of all resolved parameters used for the`estimate` run

Mutation number per genome (and its 95% bounds) in each _burden.txt is the corrected burden multiplied by Reference base number, the number of reference bases considered (one copy of each region in -r). It is therefore the expected number of mutations per haploid genome; double it for a diploid genome. The same applies to the per-channel mutation_number_genome column of the _corrected.txt files and to summarize's *_mutations_per_genome columns.

Visualization Files

File	Description
`SBS/SBS_96_plots_{sample}.pdf`	96-context SBS mutational signature plots (uncorrected and corrected counts)
`INDEL/ID_83_plots_{sample}.pdf`	83-context indel (ID83) mutational signature plots
`DBS/DBS_78_plots_{sample}.pdf`	78-context DBS mutational signature plots
`SBS/{sample}_sbs_burden_by_group_size.pdf`	SBS burden across cumulative and exact duplex group sizes
`INDEL/{sample}_indel_burden_by_group_size.pdf`	Same, for indels
`DBS/{sample}_dbs_burden_by_group_size.pdf`	Same, for DBS

Optional / Conditional Files

File	Condition	Description
`{sample}_gene_coverage.txt`	`-gb` option	Mean duplex coverage per gene for dNdScv correction
`SBS/{sample}_sbs_burden_re_estimate.txt`	`-rb` option	Re-estimated SBS burden for specific regions
`INDEL/{sample}_indel_burden_re_estimate.txt`	`-rb` option	Re-estimated indel burden for specific regions
`SBS/{sample}_sbs_96_corrected_re_estimate.txt`	`-rb` option	Re-estimated 96-context SBS counts for specific regions
`SBS/SBS_96_plots_{sample}_re_estimate.pdf`	`-rb` option	Signature plots for re-estimated regions
`SBS/{sample}_sbs_flt.vcf`	`-d` option	SNV calls kept after the`--dilute` normal-comparison filter

End-to-End Examples

The following examples assume the reference genome has already been indexed (DupCaller.py index, samtools faidx, bwa-mem2 index; Step 1) and that the germline VCF and noise mask are available. Replace all placeholder paths with your actual files.

With Matched Normal

SAMPLE=sample1
NORMAL=normal1
REF=/path/to/hg38.fa
GERMLINE=/path/to/gnomad.hg38.vcf.gz
NOISE=/path/to/noise_mask.bed.gz
THREADS=16

# 1. Trim barcodes
DupCaller.py trim -i ${SAMPLE}_R1.fq.gz -i2 ${SAMPLE}_R2.fq.gz -p NNNXXXX -o ${SAMPLE}
DupCaller.py trim -i ${NORMAL}_R1.fq.gz -i2 ${NORMAL}_R2.fq.gz -p NNNXXXX -o ${NORMAL}

# 2. Align
bwa-mem2 mem -C -T 0 -t ${THREADS} -R "@RG\tID:${SAMPLE}\tSM:${SAMPLE}\tPL:ILLUMINA" \
    ${REF} ${SAMPLE}_1.fastq ${SAMPLE}_2.fastq | samtools sort -@ ${THREADS} > ${SAMPLE}.bam
samtools index -@ ${THREADS} ${SAMPLE}.bam

bwa-mem2 mem -C -T 0 -t ${THREADS} -R "@RG\tID:${NORMAL}\tSM:${NORMAL}\tPL:ILLUMINA" \
    ${REF} ${NORMAL}_1.fastq ${NORMAL}_2.fastq | samtools sort -@ ${THREADS} > ${NORMAL}.bam
samtools index -@ ${THREADS} ${NORMAL}.bam

# 3. Mark duplicates
gatk MarkDuplicates -I ${SAMPLE}.bam -O ${SAMPLE}.mkdped.bam -M ${SAMPLE}.mkdp_metrics.txt \
    --READ_NAME_REGEX "(?:.*:)?([0-9]+)[^:]*:([0-9]+)[^:]*:([0-9]+)[^:]*$" \
    --DUPLEX_UMI --TAGGING_POLICY OpticalOnly --BARCODE_TAG DB
samtools index ${SAMPLE}.mkdped.bam

gatk MarkDuplicates -I ${NORMAL}.bam -O ${NORMAL}.mkdped.bam -M ${NORMAL}.mkdp_metrics.txt \
    --READ_NAME_REGEX "(?:.*:)?([0-9]+)[^:]*:([0-9]+)[^:]*:([0-9]+)[^:]*$" \
    --DUPLEX_UMI --TAGGING_POLICY OpticalOnly --BARCODE_TAG DB
samtools index ${NORMAL}.mkdped.bam

# 4. Call variants
DupCaller.py call -b ${SAMPLE}.mkdped.bam -f ${REF} -o ${SAMPLE} \
    -n ${NORMAL}.mkdped.bam -g ${GERMLINE} -m ${NOISE} -p ${THREADS}

# 5. Estimate mutational burden
DupCaller.py estimate -i ${SAMPLE} -f ${REF} -r chr{1..22} chrX

Without Matched Normal

When no matched normal is available, omit -n and set -maf to cap the maximum allele frequency of called variants. A value of 0.1 is appropriate for most somatic applications to exclude common germline variants.

SAMPLE=sample1
REF=/path/to/hg38.fa
GERMLINE=/path/to/gnomad.hg38.vcf.gz
NOISE=/path/to/noise_mask.bed.gz
THREADS=16

# 1. Trim barcodes
DupCaller.py trim -i ${SAMPLE}_R1.fq.gz -i2 ${SAMPLE}_R2.fq.gz -p NNNXXXX -o ${SAMPLE}

# 2. Align
bwa-mem2 mem -C -T 0 -t ${THREADS} -R "@RG\tID:${SAMPLE}\tSM:${SAMPLE}\tPL:ILLUMINA" \
    ${REF} ${SAMPLE}_1.fastq ${SAMPLE}_2.fastq | samtools sort -@ ${THREADS} > ${SAMPLE}.bam
samtools index -@ ${THREADS} ${SAMPLE}.bam

# 3. Mark duplicates
gatk MarkDuplicates -I ${SAMPLE}.bam -O ${SAMPLE}.mkdped.bam -M ${SAMPLE}.mkdp_metrics.txt \
    --READ_NAME_REGEX "(?:.*:)?([0-9]+)[^:]*:([0-9]+)[^:]*:([0-9]+)[^:]*$" \
    --DUPLEX_UMI --TAGGING_POLICY OpticalOnly --BARCODE_TAG DB
samtools index ${SAMPLE}.mkdped.bam

# 4. Call variants (no matched normal; filter by maximum allele frequency)
DupCaller.py call -b ${SAMPLE}.mkdped.bam -f ${REF} -o ${SAMPLE} \
    -g ${GERMLINE} -m ${NOISE} -maf 0.1 -p ${THREADS}

# 5. Estimate mutational burden
DupCaller.py estimate -i ${SAMPLE} -f ${REF} -r chr{1..22} chrX

Resources

Pre-built reference indexes, germline VCFs, and noise masks are available for download:

Human (GRCh38/hg38)

Resource	Link/Source
Reference genome	TCGA hg38 reference file
Reference index	Pre-built hg38 DupCaller reference
Germline VCF	af-only-gnomad.hg38.vcf.gz file from the legacy GATK resource bundle. A copy of the file can be foundhere.
Noise mask	NanoSeq noise mask

Mouse (GRCm39/mm39)

Resource	Link
Reference genome	UCSC mm39 reference file
Reference index	Pre-built mm39 DupCaller reference
Germline VCF	mgp_strains — includes VCF for SNPs in popular mouse strains. Use the vcf file for your strain.
Noise mask	NOISE.mm39.bed.gz, in-house noise mask from mouse duplex sequencing data

Citation

DupCaller: Cheng, Y. et al. Improved Mutation Detection in Duplex Sequencing Data with Sample-Specific Error Profiles. bioRxiv (2025). https://doi.org/10.1101/2025.07.13.664565

PERF: Avvaru, A. K., Sowpati, D. T. & Mishra, R. K. PERF: an exhaustive algorithm for ultra-fast and efficient identification of microsatellites from large DNA sequences. Bioinformatics 34, 943-948 (2018). https://doi.org/10.1093/bioinformatics/btx721

Copyright

BSD 2-Clause License

Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:

Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

Contact

Yuhe Cheng (yuc211@ucsd.edu)

Name		Name	Last commit message	Last commit date
.github/workflows		.github/workflows
docs		docs
nextflow		nextflow
src		src
tests		tests
.gitignore		.gitignore
.pre-commit-config.yaml		.pre-commit-config.yaml
.prettierignore		.prettierignore
CHANGELOG.md		CHANGELOG.md
Dockerfile		Dockerfile
LICENSE		LICENSE
README.md		README.md
setup.py		setup.py

Latest commit

History

Folders and files

Repository files navigation

DupCaller

Table of Contents

Prerequisites

Installation

pip Installation

Docker / Singularity

Pipeline

Step 1: Index Reference Genome

Parameters

Step 2: Trim Barcodes

Step 3: Align Reads

Step 4: Mark Duplicates

Step 5: Call Variants

Parameters

Required

Recommended

Optional

Advanced

Germline and Noise Masks

Low-Coverage Error-Profile Fallback

Strand Independence Filter

Step 6: Estimate Mutational Burden

Per-Gene Coverage

Re-estimation for Specific Regions

Parameters

Required

Optional

Step 7: Summarize Across Samples

Parameters

Optional: Aggregate Error Profiles Across Samples

Results

Core Output Files

Error Profile Files

Burden Estimation Files

Visualization Files

Optional / Conditional Files

End-to-End Examples

With Matched Normal

Without Matched Normal

Resources

Human (GRCh38/hg38)

Mouse (GRCm39/mm39)

Citation

Copyright

Contact

About

Resources

Stars

Watchers

Forks

Releases

Packages

Used by

Contributors

Languages