BAMulator

BAMulator simulates variants directly into existing BAM files while preserving the characteristics of the original sequencing data.

Supported variant types include:

• SNVs
• INDELs
• CNVs
• SVs (including balanced and whole-arm translocations)
• Subclonal variants
• Haplotype-based variants

Requirements

Input BAMs must be coordinate sorted and indexed.

Install the Python dependencies and build the simulation engine:

python3 -m pip install -r requirements.txt
make -C src

The wrapper expects Linux-compatible binaries at:

bwa/bwa
samtools/samtools

Optional phasing also requires java, bcftools, bgzip, and tabix on PATH.

Quick start

Variants can be provided in three ways:

1. Use an existing VCF, BCF or BEDPE
2. Generate a BAMulator config using the included scripts
3. Write a BAMulator TSV manually

For most use cases, the first two options are recommended.

Use an existing VCF

python3 bamulator.py \
  --variants variants.vcf.gz \
  --sample_bam sample.bam \
  --reference hg38.fa \
  --output simulated

The main output is:

simulated/sample_merged.bam

VCF, compressed VCF, BCF and BEDPE inputs are supported.

For multi-sample VCF/BCF files, use --sample_map to associate each sample with its BAM.

Generate variant configurations

BAMulator includes scripts to automatically generate variants from a BAM directory and a gene-panel BED.

Generate SNVs and indels

For example, generate 20 small variants:

python3 scripts/generate_indels_config.py \
  --indir /path/to/bams \
  --bed panel.bed \
  --genome hg38.fa \
  --num_variants 20 \
  --proportions 1 \
  --output variants.tsv

Then simulate them:

python3 bamulator.py \
  --variants variants.tsv \
  --reference hg38.fa \
  --output simulated

Exactly one of --indir or --list must be provided.

The SNV, insertion, deletion and delins rates must sum to 1.

Generate structural variants

For example, generate a set containing only deletions and duplications:

python3 scripts/generate_sv_config.py \
  --indir /path/to/bams \
  --bed panel.bed \
  --num_variants 20 \
  --proportions 1 \
  --deletion_rate 0.5 \
  --duplication_rate 0.5 \
  --inversion_rate 0 \
  --translocation_rate 0 \
  --output variants.tsv

Then run:

python3 bamulator.py \
  --variants variants.tsv \
  --reference hg38.fa \
  --output simulated

Non-zero SV rates are automatically normalized.

Manually define variants

Manual TSV files are useful when you want to simulate specific variants.

Small variant example

SAMPLE	CHROM	POS	REF	ALT	TYPE	CLONE	CLONAL_PROPORTION	GENOTYPE	ADDITIONAL_INFO
sample.bam	chr7	55259515	T	G	snv	clone1	1.0	0/1	EGFR_p.L858R

Run:

python3 bamulator.py \
  --variants variants.tsv \
  --reference hg19.fa \
  --output simulated

Supported small-variant types are:

snv
insertion
deletion
delins

POS is 1-based.

CLONAL_PROPORTION controls the fraction of the sample carrying the variant. For example:

1.0   clonal
0.5   50% of the sample
0.2   20% of the sample

An optional HAPLOTYPE column can be used for explicit allele-specific simulation.

For a small-variant pass run on an ASCN output BAM, add mutation_multiplicity=K to ADDITIONAL_INFO to mutate only the first K descendant copies of the selected parental haplotype. For example, a mutation with multiplicity 2 on the major allele of a 4+2 segment is applied to two of the four hap1 copies. K must not exceed the copy count of the selected haplotype. If mutation_multiplicity is omitted, all descendants of that haplotype retain the historical behavior and receive the mutation.

When ASCN and small variants share a CLONE, stable fragment selection makes the second pass a nested subset of the first. For tumor purity p, cancer-cell fraction f, total tumor copy number C, and mutation multiplicity K, use CLONAL_PROPORTION = p * f; the expected bulk-sample VAF is p * f * K / (p * C + 2 * (1 - p)).

Structural variant example

Structural variants use a different TSV format.

Example heterozygous deletion:

SAMPLE	CHR1	START1	END1	CHR2	START2	END2	SVTYPE	SVID	PLOIDY	SV_MECHANISM	NTINS	HOMOLOGY	CLONE	CLONAL_PROPORTION	ADDITIONAL_INFO	HAPLOTYPE
sample.bam	chr9	21970277	21975386	.	.	.	DEL	CDKN2A_del	1	random	.	.	clone1	1.0	.	1

Supported structural variant types are:

DEL
DUP
ASCN
INV
BALANCED_TRANSLOCATION
WHOLE_ARM_TRANSLOCATION

