forensic-mixture-simulations

Simulation and figure pipeline for the DNAMIXTURE paper (forensic likelihood ratios for short-read sequencing data from DNA mixtures). Every figure in the paper is produced by one Snakemake workflow plus one R script, with fixed RNG seeds throughout, so all results are exactly reproducible (on the same machine and compiler).

The likelihood program itself (DNAMIXTURE) is developed at https://github.com/ras-nielsen/DNAMIXTURE ; the copy in likelihood/ is the version used for the paper. Build it once with cd likelihood && make (the workflows also build it automatically).

Requirements

python3 (pysam, pandas, numpy), snakemake, bcftools, bgzip/tabix, plink 1.9 (LD workflow only), ped-sim (family workflows only), R with tidyverse, viridis, ggh4x, wesanderson, ggforce, ggbreak, patchwork.

Step 0: input VCFs

All inputs derive from the public 1000 Genomes Phase 3 release (GRCh37, 20130502):

snakemake -s paper_scripts/make_vcfs.smk --cores 4

This creates vcf/ with the study-sample subset (biallelic SNPs, EUR_AF in [0.01, 0.99]), a thinned ~500k-site version, and the phased GBR founder VCF for ped-sim.

Step 1: pedigree simulation (family workflows only)

snakemake -s paper_scripts/generate_pedigree.smk --cores 2 --config ped_sim=/path/to/ped-sim

The pedigree (metadata/family.def) encodes the six-generation family of Figure 1 (34 candidate profiles, 10 founders); metadata/victims.def is a separate two-founder pedigree providing the unrelated victim; metadata/relationships.csv maps simulated ids to relationship labels.

Step 2: per-figure simulation workflows

Workflow	Paper figure	Scenario
contam_models.smk	Fig 2/5	contaminant-model (mis)specification x suspect fraction
family_falsepos_L1.smk	Fig 3	L1 for all 34 pedigree candidates
family_models.smk	Fig 4	L1-L4 for all pedigree candidates
populations.smk	Fig 6	ancestry mismatch between suspect and AF panel
error.smk	Fig 7	sequencing-error robustness
depth_snps.smk	Fig 8	depth x SNP-count x mixture-proportion grid
ld.smk	Fig 9	LD pruning x SNP-selection window size
chi.smk	(calibration)	null distribution of L1

Each workflow: simulate read pileups (helpers/generate_pileup.py) -> select SNPs (helpers/selector.py; one SNP with allele frequency closest to 0.5 per 100 kb window, plus one random background SNP per window for the population-match test) -> write DNAMIXTURE JSON (helpers/make_json.py) -> run likelihood/DNAMIXTURE -> parse (helpers/parse_likelihood.py) -> paper_data/<scenario>/results.tsv.

Every simulation receives a deterministic --seed derived from its replicate/condition, recorded in (or reconstructible from) the results table.

Lambda evidence score

DNAMIXTURE computes the lambda evidence score by default on each L1 run: log Lambda = min(L1) over the two contaminant models when max(L2) over the two contaminant models exceeds log(C) (C = 10; -C to change), else 0. The parsed results include the components (l1_pop, l1_ind, l2_pop, l2_ind, gate_passed, log_lambda).

Run DNAMIXTURE without -l to compute only the lambda evidence score; helpers/parse_likelihood.py then reports only the lambda-block fields.

Step 3: figures

Rscript paper_scripts/<scenario>.R

writes the figure(s) to paper_plots/. All log likelihood ratios are natural logarithms (axes use a base-10 pseudo-log transform for display).

Real-case input: BAM + VCF (tools/dnamixture_prep.py)

tools/dnamixture_prep.py builds a DNAMIXTURE input JSON directly from standard formats: a BAM file for the crime-stain reads, VCFs for the suspect and (optionally) the victim, and a sites-only allele-frequency panel (see panels/). Requires Python 3.10+ with pysam (pip install pysam); BAM and VCFs must be indexed and on the same reference build as the panel (the tool verifies contig names and lengths and refuses on mismatch).

python tools/dnamixture_prep.py \
    --bam stain.bam \
    --suspect-vcf suspect.vcf.gz \
    --victim-vcf victim.vcf.gz \
    --panel panels/1000g.phase3.maf01.sites500k.vcf.gz \
    --af-field EUR_AF \
    -o case.json

likelihood/DNAMIXTURE -i case.json -l L1

Sites are selected with the published scheme: per 100 kb window, the biallelic panel SNP with allele frequency closest to 0.5, plus one random background SNP per window for the population-match test. Reads are filtered by mapping and base quality (--min-mq, --min-bq) and indels are skipped. Omit --victim-vcf when no victim genome is available and analyze the JSON with DNAMIXTURE --no-victim.

THE BUNDLED REFERENCE PANELS ARE SUPPLIED FOR TESTING PURPOSES ONLY. FOR FORENSICS APPLICATIONS PLEASE USE A PANEL TAILORED TO YOUR USE CASE (see panels/README.md; tools/make_af_panel.py builds custom panels).