BAMulator — Fast variant simulation on BAMs

Simulate variants directly on existing BAM files

Install

pip3 install -r requirements.txt

Usage

Usage: bamulator [OPTIONS]

Options:
  -v, --variants TEXT     Tab-separated file listing BAM files and variants
                          [required]
  -t, --threads INTEGER   Total number of threads
  -r, --reference TEXT    Genome reference in FASTA format  [required]
  -o, --output TEXT       Output directory for simulated BAMs and phased VCFs
                          [required]
  -s, --suffix TEXT       Suffix for simulated BAM files
  --input_vcf_list TEXT   File listing input VCFs or directory containing VCFs
  --genome_version TEXT   Genome version for phasing (e.g., hg38, hg19)
  --map_dir TEXT          Directory containing PLINK .map files for phasing
  --ref_dir TEXT          Directory containing .bref3 reference files for
                          phasing
  --help                  Show this message and exit.

Quick start

bamulator \
  --variants variants.txt \
  --reference /path/to/ref.fasta \
  --threads 4 \
  --suffix .simulated \
  --output /path/to/outdir

Output BAMs are named <sample><suffix>.bam (default: .simulated).

Variant config (variants.txt)

Tab-separated, one variant per line, with a header row. The header is used to detect which of the two formats below a file uses.

SNVs / indels

Header must include SAMPLE and REF.

SAMPLE	CHROM	POS	REF	ALT	TYPE	CLONE	CLONAL_PROPORTION	GENOTYPE	ADDITIONAL_INFO
/home/user/input/sample1.bam	chr7	55259515	T	G	snv	clone1	1.0	0/1	EGFR p.L858R
/home/user/input/sample2.bam	chr7	55242467	GAATTAAGAGAAGCAACA	GTTGCT	delins	clone1	1.0	0/1	.

• POS is 1-based
• TYPE: snv, insertion, deletion, or delins
• CLONAL_PROPORTION: tumor purity / VAF (0-1)
• GENOTYPE: 0/1 (heterozygous) or 1/1 (homozygous) -- combined with CLONAL_PROPORTION to set the simulated allele fraction

Structural variants (DEL / DUP / INV / translocations)

Header must include SVTYPE.

SAMPLE	CHR1	START1	END1	CHR2	START2	END2	SVTYPE	SVID	PLOIDY	SV_MECHANISM	NTINS	HOMOLOGY	CLONE	CLONAL_PROPORTION	ADDITIONAL_INFO	HAPLOTYPE
/home/user/input/sample3.bam	chr1	156082722	156089883	.	.	.	DEL	SOS1_7_to_8	1	random	.	.	clone1	1.0	.	1
/home/user/input/sample4.bam	chr7	150652435	150661275	.	.	.	DUP	TAZ_1_to_4	3	random	.	.	clone1	1.0	.	1
/home/user/input/sample5.bam	chr2	29445267	42519312	.	.	.	INV	EML4_ALK_inv	2	random	.	.	clone1	1.0	.	1
/home/user/input/sample6.bam	chr4	25668400	25668400	chr6	117645600	117645600	BALANCED_TRANSLOCATION	SLC34A2_ROS1	2	random	.	.	clone1	1.0	.	1