Common PLOIDY values:

PLOIDY	Meaning
0	Homozygous deletion
1	Heterozygous deletion
2	Balanced event
3	One-copy gain
4	Two-copy gain

Native TSV breakpoint coordinates are 1-based.

ASCN represents an allele-specific copy-number segment without asserting a particular breakpoint mechanism. Set PLOIDY to total copy number and provide the target parental copy counts in ADDITIONAL_INFO:

SAMPLE  CHR1  START1  END1  ...  SVTYPE  SVID       PLOIDY  ...  ADDITIONAL_INFO             HAPLOTYPE
sample.bam chr17 1 20000000 ...  ASCN    chr17_loh  2       ...  major_cn=2;minor_cn=0       hap1

The engine assigns every fragment consistently to haplotype 1 or 2, treats haplotype 1 as the major allele, and emits the configured number of copies in the CLONAL_PROPORTION fraction. This supports CN-neutral LOH (2+0), balanced whole-genome duplication (2+2), and asymmetric amplification such as 4+2. ASCN models read depth and phased B-allele frequency; it does not invent tandem-junction reads for copy-number segments inferred without a known breakpoint structure.

When every row for a sample is ASCN and its segments do not overlap, BAMulator automatically uses a direct-BAM implementation. It traverses the coordinate-sorted input once, applies phased background SNVs in place, and drops or duplicates aligned read pairs without changing their coordinates. With --threads N greater than one, the reader feeds bounded 4,096-record chunks to N-1 simulation workers while the coordinating thread preserves the original record order during output. This parallelizes record decoding, copy selection, phasing, and editing without chromosome seeks, temporary BAM shards, or unbounded memory. --threads 1 retains the lower-overhead serial path.

The threaded implementation deliberately has only two roles. The coordinating thread is the only thread that reads and writes BAM data. Worker threads take complete chunks, transform their records, and mark the chunks ready. The coordinator can read the next chunk while workers process earlier ones, but it always waits for and writes the oldest chunk first. A small bounded queue (at most two chunks per worker) prevents a fast reader from consuming unbounded memory. The queue mutex protects only the short act of handing out a chunk; the expensive record transformation happens outside that lock.

The final indexed <sample>_merged.bam is written directly: no uncompressed FASTQs, untouched-read BAM, BWA realignment, or genome-wide read-name table is created. Mixed ASCN/non-ASCN inputs and overlapping ASCN segments retain the general FASTQ/remapping path. Use --legacy_ascn only to force that older path for comparison or troubleshooting.

For point breakpoints such as translocations:

START1 = END1
START2 = END2

Disease presets

Ready-made variant configurations, curated from published hotspot and dosage-sensitivity data (COSMIC, ClinVar, ClinGen, OncoKB, and the primary literature), are available under:

scripts/presets/hg19/
scripts/presets/hg38/

Each build is organized by disease category:

scripts/presets/<build>/solid_tumors/
scripts/presets/<build>/hematological/
scripts/presets/<build>/germline/

For example:

python3 bamulator.py \
  --variants scripts/presets/hg19/solid_tumors/lung_cancer_small_variants.tsv \
  --reference hg19.fa \
  --output simulated

Edit the SAMPLE column in the preset before running.

Available preset groups include:

• Solid tumors — NSCLC, melanoma, colorectal cancer, breast cancer, CNS / GBM, pancreatic cancer, prostate cancer, ovarian cancer, gastric cancer
• Hematological — AML
• Germline — recurrent syndromic microdeletion/microduplication rearrangements (neurological and cardiac), e.g. 22q11.2 (DiGeorge), 7q11.23 (Williams-Beuren), 17p12 (CMT1A/HNPP), 17p11.2 (Smith-Magenis/Potocki-Lupski), 17q11.2 (NF1 microdeletion), 15q11-q13 (Prader-Willi/Angelman), 1p36, 5p15 (Cri-du-chat); germline small-variant coverage not yet included

Some panels have separate files for small variants and structural variants.

Variants that cannot be simulated together because they occupy the same position are stored in companion *_alternates.tsv files.

Each preset row carries a SOURCE column citing the database or publication backing that variant (e.g. COSMIC, ClinVar, PMID:xxxxxxxx), so entries can be traced back to published evidence rather than taken on faith.

Combining small and structural variants

Small variants and structural variants cannot currently be simulated in the same BAMulator run.

To simulate both, run BAMulator twice.

First simulate the structural variants:

python3 bamulator.py \
  --variants svs.tsv \
  --reference hg38.fa \
  --output sv_out

Then use the resulting BAM as the input for the small-variant configuration:

sv_out/sample_merged.bam

and run BAMulator again:

python3 bamulator.py \
  --variants snvs.tsv \
  --reference hg38.fa \
  --output final

Reproducible simulations

Use --seed to reproduce the same stochastic simulation:

python3 bamulator.py \
  --variants variants.tsv \
  --reference hg38.fa \
  --output simulated \
  --seed 42

The effective seed is logged for every run.

Outputs

A typical run produces:

sample_merged.bam
sample_merged.bam.bai

sample_simulated.bam
sample_simulated.bam.bai

sample_unsimulated.bam

sample_simulated_R1.fq
sample_simulated_R2.fq

sample_coverage_warnings.tsv

bamulator.truth.vcf

The main downstream file is:

sample_merged.bam

bamulator.truth.vcf contains the variants requested for the simulation.

sample_coverage_warnings.tsv reports variants with insufficient source coverage.

Check this file whenever a configured variant does not appear in the output.

Variants with zero or fewer than 10 overlapping read pairs are reported as coverage warnings.

Phasing

Optional phasing can be enabled with:

--phase_snvs

Example:

python3 bamulator.py \
  --variants deletions.tsv \
  --reference hg38.fa \
  --output simulated \
  --phase_snvs \
  --input_vcf_list vcfs/ \
  --genome_version hg38 \
  --map_dir maps/ \
  --ref_dir reference-panels/

Phasing allows nearby heterozygous SNVs to be retained on consistent parental haplotypes.

Supported genome labels are:

hg19
b37
hg38
b38

The optional HAPLOTYPE field can explicitly assign variants to:

maternal
paternal
both

Nearby variants assigned to the same parental haplotype remain in cis, while variants assigned to opposite haplotypes remain in trans.

Docker variant-caller runner

The variant-caller pipeline now lives at /media/bdelolmo/BERNAT/BAMulator_paper/run_variant_callers.py. The legacy scripts/run_variant_callers.py path remains as a compatibility launcher. It provides one Docker-based interface for:

• GATK Mutect2 (tumor-only or tumor/normal, orientation-bias learning, optional contamination estimation, and FilterMutectCalls)
• GATK HaplotypeCaller (gVCF generation, genotyping, and standard SNP/indel hard filters)
• Manta (germline or paired somatic SV calling)
• Strelka2 (germline or paired somatic SNV/indel calling)
• DeepVariant (germline small-variant calling)
• GRIDSS (germline or paired somatic SV calling, with gridss_somatic_filter in somatic mode)
• CNVkit (WGS, WES, or targeted-panel CNV calling; matched normal, reusable reference, or automatic flat-reference tumor-only mode)
• GATK ModelSegments/CallCopyRatioSegments (WGS, WES, or panel somatic CNA calling; optional read-count PoN and optional tumor/normal BAF modeling)
• CNVpytor (single-sample WGS read-depth calling with optional phased-SNP/BAF analysis and no control cohort)

Each caller writes to its own subdirectory under --output-dir. The script checks FASTA/BAM indexes, records the exact commands, captures one log per step, and refuses to overwrite a non-empty result. The layout is consistent across callers:

caller_results/
└── CALLER/
    ├── results/                 final VCFs and indexes
    ├── work/                    intermediate files and native workflows
    ├── logs/                    one log per pipeline step
    └── manifest/
        └── run_manifest.json    image, commands, and layout

Inputs are bind-mounted read-only, except that GRIDSS needs write access to the reference directory for its standard setupreference cache/index step. Use --dry-run to inspect all commands before running them.

Use --only to select one or more comma-separated callers from a shared pipeline invocation. Common options are passed to each selected caller, while caller-specific options are routed only to callers that support them:

python3 scripts/run_variant_callers.py --only cnvpytor,cnvkit \
  --bam germline.bam --tumor-bam tumor.bam \
  --chromosome chr22 --sequencing-mode WGS \
  --reference reference.fa --output-dir caller_results --threads 8

# Somatic SNVs/indels
python3 scripts/run_variant_callers.py mutect2 \
  --tumor-bam tumor.bam --normal-bam normal.bam \
  --tumor-sample TUMOR --normal-sample NORMAL \
  --reference reference.fa --output-dir caller_results --threads 8

# Paired somatic SVs
python3 scripts/run_variant_callers.py manta \
  --tumor-bam tumor.bam --normal-bam normal.bam \
  --reference reference.fa --output-dir caller_results --threads 8

# Paired somatic SVs with GRIDSS
python3 scripts/run_variant_callers.py gridss \
  --tumor-bam tumor.bam --normal-bam normal.bam \
  --tumor-sample TUMOR --normal-sample NORMAL \
  --reference reference.fa --output-dir caller_results --threads 8

# Germline small variants
python3 scripts/run_variant_callers.py deepvariant \
  --bam sample.bam --model-type WES --intervals targets.bed \
  --reference reference.fa --output-dir caller_results --threads 8

# Tumor-only CNVkit WGS call. With no normal/reference, a flat reference is built.
python3 scripts/run_variant_callers.py cnvkit \
  --tumor-bam simulated.bam --sequencing-mode WGS \
  --reference reference.fa --output-dir caller_results --threads 8

# Tumor-only GATK CNA call, restricted to chr22 and with optional BAF sites.
python3 scripts/run_variant_callers.py gatk-cnv \
  --tumor-bam simulated.chr22.bam --sequencing-mode WGS \
  --intervals chr22.bed --snp-sites common_snps.vcf.gz \
  --reference reference.fa --output-dir caller_results --threads 8

# Single-sample CNVpytor chr22 call with BAM-recounted BAF values.
python3 scripts/run_variant_callers.py cnvpytor \
  --bam simulated.chr22.bam --chromosome chr22 \
  --snp-vcf phased_snps.vcf.gz --sample-id HG002 --use-phase \
  --include-nonpass-snps \
  --reference-genome hg38 --bin-size 1000 --bin-size 10000 \
  --plot \
  --reference reference.fa --output-dir caller_results --threads 8

The default images are broadinstitute/gatk:4.2.2.0 for GATK workflows, bdolmo/manta:1.6.0, aokad/strelka2:latest, google/deepvariant:latest, gridss/gridss:2.13.2, etal/cnvkit:latest, and hydragenetics/cnvpytor:1.3.1; override any default with --image. Manta and Strelka2 require --intervals to be a bgzip-compressed BED with an adjacent .tbi or .csi index. GRIDSS processes the BAM genome-wide and does not accept --intervals. Mutect2 contamination filtering is enabled by providing --contamination-sites with a bgzip/tabix common-sites VCF. Caller-specific options can be appended with repeated --caller-arg=VALUE arguments.

CNVkit and GATK CNV require --intervals in WES/PANEL mode; in WGS mode this argument is optional and can restrict a chromosome-only BAM. Without controls, CNVkit builds a flat copy-neutral reference. GATK uses GC-annotated intervals without a read-count PoN; in that mode its standardized and "denoised" copy ratios are identical because PCA denoising cannot be performed. This is useful for large controlled CNV simulations, but real WES/panel samples will generally have more false segments than runs using process-matched normals. Add --read-count-pon when one becomes available. --snp-sites enables BAF/LOH modeling, and an optional --normal-bam then supplies matched-normal allele counts without requiring a multi-sample PoN.

CNVpytor is the easiest no-cohort baseline for WGS and chromosome-only BAMs. It produces results/cnvpytor.pytor plus a tab-separated results/cnvpytor.calls.tsv. Repeating --bin-size runs several resolutions. When --snp-vcf and --sample-id are supplied, the workflow imports the SNPs, uses -pileup to recalculate allele depths from the simulated BAM, applies the built-in SNP mask, computes BAF tracks, and enables CNVpytor's joint read-depth/BAF caller. Thus it does not rely on stale or placeholder AD values in the input VCF. Site/truth VCFs without AD are supported: positions are imported with -noAD before pileup recounting. Add --plot to write headless PNG plots for every requested bin size; chromosome-restricted plots combine RD, BAF, and BAF-likelihood panels. When phased analysis is requested, both portable unphased plot signals and phased signals/calls are retained. Trusted truth/site VCFs that use FILTER=. instead of PASS require --include-nonpass-snps; avoid that option for raw unfiltered discovery calls. Use an explicit --reference-genome hg19 or hg38 when possible; auto uses CNVpytor's contig-length detection. For phased truth genotypes, --use-phase passes their haplotype orientation into BAF calculation and joint calling. If CNVpytor cannot form eligible BAF normalization bins, the final call automatically falls back to read-depth-only mode and records the effective mode in results/cnvpytor.calls.mode.txt. CNVpytor is not the recommended primary caller for sparse WES or diagnostic-panel data.

Current limitations

• Small variants and structural variants cannot be mixed in one run.
• Mixed small/SV VCF or BCF inputs are rejected.
• Background phasing currently supports SNVs but not background indels.
• HOMOLOGY is parsed but is not currently applied to simulated sequences.
• Phasing currently processes chromosomes 1–22 and X.
• --suffix and --max_reads are accepted for compatibility but currently do not modify the corresponding engine behavior.
• --minins and --maxins are accepted by the SV generator but currently unused. Breakpoint NTINS length is controlled with --ntmin and --ntmax.