• SVTYPE: DEL, DUP, INV, BALANCED_TRANSLOCATION, or WHOLE_ARM_TRANSLOCATION
• CHR2/START2/END2: partner breakpoint locus, . for single-locus SVs (DEL/DUP/INV); required for the two translocation types
• PLOIDY: copy number at the locus after the event, against a diploid (CN=2) background -- e.g. 1 = heterozygous deletion, 0 = homozygous deletion, 3 = single-copy gain, 4 = two-copy gain. Combined with CLONAL_PROPORTION (tumor purity), this sets how far simulated depth moves from baseline for DEL/DUP. INV/translocations are balanced rearrangements (no depth change), so use 2 here -- their penetrance (what fraction of reads show the junction) is instead driven by CLONAL_PROPORTION together with HAPLOTYPE cardinality (see below).
• NTINS: non-template bases inserted at the breakpoint junction (e.g. from NHEJ repair), . for none
• HAPLOTYPE: which copy carries the variant -- a single value (0, 1, or 2) means only one chromosome copy does (heterozygous-equivalent); a comma-separated list (e.g. 1,2) means every copy does. For INV/ translocations this is combined with CLONAL_PROPORTION the same way GENOTYPE zygosity combines with VAF for SNVs (DEL/DUP don't use it this way -- PLOIDY already encodes their zygosity).

Combining SV + SNV/indel variants

A single config file can't mix both formats (see format detection above). To get both kinds of variant into one output BAM, run bamulator twice, pointing the second run's SAMPLE column at the first run's <sample><suffix>_merged.bam instead of the original input BAM:

bamulator -v svs.txt -r ref.fasta -o out1  # produces out1/sample.simulated_merged.bam
# edit snvs.txt's SAMPLE column to out1/sample.simulated_merged.bam
bamulator -v snvs.txt -r ref.fasta -o out2  # final BAM has both SV and SNV/indel variants

Phasing (optional)

Provide --input_vcf_list plus:

• --genome_version (e.g., hg38, hg19)
• --map_dir with PLINK .map
• --ref_dir with .bref3

Phased VCFs are written to --output. For DEL variants, nearby phased heterozygous SNVs are also painted onto the simulated reads (REF or ALT depending on which haplotype the deletion's HAPLOTYPE matches), so the deletion sits in a realistic phased-heterozygosity background rather than a suspiciously clean flanking region.

Features

• SNVs, indels (including complex delins)
• Structural variants: deletions, duplications, inversions, and balanced / whole-arm translocations, with breakpoint-junction reads (including non-template insertions, NTINS)
• Tumor purity (CLONAL_PROPORTION) drives simulated allele fraction throughout: with GENOTYPE for SNVs/indels, with PLOIDY for DEL/DUP depth, and with HAPLOTYPE cardinality for INV/translocation penetrance
• Phased background heterozygous SNVs are painted onto reads spanning a simulated deletion, for realistic surrounding haplotype context
• Multithreaded execution
• Every run writes <sample>_coverage_warnings.tsv to --output, flagging any variant whose target locus had zero or too few (<10) overlapping read pairs in the source BAM -- such a variant is silently absent from the output otherwise, indistinguishable from "successfully simulated but invisible". Check this file when a configured variant doesn't show up.

Limitations

• No SVs requiring sequenced breakpoints as input (coordinates only)
• HOMOLOGY is parsed but not yet applied to simulated sequences

Testing

python3 -m pytest tests/

An integration suite that builds a small synthetic reference + BAM once per session and drives the real bamulator binary plus the bwa/samtools post-processing pipeline for every variant type, checking coverage, breakpoint sequences, and NTINS content against ground truth rather than mocking anything. Takes well under a minute; no fixture files are checked in.

Presets

scripts/presets/ has ready-made variant configs for common tumor hotspot panels, one set of files per tumor type. A single variant config can't mix SNV/indel and SV rows (see format detection above), and some hotspots occupy the same position as another and so can't be run together -- those are split into a companion *_alternates.tsv file:

• NSCLC: lung_cancer_small_variants.tsv (EGFR, KRAS, BRAF, MET, ALK/ROS1/RET resistance mutations, TP53, etc.), lung_cancer_SV.tsv (EML4::ALK, KIF5B::RET, CCDC6::RET, *::ROS1 fusions; CDKN2A/B deletions), lung_cancer_small_variants_alternates.tsv (RET p.G810R vs p.G810S, PTEN p.R130Q vs p.R130fs)
• Melanoma: melanoma_small_variants.tsv (BRAF, NRAS, KIT, NF1), melanoma_small_variants_alternates.tsv (BRAF p.V600E vs p.V600K)
• Colorectal: colorectal_small_variants.tsv (KRAS, NRAS, BRAF, PIK3CA, TP53, APC)
• Breast: breast_cancer_small_variants.tsv (PIK3CA, TP53, ESR1, BRCA1/2, CDH1), breast_cancer_SV.tsv (ERBB2/HER2 amplification)
• AML: aml_small_variants.tsv (NPM1, FLT3-ITD, IDH1/2, DNMT3A, KIT, TP53)
• Cerebral (glioma/GBM): cerebral_tumor_small_variants.tsv (IDH1/2, TERT promoter, TP53, PTEN, BRAF, H3F3A K27M/G34R, ATRX), cerebral_tumor_SV.tsv (CDKN2A/B deletions, EGFR amplification, EGFRvIII exon2-7 deletion, KIAA1549::BRAF fusion)

All genomic coordinates are GRCh37/hg19, individually verified against the reference sequence (not just copied from a table) -- codon-level where a plain REF-base match could coincide with the wrong position. A few entries represent a real, verified example variant standing in for a gene whose pathogenic mutations are diffuse rather than a single hotspot (e.g. NF1, CDH1, FLT3-ITD) -- noted in each row's ADDITIONAL_INFO. Edit the SAMPLE column to point at your BAM(s) before running. To combine SV and small-variant hotspots in one output BAM, chain two runs -- see "Combining SV + SNV/indel variants" below.

hg38

scripts/presets/hg38/ mirrors every file above, lifted from hg19 with liftOver (UCSC hg19ToHg38.over.chain) and then independently re-verified -- every lifted position's reference base was re-fetched from the real GRCh38 sequence and compared against the hg19 REF (all matched with zero mismatches for this panel; any that hadn't would have needed a manual re-derivation instead of a straight coordinate carry-over). The three insertions built by duplicating adjacent reference sequence (NPM1 p.W288fs, FLT3-ITD, ERBB2 p.Y772_A775dup) had their duplicated segment independently re-derived from GRCh38 too, not just carried over. Use these files with --reference pointed at a GRCh38 FASTA.

Helper scripts

Generate structural variant configs (DEL/DUP/INV/translocations) over a gene panel. Set any rate to 0 to exclude that type -- e.g. for CNVs only:

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

Generate SNV/indel configs:

python3 scripts/generate_indels_config.py \
  --indir /path/to/bam_dir \
  --bed /path/to/gene_panel.bed \
  --genome /path/to/ref.fasta \
  --num_variants 20 \
  --proportions 1 \
  --output variants.tsv