The eidolon project

Formerly rusty-neat / rneat — renamed to eidolon in v2.0.0. Same tool, same NEAT lineage. The deprecated rneat alias was removed in v4.0.0; invoke eidolon. See CHANGELOG.md.

Upgrading from 2.0.0 → 3.0.0? The names of emitted output tokens changed. See Upgrading from 2.0.0 below — read it if you have any script that parses eidolon VCFs or FASTQ/BAM read names.

Welcome to eidolon, a Rust port of NEAT (https://github.com/ncsa/neat), a genetic simulation program that creates fastq that appear to be from sequencers, carry the same statistical properties as your data, and generate a golden bam and fastq that gives you ideal alignments and what variants were inserted. In addition, eidolon generates “noise” in the form of sequencing errors as it is writing out files. These features can help you hone in your alignment and variant calling software to your data. Training models on your data will allow eidolon to faithfully reproduce the statistical properties of your dataset.

We have spent some dedicated time toward gearing the current version of eidolon to simulate cancer genetics, including adding structural variant simulations (CNV, BND, SVs, and others), and creating a wrapper that simulates purity levels and then stitches the results back together. We’ve geared this software with an aim at keeping memory usage as low and CPU time as short as possible. Let us know your real world experience by creating a Feedback issue, if you have something that’s not quite a bug, or have a positive experience to share. As always, let us know if you find a bug and give as many details as you can to help us troubleshoot.

eidolon trains reusable models from your own data (mutation, sequencing-error, fragment-length, GC-bias, and BAM/alignment models), outputs a golden BAM with ideal alignments alongside the FASTQ and truth VCF, accepts a BED to target read creation to regions, and can read custom variants from an input VCF — including per-variant allele frequencies, to reproduce a continuous AF spectrum for pooled or somatic data. More recent releases add native tumor/normal simulation (eidolon gen-cancer-reads), structural variants and copy number, per-tissue somatic models, and trinucleotide-context-aware SNP placement so context-specific mutational signatures reproduce. A file-streaming writer keeps disk I/O and footprint low. See CHANGELOG.md for the full release history, and please open a Feedback issue with your real-world experience.

Find us on Zenodo:

Scope: germline is general, somatic is human

These two halves of eidolon have deliberately different reach, and it is worth knowing which one you are using.

Germline simulation is species-agnostic. gen-reads needs a reference and, optionally, models trained from your own data — nothing in that path assumes a human genome. It has been exercised on H1N1 (13.5 kb, eight contigs), E. coli, yeast, soy, and human references from a single chromosome up to whole-genome GRCh38. If you are simulating germline variation for any organism, that is a supported use.

Somatic / cancer simulation is specialized to human. The bundled mutational models are derived from human cancer cohorts — COSMIC signatures and PCAWG structural-variant distributions — and the per-tissue models (BRCA, skin, lung) are human tissue types. SV size and type distributions, mutation rates, and signature weights all carry that provenance.

gen-cancer-reads will run against a non-human reference and produce output, but the model driving it has no meaning there: a yeast genome given a human pan-cancer SV spectrum produces numbers, not biology. We do not claim cancer simulation for non-human references, and the validation campaigns behind those models are human-only. Use another organism’s reference for cancer work only if you have trained your own models against matched data, and treat the bundled ones as human-specific.

This split is also why the two halves are validated differently: germline correctness is checked across deliberately varied contig shapes and sizes, while somatic realism is measured against real human sequencing data.

Cancer simulation

eidolon simulates tumor / normal sequencing data end-to-end. The eidolon gen-cancer-reads -c <config.yml> subcommand runs two gen-reads passes — one normal-genotype, one tumor — over the same reference and merges them at a configurable purity into a single “tumor biopsy” FASTQ that downstream somatic callers (Mutect2, Strelka, Manta, …) consume directly, plus an origin-tagged truth VCF (INFO/EIDOLON_ORIGIN ∈ {germline, somatic, shared}) for scoring. It also generates the foundational structural-variant types cancer SVs depend on — <BND> translocations, <INV> inversions, and de novo <INS> — and ships bundled pan-cancer and per-tissue (BRCA / skin / lung) models plus Docker-based benchmark pipelines.

See docs/cancer_howto.md for a copy-paste guide with worked examples, output reference, benchmarking, and model training. Design rationale and calibration caveats live in docs/cancer_simulator.md.

How eidolon compares to NEAT

eidolon is a Rust port of NEAT that tracks the NEAT feature set while adding a native cancer workflow, a low and flat memory footprint, and reproducible output. The table below compares the original Python 2 NEAT (the 2.x “genReads” line), the current Python 3 NEAT 4.x, and eidolon.

Speed and memory are single-thread medians (n=3, 10× coverage) on one exclusive Delta node against NEAT 4.6.1, measured on v4.0.0 (#840); see docs/access_report_draft.md §3.2. The whole human genome runs single-threaded in 80 minutes at 12.3 GB. Through an identical GATK pipeline on chr22 (§3.3), SNP recall is equivalent (0.990 vs 0.989), indel recall is higher (0.978 vs 0.969), and eidolon’s precision is higher for both.

On NEAT 4 and structural variants. NEAT 4.6.1 ships Inversion, Duplication, Translocation, Transposition and CopyNumberVariant classes under neat/variants/, but none is reachable: none is exported from that package’s __init__.py, none is referenced anywhere outside its own file, and VariantTypes.types lists only Insertion, Deletion, SingleNucleotideVariant and UnknownVariant. The mutation model can emit only SNVs, insertions and deletions. The input-VCF path classifies records by REF/ALT length arithmetic, so a symbolic <INV>/<DUP> or a breakend ALT becomes an UnknownVariant — whose get_alt() and get_ref_len() both raise, since its __init__ sets no alt and its metadata is nested one level deeper than those accessors expect. NEAT 4 does not generate structural variants and does not render them from an input VCF either. NEAT 2.1 has no structural-variant code at all. Verified against tags 2.1 and 4.6.1 on 12 Aug 2026.

Citations. NEAT: Stephens et al. (2016), PLOS ONE 11(11):e0167047, doi:10.1371/journal.pone.0167047; and Allen et al. (2026), Journal of Open Source Software 11(121):9056, doi:10.21105/joss.09056. eidolon: doi:10.5281/zenodo.20100558.

Where eidolon fits. NEAT 4.x is a capable, actively developed simulator, and the two tools share most of their core feature set. eidolon is the right choice when you want:

• Cancer simulation — a native tumor/normal workflow (gen-cancer-reads) with configurable purity and an origin-tagged truth VCF, plus CNVs and per-tissue cancer models. This is eidolon-only.
• Intra-tumor heterogeneity — a tumor as a mixture of subclones at distinct cancer-cell fractions, either specified inline, fitted from PyClone-VI / DPClust output, or replayed from a real tumor’s observed VAF. Validated on SEQC2 HCC1395 at r = 0.99 against the intended per-site VAF. Also eidolon-only — the axis the NEAT lineage cannot represent at all.
• Low, flat memory — streaming FASTQ writes keep peak RSS small and roughly constant across genome size and thread count, which matters on shared HPC allocations.
• Reproducibility — the same seed yields byte-identical FASTQ regardless of the thread count.
• Easy deployment — a single self-contained binary with no Python environment to manage. A Bioconda package is in review (bioconda-recipes#67420).

Installing eidolon

Two routes today: a release binary, which needs no Rust toolchain, or building from source. A Bioconda package is pending (below).

A release binary

Each release publishes binaries on the releases page. Pick the one matching your system:

The binary is self-contained — it needs nothing beyond a standard C library, so there is no environment to set up. Download it, make it executable, and run it. Renaming it to eidolon is the usual convenience:

mv eidolon-x86_64-unknown-linux-gnu eidolon
chmod +x eidolon
./eidolon --version

On Windows, rename it to eidolon.exe.

Before v3.5.1 the rename is required. Those releases stop with unrecognized subcommand 'eidolon-x86_64-unknown-linux-gnu' when run under the download name, because the binary reads the name it was invoked by to choose its command (#794). From v3.5.1 any name works. A binary named after a subcommand, such as gen-reads, still runs that subcommand directly.

Let us know if one of these does not work on your system.

Bioconda (not yet available)

A Bioconda package is in review at bioconda-recipes#67420. Until it is merged, conda install -c bioconda eidolon finds nothing. The package Bioconda does carry, rneat, is this project under its earlier name, frozen at 1.20.1; use a release binary instead.

Building from source

You will need the Rust toolchain including cargo; see the cargo installation guide. You will also need git.

git clone git@github.com:ncsa/eidolon.git
cd eidolon
cargo build --release

Compilation may stop on a missing system dependency — cmake is the usual one, since some of the crates eidolon depends on need it. On Debian/Ubuntu that is sudo apt install cmake; on RHEL/Rocky, sudo dnf install cmake. There may be others depending on your system. Drop a comment on an issue if you need specific help.

Build from source when you are testing a change: a released binary is, by construction, not the code you are working on. See Running on HPC for building on a cluster.

$ eidolon --help
Usage: eidolon [OPTIONS] [SUB-COMMAND]

SUB-COMMANDS:
  gen-reads              Generates reads for an input dataset
  filter-reads           Filters the output of gen-reads
  gen-mut-model          Generates a mutation model from real VCF data
  gen-seq-error-model    Generates a sequencing error model from real FASTQ data
  gen-frag-length-model  Generates a fragment length model from a BAM or SAM file
  gen-gc-bias-model      Generates a GC bias model from a reference FASTA and aligned BAM
  gen-bam-models         Builds multiple models (frag-length, GC bias) from one BAM in a single pass
  compare-vcfs           Compares a NEAT-simulated golden VCF against a downstream variant-caller VCF
  gen-cancer-reads       Simulates a tumor/normal mixture (two gen-reads passes + merge)
  validate               Checks an emitted FASTQ/VCF/BAM against what its consumers accept
  compare-af             Per-allele AF correlation between a truth VCF and a simulated one
  help                   Print this message or the help of the given subcommand(s)

Options:
      --log-level [<log_level>]  Verbosity of the written .neat.log. The on-screen log is always
                                 info. [default: info] [possible values: trace, debug, info,
                                 warn, error, off]
      --log-dest <log_dest>      Sets the log destination (full path with full filename) for the written log
  -h, --help                     Print help

validate and compare-af take file paths and flags directly rather than a configuration YAML — both are diagnostic/measurement runs invoked ad hoc and from harness scripts, where writing a YAML per check is friction.

To check options for a subcommand:

$ eidolon gen-reads --help
Generates reads for an input dataset

Usage: eidolon gen-reads [OPTIONS]

Options:
  -c, --configuration-yaml <configuration_yaml>  Path to configuration file.
  -h, --help                                     Print help

To run filter reads, check the help menu.

$ eidolon filter-reads --help
Filters the output of gen-reads

Usage: eidolon filter-reads --configuration-yaml <configuration_yaml>

Options:
  -c, --configuration-yaml <configuration_yaml>  Path to configuration file.
  -h, --help                                     Print help

To run gen-mut-model:

$ eidolon gen-mut-model --help
Generates a mutation model from input VCF data

Usage: eidolon gen-mut-model [OPTIONS]

Options:
  -c, --configuration-yaml <configuration_yaml>  Path to configuration file.
  -h, --help                                     Print help

Use the help menu to see the available options and leave an issue if you find something bad happening.

To compile and run eidolon yourself, besides the rust toolchain, you will need git installed for your operating system. You will then need to git clone and cd into the repo directory. From your home directory in Linux the process might look something like:

~/$ git clone git@github.com:ncsa/eidolon.git
~/$ cd eidolon
~/eidolon/$

For Windows and Mac users, try the binary packaged with the latest version of eidolon.

Once in the repo, you can build the program either in debug (default) or release mode. The main difference is how much info it gives you if there is an error. Release mode also has some optimizations to run it faster.

~/eidolon/$ cargo build --release

If you prefer to run the package directly without using the binary, you can also use

~/eidolon/$ cargo run -- gen-reads -c my_config.yml

Rust will download any required packages. Compiling Rust code is the slowest part of the process. The final binary will be built and the program run immediately after in the second case. To run the program manually, from the repo main dir, run

~/eidolon/$ ./target/release/eidolon -h

eidolon uses a configuration file to read values it needs for the run. A command line execution might look like this:

~/eidolon/$ ./target/release/eidolon -c /path/to/filled/in/config.yml

Every key in a config file must be one the subcommand reads. An unknown key stops the run before it starts, with an error naming the key, the nearest accepted key when one is close, and the full list of accepted keys. That covers a typo (tumor_mutation_model: for tumor_model:), a key belonging to a different subcommand, and keys from a NEAT config, for which gen-reads names the eidolon equivalent where one exists. An unread key used to be ignored, leaving its setting at the default without warning. The templates in template_config/ list every key each subcommand accepts. If you record the output in the logs of Seed string to regenerate these exact results: XXXXXXX, you should be able to use that string as input with rng_seed and reproduce your results.

Fastq Output

The fastq output will have a key name that identifies the block where the read was drawn from, for quick comparisons in alignments. The output BAM file will contain the original sequence and cigar string. The name will have the format EIDOLON_generated_<contig short name>_<fragment_start>_<fragment_end>_<uniq>/1 (or /2 for the second read in a pair), where start and end are zero-padded to 10 digits and <uniq> is a 16-digit hex per-fragment tag that keeps same-position fragments from colliding (see #210).

@EIDOLON_generated_Chromosome_0000000000_0000000353_000000000000002a/1
CTTTTCATTCTGACTGCAACGGGCAATATGTCTCTGTGTGGATTAAAAAAAGAGTGTCTGATAGCAGCTTCTTAACTGGTTACCTGCCGTGAGTAAATTAAAATTTTATTGACTTAGGTCACTAAATACTTTAACCAATATAGGCATAGC
+
>AC7<GDEGGGGEFGA<GFCGG;GGGGGF>GEGGEGGGFGFEFGCEGGGGGGCG:AEFGFFGG>FG;GDGA9$GGAGF=GFG=EFFCGGGFGGGGGC$BFGEFFAGGG9F7E@>?GFGGGG>EBFGFDG)DGC6DEDFA2EG:EGG%FFB

The above example comes from the chromosome in the reference (which was simply named “Chromosome”) between 0 and 353, which is the size of a fragment in the simulated DNA. If there is a pair with this read, it will have the same coordinates, though it started at index 352 instead.

BAM Output

eidolon can write a golden BAM file alongside the FASTQ output. The BAM contains the same reads as the FASTQ — same sequences, same quality scores, same variants and sequencing errors applied — with alignment information included. To enable it, set produce_bam: true in your gen-reads config. The output path is derived automatically from output_filename (e.g. output_filename: my_run → my_run.bam). Please note that the BAM has a higher overhead than the fastq, and may take longer to produce.

produce_bam: true

The CIGAR strings in the BAM reflect the full ground-truth alignment:

• Genomic variants (SNPs, insertions, deletions) are encoded as M, I, and D ops.
• Sequencing error indels are also encoded: deletion errors add D ops for skipped reference bases; insertion errors add I ops for inserted bases.

Every mapped record carries an NM tag, its edit distance to the reference: mismatched aligned bases plus inserted plus deleted bases. It is counted the way samtools calmd counts it: soft-clipped adapter bases are not edits, and an N on either side counts as a mismatch. samtools stats and MultiQC take their mismatch count and error rate from NM, so they report the simulated error rate. Unmapped records carry no NM, and no MD tag is written; samtools calmd can add one.

Reads spanning an SV junction (BND, INV, and the junctions of DEL and DUP) are written as one alignment at the position of their first piece, although part of each read comes from elsewhere. Their NM is correct for that alignment, so it is high, and when SVs are simulated these reads raise the error rate samtools stats reports above the sequencing error rate alone.

The BAM is written in coordinate-sorted order. No post-processing sort is required before indexing:

samtools index my_run.bam

Note: heterozygous variants are applied probabilistically, identical to the FASTQ. A read drawn to the reference allele will not show the variant in either the FASTQ or the BAM.

Targeted Generation with a BED File

For exome, targeted-panel, or any run where you only want reads over specific regions, you can supply a BED file at generation time:

target_bed: /path/to/targets.bed

When target_bed is set, gen-reads skips contigs absent from the BED entirely — no blocks are read, no variants are placed, and no reads are generated for those contigs. Within covered contigs, reads and variants are generated only over the intersecting non-N regions. This is the recommended approach for targeted runs on large genomes; it is far more efficient than generating genome-wide reads and post-filtering with filter-reads.

The BED contig names must match the short names derived from the reference FASTA (the text after > up to the first whitespace character). Both .bed and .bed.gz inputs are accepted.

coverage on a targeted run is a fragment budget, not on-target depth

On a targeted run, expect on-target depth to come in below the coverage you asked for, and scale your request up to compensate. This is a real property of targeted sequencing, not a defect, but the arithmetic is not obvious and coverage does not currently correct for it.

A fragment is owned by the target its start falls in, but its far end runs into the flanking sequence beyond the target — as it must: a 250 bp fragment cannot fit inside a 178 bp exon. Those spilled bases are real reads (they are exactly the off-target reads a real capture kit produces), but they are not on-target depth. So:

delivered on-target depth  ≈  coverage × on-target base fraction

and that fraction is governed by target width against fragment length. Measured on a 436-target BED (mean width 178 bp) over real GRCh38 chr22 sequence, read_len: 100, fragment_mean: 250, coverage: 60:

To get ~60× on target with those parameters, ask for roughly 60 / 0.405 ≈ 148×. The fraction rises quickly as targets get wider relative to fragments, so a panel with 1 kb targets loses far less than a 178 bp exome, and a whole-genome run (no BED) loses nothing.

Two related notes:

• Do not read a 100%-on-target simulation as “better”. eidolon before v3.1.x confined every fragment inside its target, which yields 100% on-target and is physically impossible for targets narrower than the fragment length. It also under-delivered the fragment budget itself (60% of requested, in the run above) and collapsed the realized fragment-length distribution near target edges.
• Making coverage mean on-target depth directly — i.e. having eidolon scale the budget for you — is tracked in #578. Until then, scale it yourself.

Custom Mutation Rate Regions with a BED File

In addition to targeting certain regions, you can use a BED file with a custom field entered in any column after the third.

mutation_regions: /path/to/mutation_regions.bed

The text pattern is “mut_rate=0.0001” followed by a delimiter (end of line, space, semicolon, comma, bar) where the number after the equal sign is any float, and will be the mutation rate for the region defined in columns 1-3.

chr1    39930   39957   Bed_info    mut_rate=0.002
chr2    0   111199  Other_bed_info  mut_rate=0.02

Any region outside of the regions defined in the mutation regions BED will be assigned the default mutation rate, set by the mutation model, which can be overridden a custom default defined in the config file. For example, if you only want to have mutations in exomes, but you want reads from the full genome, you could set the mutation_rate to 0.0 and use a bed file defining all exomes, with a column appende mut_rate=0.0011 or whatever mutation rate you desire. This BED works in conjunction with the target_bed, and any region outside the target bed will be excluded from reads anyway, and thus will have no variants.

Input Variants VCF

You can supply a VCF of variants to force into the simulation:

input_vcf: /path/to/variants.vcf.gz

eidolon will place every variant from the VCF into the corresponding position in the simulated reads and the output VCF. Random variants are still generated at mutation_rate for all positions not covered by the input VCF; set mutation_rate: 0.0 to disable random variants entirely and output only the provided set of variants.

Requirements

• The VCF must be single-sample (one sample column).
• Every record must include GT in the FORMAT field. eidolon uses the genotype to determine whether to apply the variant to all reads covering the position (homozygous, e.g. 1/1) or only a probabilistic subset (heterozygous, e.g. 0/1). Records without GT are rejected.
• Contig names must match the short names derived from the reference FASTA (text after > up to the first whitespace character). Variants on unrecognised contigs are skipped with a warning.
• Both .vcf and .vcf.gz files are accepted.

Supported variant types

Lowercase bases

Bases are case-insensitive. A soft-masked (lowercase) reference is read as plain A/C/G/T, so it simulates exactly as its uppercase copy does. REF and ALT bases from the input VCF are written to the output VCF in uppercase, breakend ALTs included: t]chr2:300] is written T]chr2:300]. The mate’s contig name keeps its case, since contig names are case-sensitive, and <...> IDs such as <INS:ME:Alu> are written as given. VCF 4.2 treats lowercase bases as uppercase, so this changes no variant, only how it is spelled.

Symbolic / structural variants

Symbolic ALTs (VCF 4.2 §1.4) are accepted and round-tripped to the output VCF verbatim apart from base case (see “Lowercase bases” below), with INFO/END, INFO/SVLEN, and INFO/CN preserved. As of v1.10, eidolon can also generate symbolic SVs de novo from a learned model — opt in by setting sv_rate_scale: 1.0 (or higher) in your gen-reads YAML; see “De novo SV generation” below.

DEL anchor convention: for <DEL>, POS is the unaffected base immediately before the deletion (per VCF 4.2), so the modulated span is [POS+1, END] in 1-based coords. For <DUP> / <CNV> / <INV>, POS is the first base of the affected region, so the span is [POS, END].

When an SV zeroes out coverage (hom <DEL> or INFO/CN=0), the mutation rate over the same span is also zeroed so de-novo SNPs don’t pollute the output VCF with variants that never appear in reads.

Current caveats

• Multi-allelic records: only the first ALT allele is used; additional alleles are silently ignored. Split multi-allelic records with bcftools norm -m - before passing to eidolon. (This is true for both literal and symbolic ALTs — <DEL>,<DUP> on the same line is treated as the first ALT only.)
• REF allele verification: eidolon does not check that the REF field matches the reference sequence at that position. Mismatches will produce biologically incorrect output without any warning.
• Literal complex variants: records whose REF and ALT are both multi-base strings (and the ALT is literal bases, not a <TAG>) are skipped with a logged warning and do not appear in the output.
• <INV> interior is forward-strand: inversion junction reads are emitted (the chimeric signal callers use to detect the inversion), but reads from the interior of an inverted span are still transcribed from the forward-strand reference. The <INV> record round-trips to the output VCF.
• Breakends: <BND> junction reads are generated (chimeric reads spanning the mate locus; validated against Manta — see the ACCESS report), and the record round-trips to the output VCF.

De novo SV generation

eidolon can sample SVs directly from a learned SvModel rather than relying on user-supplied records: <DEL>, <DUP>, <CNV>, <INV> and literal insertions are drawn per contig, and <BND> translocations are drawn separately as mate pairs across two contigs. <BND> therefore needs a reference with at least two contigs — on a single-contig reference the correct yield is zero rather than a same-contig junction, so a chr22-only run reports truth BND: 0 by construction. Generation is off by default (opt-in via sv_rate_scale) so v1.9 pipelines remain unchanged.

To enable, add a single line to your gen-reads YAML:

sv_rate_scale: 1.0

• 0.0 (the default) disables de novo SV generation entirely — only input_vcf records flow through.
• 1.0 reproduces the rate from the trained model.
• Larger values scale the rate proportionally for stress testing.

When enabled, eidolon consults MutationModel.sv_model:

1. If you trained your own model with eidolon gen-mut-model -c <yaml> against an SV-rich VCF (e.g. a gnomAD-SV slice), the model file includes a fitted SV component covering type / length / copy-number / homozygous frequencies. Pass it via mutation_model: /path/to/your.json.gz in the gen-reads YAML.
2. If you don’t supply a model, gen-reads loads the bundled default. The default carries approximate gnomAD-SV v2.1 parameters — useful for kicking the tires, but not a substitute for retraining on data that matches your downstream use case. See eidolon-core/src/models/sv_model_defaults.rs for the parameter sources.

Per-contig sampling: Poisson count from per_base_rate × contig_len × sv_rate_scale → weighted type pick → log-normal length (rejected outside [50bp, contig_len / 4]) → uniform anchor with overlap and N-gap rejection → INFO/CN draw for <CNV> → Bernoulli for genotype.

De novo records merge with input_vcf SVs and flow through the same depth-modulation path, so simulated coverage and the golden VCF round-trip behave identically for both sources. compare-vcfs skips both into the skipped.symbolic bucket — symbolic ALTs aren’t byte-comparable.

Caveats:

• The bundled default is a literature-derived approximation, not a refit on the actual gnomAD VCFs. Don’t use it for distribution-faithful benchmarking — retrain.
• <INV> and breakends are not generated (the read content for either isn’t modeled yet).
• A trained model’s sv_model field is None if the training VCF lacked sufficient SV observations (< 2 per type after filtering). Loading that model with sv_rate_scale > 0 is harmless — generation just no-ops.

FASTQ Shuffling

eidolon writes reads in contig order. To shuffle the output, use seqkit shuffle as a post-processing step:

# single-ended
seqkit shuffle sample.fastq.gz -o sample_shuffled.fastq.gz

# paired-ended (keeps mates in sync)
seqkit shuffle -2 sample_R1.fastq.gz sample_R2.fastq.gz \
    -o sample_R1_shuffled.fastq.gz -o sample_R2_shuffled.fastq.gz

seqkit uses reservoir sampling and streams from disk, keeping memory use bounded regardless of file size (seqkit is an open source toolkit for FASTQ files: https://github.com/shenwei356/seqkit).

3′ Adapter Readthrough

Real Illumina libraries whose insert is shorter than the read length “read through” the fragment into the 3′ sequencing adapter — the read’s tail is adapter sequence, not genome. eidolon can simulate this so short-insert data looks realistic and exercises downstream adapter-trimming QC. Disabled by default — output is byte-identical to prior versions when the adapters block is omitted.

adapters:
  enabled: true
  preset: truseq        # truseq | nextera | custom
  # for preset: custom, give explicit 5′→3′ uppercase-ACGT sequences:
  # r1: AGATCGGAAGAGC...
  # r2: AGATCGGAAGAGC...

When enabled, any fragment whose insert is shorter than read_len is kept (instead of being resampled away) and the read is padded to read_len at its 3′ end with adapter sequence — R1 gets r1, R2 gets r2 (appended after R2 is reverse-complemented, matching real read orientation). Adapter bases carry the same quality/error models as genomic bases, are marked soft-clipped (S) in the golden BAM, and never introduce variants. How much readthrough you get is controlled by how many inserts fall below read_len — i.e. by fragment_mean / fragment_st_dev relative to read_len. With a long-insert setting (fragment_mean ≫ read_len) essentially no reads reach the adapter and enabling this changes nothing.

Presets:

• truseq — Illumina TruSeq (AGATCGGAAGAGC…), the most common.
• nextera — Nextera / transposase (CTGTCTCTTATACACATCT).
• custom — supply your own r1 / r2.

Not compatible with long_reads: true (adapter readthrough is a short-read phenomenon).

When to enable it

• Validating adapter-trimming pipelines — generate reads with known adapter content to confirm fastp / cutadapt / Trim Galore settings actually detect and remove it.
• Short-insert libraries — small-RNA / cfDNA / degraded-DNA style data where inserts sit near or below the read length and readthrough is expected; the simulated reads then match what the sequencer really produces.
• Aligner soft-clip benchmarking — untrimmed adapter tails should be soft-clipped by the aligner; use this to check that behaviour end-to-end.
• Realistic FASTQ→VCF tests — exercise the whole trim → align → call path on data that carries adapters instead of idealized adapter-free reads.

To keep short inserts as genomic reads without modeling adapters (e.g. to study the coverage behaviour of a short-insert library on its own), set keep_short_fragments: true instead of enabling adapters — short fragments are then emitted as insert-length reads with no adapter padding.

Validated behaviour (paired germline benchmark: chr22, 30×, TruSeq, 3 reps/arm, BWA-MEM2 → GATK HaplotypeCaller vs the golden VCF): adapter readthrough is detectable (fastp adapter content and aligner soft-clip fraction rise only when it is enabled) and fully callable — SNP and indel recall/precision/F1 are unchanged within replication noise whether the adapters are fastp-trimmed or left for BWA-MEM2 to soft-clip. An adapter-free short-insert control showed the only fidelity difference versus a long-insert baseline is the reduced effective coverage of short inserts (mate overlap), not the adapters.

Parallel Processing

eidolon gen-reads processes contigs in parallel by default using rayon’s work-stealing thread pool. Each contig is an independent unit of work — variant generation, fragment sampling, and FASTQ/BAM writing all happen concurrently across contigs — so references with many contigs scale well across cores.

Output is byte-identical regardless of num_threads (the same seed always produces the same reads in the same order), so you can change the thread count freely without affecting results.

A note on scaling: read generation is largely memory-bandwidth bound, so on references dominated by one or a few large chromosomes the wall-time gain from extra cores is limited — the bottleneck is memory throughput, not CPU. An experimental sub-contig chunk size knob (below) can split a single chromosome across cores, but it is disabled by default because in our benchmarks it did not improve wall time and occasionally regressed it.

Thread count

By default eidolon uses all available logical cores. You can cap the thread count with the num_threads config key:

# use 4 threads instead of all available cores
num_threads: 4

or disable parallelism entirely (useful for debugging or reproducibility testing)

# disable parallelism entirely (useful for debugging or reproducibility testing)
num_threads: 1

Omit num_threads (or set it to .) to restore the default all-cores behaviour.

Thread count and hardware — fewer threads can be faster

Read generation moves a lot of data relative to the arithmetic it does, so it is largely memory-bandwidth bound. On a typical desktop or laptop (which has only a couple of memory channels), a single eidolon thread can already saturate much of the available memory bandwidth — so adding cores yields little speedup and, past a point, can even run slower as threads contend for the memory bus. In our desktop benchmarks, large references ran about as fast (sometimes faster) on 1–4 threads as on all 8.

Practical guidance:

• Desktop / laptop: try num_threads: 1 (or a small number) and compare — it is often as fast or faster than all-cores for a single large run, and leaves cores free for other work. If you are simulating many samples, running several single-threaded eidolon jobs in parallel typically beats one many-threaded job.
• HPC nodes with many memory channels (and multiple sockets) have far more aggregate bandwidth, so higher num_threads scales better there.
• Output is byte-identical regardless of num_threads, so it is safe to tune the thread count purely for speed on your hardware.

Chunk size (experimental, opt-in)

By default each contig is one unit of work. Setting chunk_size splits contigs into sub-contig chunks so a single large chromosome can be worked by several cores at once. It is off by default: read generation is memory-bandwidth bound, so chunking did not improve wall time in our benchmarks (and added a little overhead). It is kept for CPU-bound or very-many-core scenarios where it may help. The size is in base pairs and independent of the thread count, so output stays byte-identical regardless of num_threads.

# default — disabled (one chunk per whole contig); omit or set to . or 0
chunk_size: 0

# opt in: fixed chunk size in base pairs
chunk_size: 5000000

When might it help? Only when read generation is CPU-bound rather than memory-bandwidth bound — the opposite of what we measured on a typical desktop. Consider trying it when both are true:

1. Your reference is dominated by one or a few large contigs, so the default per-contig parallelism leaves most cores idle (e.g. a single-chromosome assembly, or a genome where one chromosome dwarfs the rest).
2. You have memory bandwidth to spare relative to cores — a multi-socket or many-memory-channel HPC node rather than a commodity desktop — and/or you run compute-heavier settings (GC-bias-weighted coverage, long-read mode, very high coverage) where per-read CPU work dominates memory traffic.

It is not worth enabling on a typical workstation, or on references with many contigs (those already parallelize per contig). Always benchmark on your own hardware: start with chunk_size ≈ longest_contig_bp / (4 × cores) (a few chunks per core), compare wall time against the default, and keep it only if it is actually faster.

BAM output is fully parallel

eidolon uses a per-contig temp-file strategy: each contig worker writes its alignment records to a private temporary BAM body file, then a single concatenation pass assembles them in reference order into the final coordinate-sorted BAM.

Reproducibility

Each contig’s random number generator is derived deterministically from the parent seed and the contig’s position in the reference, so output is identical across runs with the same seed even when the number of threads changes.

Reading Bed Data

eidolon can now read bed data and filter the reads based on regions specified. This comes with some caveats. First is that rust can’t handle any header lines or non-header rows in the bed file, because of potential range of possibilities. Second, the names must match what is in the fasta. Third, eidolon can only filter, and does not use the rest of the bed information. It treats each record as a region of interest only.

Bed data will have some challenges. For example, if the contig names in the bed file don’t match the assumed contig name from the fasta file, how will rust be able to know which is what? To make things work, we are going to assume, for now, that the bed file contig names match the names as derived by eidolon. To determine this:

First, in a linux terminal, use your input fasta to find the names of the contigs like this:

$ grep "^>" input.fasta
>NC_001133.9 Saccharomyces cerevisiae S288C chromosome I, complete sequence
>NC_001134.8 Saccharomyces cerevisiae S288C chromosome II, complete sequence
>NC_001135.5 Saccharomyces cerevisiae S288C chromosome III, complete sequence
>etc

This will give you the full fasta name. eidolon determines the short name of the input fasta by skipping the initial ‘>’ character, then taking the string up to the first whitespace delimiter. In the above example, the short names are “NC_001133.9”, “NC_001134.8”, etc.

Next, we check the bed. This example awk command will look for unique values in the first column of your bed file and print out what it finds.

$ awk '!seen[$1]++ {print $1}' file.bed
1
2
3
etc

In this case, the bed file uses a simple numbering scheme to number the chromosomes. We can see a mapping with the roman numeral chromosomes in the name, but it’s tricky to handle consistently. So, the following command will allow you, the user, to map these chromosomes and thus instruct eidolon. Yours will vary based on the inputs.

awk -F'\t' '{$1=($1=="1"?"NC_001133.9":$1); $1=($1=="2"?"NC_001134.8":$1); print}' OFS='\t' file.bed > file_renamed.bed

You will have to tailor this command to your dataset, but this is one way to map the names from the bed to the fasta in a way eidolon will understand.

Cancer simulation how-to

The following guide to eidolon gen-cancer-reads should give you some ideas on how to use eidolon to simula cancer reads and test your cancer pipelines. This is based on previous work done by NCSA and ICGC/ARGO and needs real testing to validate it works as expected. Design rationale and the mutation-rate / SV calibration details live in cancer_simulator.md.

Pipeline

gen-cancer-reads runs two gen-reads passes over one reference and merges them:

  normal pass   cov = (1 − purity)·C    germline only
  tumor  pass   cov = purity·C          germline (shared) + somatic (de novo)
        │
        merge: tag reads N_/T_, concatenate; combine the two golden VCFs with origin tags
        ▼
  <prefix>_merged_r1.fastq.gz    → aligner + somatic caller
  <prefix>_merged_truth.vcf.gz   → INFO/EIDOLON_ORIGIN ∈ {germline, somatic, shared}

The tumor pass consumes the normal pass’s golden VCF as its germline (input_vcf), so both populations carry the same germline. No bcftools/awk runtime dependency — the merge is native.

Quick start

Smoke test on the bundled H1N1 reference (seconds; proves the build works):

cat > cancer.yml <<'YAML'
reference: eidolon/test_data/references/H1N1.fa
output_dir: ./cancer_out
output_prefix: smoketest
total_coverage: 30
purity: 0.6
read_len: 70
paired_ended: true
fragment_mean: 250
fragment_st_dev: 30
rng_seed: demo
overwrite_output: true
YAML

eidolon gen-cancer-reads -c cancer.yml

template_config/gen_cancer_reads_template.yml is the fully-commented config.

Output files

The N_/T_ read-name tags prevent same-coordinate QNAME collisions between the two passes (which MarkDuplicates would otherwise drop).

Train your own model

gen-cancer-reads is model-driven: each pass takes a .json.gz mutation model (tumor_model: / normal_model:). The bundled models are starting points — for real work, train from your own somatic calls.

eidolon gen-mut-model fits a model from a reference plus a single-sample VCF:

cat > my_tumor.yml <<'YAML'
reference: /path/to/GRCh38.fa
vcf_file: /path/to/your_somatic_calls.vcf.gz
output_file: my_tumor_model.json.gz
overwrite_output: true
YAML
eidolon gen-mut-model -c my_tumor.yml

Input VCF expectations:

• single-sample, GT in FORMAT; contig names match the reference;
• SNPs and indels are fit by REF/ALT length class — multi-base REF and ALT (complex) records are skipped;
• symbolic SV records (<DEL> / <DUP> / <CNV> / <INV> / <BND> with SVTYPE / END / SVLEN) are fit into an sv_model when present, which is what sv_rate_scale draws from at simulation time;
• bed_file: restricts the fit to regions. See template_config/gen_mut_model_template.yml.

The fitted mutation_rate is variant_count / reference_length — corpus-aggregated if the VCF pools many tumors, so treat it as a spectrum descriptor, not a per-tumor rate. Set per-tumor somatic burden at simulation time with tumor_mutation_rate (see config knobs). Then:

tumor_model: my_tumor_model.json.gz
# normal_model: my_germline_model.json.gz   # optional; default = built-in germline model

Public-corpus adapters

If you don’t have your own calls, these convert public corpora into a trainable VCF and chain into gen-mut-model (--train --reference):

Pre-bundled, ready to use as tumor_model: without any download: tools/cosmic_per_tissue_{BRCA,skin,lung}.json.gz (per-tissue SNP/indel + SV) and tools/cosmic_v104_pancancer_model.json.gz (pan-cancer).

Worked examples

Assume eidolon on PATH and a reference at ~/code/data/GRCh38.fa.

Tumor/normal with your own model

cat > run.yml <<'YAML'
reference: ~/code/data/GRCh38.fa
output_dir: ./out
output_prefix: tumor70
total_coverage: 60
purity: 0.7
read_len: 250
paired_ended: true
fragment_mean: 350
fragment_st_dev: 50
tumor_model: my_tumor_model.json.gz
tumor_mutation_rate: 1e-5
rng_seed: run1
overwrite_output: true
YAML
eidolon gen-cancer-reads -c run.yml

With structural variants

Enable de novo SVs with sv_rate_scale (requires an sv_model in the tumor model — present in your fit if the training VCF carried symbolic SVs, and in the bundled per-tissue models):

tumor_model: tools/cosmic_per_tissue_BRCA.json.gz
sv_rate_scale: 1.0       # 1.0 = the model's nominal rate; higher stress-tests SV callers

Subclonal architecture (somatic VAF spectrum)

By default the somatic burden sits at ~one effective VAF, set by purity alone. Real tumors are mixtures of subclones at distinct cancer-cell fractions (CCF). Add a subclones: list to spread de-novo somatic variants across those fractions — each variant is assigned a subclone (share ∝ weight) and takes its ccf:

purity: 0.8
subclones:
  - {ccf: 1.0, weight: 0.6}   # clonal / truncal — present in every tumor cell
  - {ccf: 0.4, weight: 0.3}   # major subclone
  - {ccf: 0.15, weight: 0.1}  # minor subclone

CCF is a cellular-fraction factor that composes with the variant’s dosage and with purity — it does not replace them:

observed VAF = purity · dosage · CCF

dosage is the alt-copy fraction from the genotype (0.5 for a heterozygous SNV, 1.0 for homozygous; alt_copies/ploidy in general). So a heterozygous somatic variant at CCF f is observed at ~purity·f/2 — the value a subclonal-deconvolution tool (PyClone, SciClone, …) inverts back to f. These are orthogonal axes (purity = normal contamination, dosage = per-copy multiplicity, CCF = subclonal fraction). Germline (shared) variants are unaffected. ccf ∈ (0, 1]; weight defaults to 1.0 (equal share). Omit the block for the dosage-only default (output is byte-identical to pre-subclone runs).

Each somatic record in the golden/merged-truth VCF carries two ground-truth INFO tags:

• EIDOLON_CCF — the intended cellular fraction (subclone CCF).
• EIDOLON_VAF — the intended observed VAF after tumor/normal mixing (purity × dosage × CCF; for a reproductive somatic_vcf replay, the original input VAF). This is the number a somatic caller measures on the merged reads.

Mind the difference from FORMAT/AF: that field is measured per-pass (tumor-only), so it reads dosage × CCF (≈ EIDOLON_VAF ÷ purity) and carries sampling noise. For scoring a caller’s VAF against ground truth, compare to EIDOLON_VAF; use FORMAT/AF only if you want the tumor-cell fraction. Germline/shared records carry neither tag.

# planted observed VAF vs caller VAF: score against EIDOLON_VAF, not FORMAT/AF
bcftools view -H -i 'INFO/EIDOLON_ORIGIN="somatic"' merged_truth.vcf.gz \
  | awk -F'\t' 'match($8,/EIDOLON_VAF=[0-9.]+/){v=substr($8,RSTART,RLENGTH); sub(/^EIDOLON_VAF=/,"",v); print $2, v}'

Building the architecture from real data

Instead of hand-authoring subclones:, point at a subclonal-deconvolution tool’s cluster table with subclones_file: (tab-separated, header required; mutually exclusive with the inline list). eidolon folds the two shapes real tools emit into {ccf, weight}:

subclones_file: pyclone_clusters.tsv

Other tools convert with a one-liner — emit a header plus the two/three columns above. CCF > 1.0 (noisy clonal clusters) is clamped to 1.0 and non-positive-CCF rows are dropped, both warned with a count. Extra columns (cellular_prevalence_std, cluster_assignment_prob, …) are ignored.

Round-trip validation. Ingest a real tumor’s clusters → simulate → the golden VCF’s EIDOLON_CCF carries exactly those planted CCFs → run the deconvolution tool on the simulated reads → confirm it recovers the architecture you fed in. scripts/delta/run_subclonal_vaf_validation.sh runs the read-level half on real data (Delta): it aligns the merged reads and checks the observed VAF at each somatic site tracks EIDOLON_VAF (fidelity gated on bias + MAE-vs-noise-floor). On a soybean scaffold at 200× with three subclones spanning 4–40% VAF, the aligned reads reproduce the planted spectrum unbiased (mean err −0.003) and to the sampling-noise floor (MAE 0.026, Pearson r 0.95) — see report §3.12.

Structural variants use the same CCF composition. De-novo DEL/DUP/INV/INS/CNV records receive dosage × CCF; reciprocal translocation BND mates receive one shared CCF so the two ends of one event remain balanced. The CCF-adjusted fraction drives interval depth, chimeric/junction reads, and breakpoint double-count suppression. This was validated on Delta with a multi-contig yeast smoke (12 SV records spanning CCF 0.2, 0.5, and 1.0; all records carried EIDOLON_CCF/EIDOLON_VAF; 199 chimeric alignments survived alignment with 100% mapping). The full-GRCh38 stress run then produced 140 tagged SVs across CCF 0.2/0.5/1.0, 2,837 chimeric alignments, 100% mapping, and a 71 GB BAM. The no-subclone path remains unchanged.

SV limitations remain explicit (updated 2026-08-25): inserted novel sequence in input BND ALT alleles is not yet carried through junction reads (#498), symbolic <INS> is not a supported read-level path (#500), literal deletions past ~a read length are only partially realized (#590), and the golden BAM does not represent an insertion longer than a read in any CIGAR (#589). Literal insertions longer than read_len - 1 are now realized — #516 is fixed and its interim cap removed. Non-binary SV depth multipliers have a small over-delivery bias (#499). These limitations affect specific SV representations; they do not remove the CCF tags or the CCF scaling of the supported paths.

Reproductive replay (from a real somatic VCF)

The subclonal options above are generative — they invent somatic variants matching a target architecture. To instead replay a specific real tumor’s somatic calls, point at a somatic VCF:

purity: 0.6
tumor_mutation_rate: 0.0      # pure replay — no de-novo somatic on top
somatic_vcf: tumor_somatic.vcf.gz

Each variant is honored at its observed VAF (INFO/AF, else derived from FORMAT/AD), divided by purity so the merged reads reproduce that VAF after tumor/normal mixing (a raw Mutect2/Strelka VCF works directly). Replayed variants are tagged EIDOLON_ORIGIN=somatic / EIDOLON_PROVENANCE=somatic_input in the truth — distinct from germline even though both come from files. Notes:

• Composes with generation: leave tumor_mutation_rate > 0 (and optional subclones) to layer de-novo somatic variants on top of the replayed set.
• A variant with no AF falls back to its genotype dosage (het ≈ 0.5, hom = 1.0).
• Observed VAF above purity is physically impossible for a somatic variant; it’s clamped to a tumor-cell fraction of 1.0 with a warning.

Validated on real cancer data: replaying SEQC2 HCC1395’s high-confidence somatic SNVs (GRCh38 chr1, 200×) reproduces the tumor’s full empirical VAF distribution — 2,948 sites across every decile — at Pearson r = 0.99, unbiased (−0.004), MAE 0.024 (see report §3.12; harness scripts/delta/run_hcc1395_reproductive.sh).

One germline, many tumor scenarios

Fix the germline once and sweep purity/depth by pointing each run at the same germline_vcf: (any normal-pass golden, or a real germline VCF):

for p in 0.3 0.5 0.8; do
  sed "s/^purity:.*/purity: $p/; s/^output_prefix:.*/output_prefix: p$p/" run.yml \
    | sed "/^tumor_mutation_rate:/a germline_vcf: ./out/tumor70_normal.vcf.gz" \
    > run_$p.yml
  eidolon gen-cancer-reads -c run_$p.yml
done

Benchmarking

Docker-based scoring pipelines, reads → scored caller in one command:

# SNV/indel: BWA-MEM → Mutect2 → som.py
tools/cancer_benchmark.sh \
    --reference ~/code/data/GRCh38.fa \
    --normal-fastq ./out/tumor70_normal_r1.fastq.gz \
    --tumor-fastq  ./out/tumor70_merged_r1.fastq.gz \
    --truth-vcf    ./out/tumor70_merged_truth.vcf.gz \
    --output-dir   ./bench

# SV: BWA-MEM → Manta → truvari
tools/cancer_sv_benchmark.sh \
    --reference ~/code/data/GRCh38.fa \
    --normal-fastq ./out/tumor70_normal_r1.fastq.gz \
    --tumor-fastq  ./out/tumor70_merged_r1.fastq.gz \
    --truth-vcf    ./out/tumor70_merged_truth.vcf.gz \
    --output-dir   ./sv_bench

INFO/EIDOLON_ORIGIN lets you filter the truth to somatic-only before scoring (--truth-filter).

Config knobs

Relationship to tools/cancer_simulate.sh

The native subcommand is a drop-in replacement for the original shell orchestrator, verified equivalent by eidolon/tests/cancer_parity.rs (identical merged-FASTQ record multisets, per-pass golden VCFs, and origin classifications). The script is retained for the Docker benchmark wiring; prefer the subcommand for new work.

Cancer simulator design

Status: shipped — the cancer epic (#129) is closed. The two-pass orchestration, the native gen-cancer-reads subcommand, SV/CNV support, per-tissue models, and the PCAWG SV refit all landed across v1.12.0–v1.16.0 (see CHANGELOG); #183–#192 are closed. This document is retained as the design/rationale record — read the “Update” notes inline for current state.

Goal

Produce simulated short-read sequencing data with the statistical and structural properties of a real cancer sample, including:

1. A mixture of tumor + normal DNA at a configurable purity (real tumor biopsies are typically 30–80% tumor cells; the rest is normal-tissue contamination).
2. Germline variation shared between the two populations (the tumor and the normal both came from the same patient).
3. Tumor-only somatic variation in the tumor population — elevated SNV / indel rates, cancer-specific mutation signatures, and a richer SV catalog than the gnomAD-SV-derived germline default.

The output must be benchmarkable — a downstream somatic variant caller (Mutect2, Strelka, GRIDSS, Manta, …) should be able to run on eidolon’s simulated reads and have its calls scored against a ground-truth VCF that distinguishes germline from somatic.

Approach: orchestration first, native subcommand later

The recommended v1 path is a thin orchestration layer in tools/ that drives eidolon gen-reads twice — once per population — and concatenates the resulting FASTQs. No core eidolon changes for the MVP; we lean on the existing primitives (input_vcf, mutation_model, sv_rate_scale, paired-end fragment generation, depth modulation).

Why this approach:

• Validated by precedent. NEAT 2.1 shipped this exact workflow as its genReadsTumorTutorial (~30-line shell script, archived at ~/code/neat2/models/genReadsTumorTutorial.zip). It worked then; it’ll work now.
• Minimum new code. Until the simulation produces credible output end-to-end, every line of new core code is risk. The orchestration path gets us to a working artifact with zero eidolon-core diffs.
• Reversible. If we later decide a native eidolon gen-cancer-reads subcommand is worth the engineering, the orchestration script is the design spec for that subcommand’s behavior — it codifies the data flow.

A native subcommand was not an explicit non-goal — it was a v2+ option pending validation that the orchestration produces samples real somatic callers handle correctly. That validation passed, and the native subcommand shipped (v1.15.0, #239) — see the Update below.

Update (#239): the native subcommand eidolon gen-cancer-reads is now implemented — the orchestration validated cleanly (Mutect2 SNV + Manta SV recall), so the two-pass flow is now a first-class subcommand: eidolon gen-cancer-reads -c <yaml> runs both passes and merges them (tagged FASTQs + an origin-tagged truth VCF), with no bcftools/awk runtime dependency. Config template: template_config/gen_cancer_reads_template.yml; design + decisions: docs/archive/cancer_simulator_native_plan.md (archived — shipped). tools/cancer_simulate.sh is retained; the same-seed parity test (eidolon/tests/cancer_parity.rs) landed in v1.15.1.

Data flow

                  ┌─────────────────┐
reference FASTA ──┤                 │
                  │  Pass 1: normal │── normal_reads.fq.gz    ─┐
                  │                 │── normal_golden.vcf.gz   │
                  │  cov = (1−p)·C  │      ↓                   │
                  └─────────────────┘      │ (germline truth)  │
                                           │                   │
                  ┌─────────────────┐      │                   ├──► concat ──► merged.fq.gz
reference FASTA ──┤                 │      │                   │
        normal_golden.vcf.gz   ─────┤  Pass 2: tumor  │── tumor_reads.fq.gz   ─┘
                                    │                 │── tumor_golden.vcf.gz
                                    │  cov = p·C      │      ↓
                                    │  + tumor model  │      │ (germline + somatic)
                                    └─────────────────┘      │
                                                             ↓
                                        normal_golden + tumor_golden ──► merge ──► truth.vcf.gz
                                                                          (origin tags:
                                                                           germline / somatic / shared)

Definitions:

• p = purity, fraction of merged reads coming from the tumor pass. Real tumor purity 0.3–0.8 typical; biopsies of poor-yield tissues can go lower.
• C = total per-base coverage of the merged output (e.g. 30× for WGS).
• The germline VCF from pass 1 feeds pass 2 as input_vcf:, guaranteeing the tumor cells carry the same germline variants as the normal cells (this is biologically correct — both populations descend from the same zygote).
• Somatic variants are added only in pass 2 via the tumor’s mutation_model: and sv_rate_scale:. They appear in tumor_golden.vcf.gz but not in normal_golden.vcf.gz.

Component breakdown

Each component maps to a sub-ticket of #129:

#184 can develop in parallel with #183 using a hand-rolled or NEAT-2.1-recovered model as a stand-in. #186 blocks on #183 (data) but not on anything else.

Source-data decisions

Chosen for v1

Two open-tier sources are supported, each with its own adapter in tools/. The two are complementary, not competing — pick by access constraints and indel needs.

Both adapters can chain into eidolon gen-mut-model via --train. The output mutation model is structurally identical regardless of which adapter produced the input VCF — “cancer-ness” is in the training data, not the schema.

A pre-trained pan-cancer model fit from the COSMIC adapter against v104 is checked into the repo at tools/cosmic_v104_pancancer_model.json.gz (9.2 KB). It’s a temporary placement so the cancer simulator can be exercised end-to-end without requiring users to download the ~1 GB COSMIC corpus first. Final placement (bundled / hosted with a fetch helper / user-trains-their-own) is tracked at #186.

A note on the fitted mutation_rate. gen-mut-model normalizes variant count by total reference length, so a corpus aggregated across many donors (COSMIC, MC3) will produce a mutation_rate that reflects “fraction of bp with any observed mutation across the catalog,” not per-tumor burden. The v104 COSMIC fit lands at ~5.5e-3 — high if mistaken for a per-tumor rate. Users simulating a specific tumor purity should treat this as a corpus-level constant and shape the per-sample burden via coverage and purity in the orchestration layer, not by rescaling the model’s rate.

Considered and ruled out (for v1)

Mutation model: no new schema needed

NEAT 2.1’s three tumor models (BRCA, CLLE, SKCM) used the same pickled-dict schema as its germline NA12878 model. The “tumor-ness” came purely from the input training data — the model schema captures whatever statistics the input VCF supplies.

eidolon’s existing eidolon_core::models::mutation_model::MutationModel (mutation_rate, snp_trinuc_model, insertion_lengths, deletion_lengths, sv_model, homozygous_frequency, variant_dist) is structurally sufficient to serve as a tumor model. No new CancerModel type is needed for v1. The cancer-ness is a property of the training corpus, the orchestration, and the SV gap closures — not of the model schema.

The one feature NEAT 2.1 had that eidolon does not is HIGH_MUT_REGIONS — a per-region rate scaler (~1.7k entries in their BRCA model, used to model mutation-rate hotspot beds). This is a v2 extension (tracked in #413); the v1 MVP does not need it.

SV gap analysis

The cancer MVP’s credibility depends on eidolon being able to generate the SVs that define cancer karyotypes. v1.10.1 generated <DEL> / <DUP> / <CNV> (covering the coverage-altering axis); v1.12.0 closes the structural axis by adding <INS> / <INV> / <BND> — which is what most clinically-recognizable cancer SVs depend on.

⚠️ Retraction (2026-08-01): #187 was marked shipped; translocations were never generated

The de novo SV sampler cannot place a breakend mate on another contig. Both sv_model.rs:810 and :936 hardcode the mate contig to Some(contig_name.to_string()), and sample_variants (sv_model.rs:529-539) receives one contig’s name, length and sequence — it has no view of any other contig. Every BND eidolon has ever generated from a model is intra-chromosomal, by construction rather than by chance. Job 20719077 emitted 466 junctions and 466 of 466 were same-contig.

Four of the five exemplars in the row above (BCR-ABL, PML-RARA, EWSR1-FLI1, MYC-IGH) are inter-chromosomal fusions. None of them is reachable by the shipped generator.

What did ship, and is correct: chimeric junction-read generation and all four VCF 4.2 breakend orientation forms. The read path even supports off-contig mates — runner.rs:2360 resolves reference.get(&mate_contig) — so an inter-chromosomal BND supplied via input_vcf works today. Only the de novo sampler, which is what produces every cancer benchmark, cannot emit one.

How the mislabel arose. The BND share (23.2%) is the PCAWG per-donor mean of svclass == "TRA" rows (normalize_pcawg_sv_model.py:108-114) — TRA being precisely the inter-chromosomal class. But build_pcawg_sv_vcf.py:127-128 reads only f[0],f[1],f[3],f[4],f[10], and its BND branch (:155-159) emits a bare <BND> anchor and continues — discarding the mate chromosome, and never reading strand1/strand2 (fields 8, 9) at all. So we inherited the count of translocations while dropping the two facts that make one a translocation.

Consequently the shipped model carries no chromosome-pair, strand, or geometry field, and geometry falls back to a uniform 25% per form (sv_model.rs:101-103). About a quarter of emitted BNDs are therefore deletion-like direct joins on one contig — the reads carry a plain deletion while the truth VCF types it SVTYPE=BND.

A further conflation worth stating: BND in a VCF is a notation for an adjacency a caller did not resolve, not an event class. Manta buckets inversion-oriented junctions as BND; Delly types the same junction <INV>; both were observed on identical data in job 20745149. Treating BND as an event type with its own 23.2% share encodes which caller and chemistry produced the source calls, not tumour biology.

Remediation is planned as a measurement pass over PCAWG first (intra/inter split, strand spectrum, partner-distance distribution) with the implementation scope decided by what that finds, rather than assumed now.

Roadmap

Stage 0 — foundation (already shipped, v1.10.1)

• ✅ <DEL> / <DUP> / <CNV> de novo generation
• ✅ Lean VCF reader (enables training on whole-genome corpora)
• ✅ Full-genome gnomAD-SV-derived germline SV defaults

Stage 1 — cancer MVP (orchestration validation)

1. ✅ #189 — aneuploidy verification (resolved). Depth-modulation math (build_coverage_multipliers + coverage_multiplier_for) verified to handle whole-contig <DEL> / <DUP> / <CNV> events correctly at chromosome scale via input_vcf: — see pipeline_e2e::whole_contig_* tests. De-novo aneuploidy via SvModel::sample_variants is bounded by a configurable length cap: the sampler rejects lengths above contig_len * sv_max_length_fraction to keep its overlap-rejection search space tractable. The default fraction is 0.25 (contig_len / 4), which reproduces the historical behavior bit-for-bit and is correct for human chromosomes — whole-chromosome and most whole-arm events still can’t be generated by random draws under the default. #229 made this cap a gen-reads config knob: set sv_max_length_fraction: toward 1.0 to allow large de-novo events (the sampler also scales its retry budget proportionally and warns when it starves). This is most useful on small contigs (bacteria, viruses), where the default 0.25 sits below the trained DEL median and starves DEL/DUP/INV into BND; for chromosome-scale aneuploidy with realistic per-arm frequencies, the supported workflow remains input_vcf: — supply segmentation-call output (GATK-CNV, ASCAT, PURPLE, etc.) and the depth-modulation primitive handles it correctly. A dedicated “aneuploidy event” sampler (with per-arm frequencies rather than a global rate) is a v2 follow-up, intentionally deferred.
2. #183 — survey + commit to a tumor training corpus. Shipped: TCGA MC3 MAF adapter (tools/fetch_tumor_corpus.sh) and COSMIC GenomeScreensMutant VCF adapter (tools/fetch_cosmic_corpus.sh).
3. #184 — orchestration script. Can develop in parallel with #183 using a stand-in tumor model.
4. #186 — train initial tumor models from the chosen corpus.
5. #185 — variant-origin annotation in the merged truth VCF. Shipped: every golden VCF now carries INFO/EIDOLON_PROVENANCE (denovo for sampled variants, input for those carried through from input_vcf:). tools/cancer_simulate.sh resolves these to INFO/EIDOLON_ORIGIN (germline/somatic/shared) via a bcftools isec post-merge — see the script’s “Merge golden VCFs” section.

Exit criterion: running the orchestration script with reasonable defaults produces a merged FASTQ that a current somatic SNV caller (Mutect2 or Strelka) calls into a VCF that scores reasonably well against the origin-tagged truth.

Stage 2 — cancer SV credibility (shipped in v1.12.0)

6. ⚠️ #187 — <BND> junction reads (not translocations; see the retraction above). Chimeric-pair generation via process_chimeric_variants → generate_chimeric_pair → get_bnd_pieces + get_stitched_sequence, correctly handling the four VCF 4.2 breakend orientation forms (t[p[, t]p], [p[t, ]p]t). Paired-BND deduplication via canonical-ID HashSet. All of that is real and verified — but orientation is not translocation, and the de novo sampler emits only same-contig junctions.
7. ✅ #188 — <INV> inversions with strand flipping. Two-junction model (start of inversion + end of inversion), each independently sampling reads. generate_inv_pair + get_inv_pieces handle orientation flipping per junction. Smarter fragment-length sampling + graceful TruncatedRead handling shipped alongside, applied to BND too.
8. ✅ #190 — de novo <INS> with novel sequence. Emitted as LITERAL Insertion records (REF="A", ALT="A<novel-bases>") — same shape Manta / DELLY produce when they resolve the inserted sequence. Routes through the existing literal-insertion machinery. Novel sequence drawn from single-base composition of a ±250 bp window around the anchor.

Exit criterion: a current somatic SV caller (GRIDSS or Manta) calls translocation breakpoints from the simulated data with high enough recall that the simulator could plausibly be used as a benchmark fixture.

Status: ⚠️ NOT MET — previously recorded as “validated” (retracted 2026-08-01). The criterion names translocation breakpoints, and no run could have tested it: the de novo sampler cannot place a mate off-contig, so no translocation was ever present in the data the callers were scored against. The evidence cited for “validated” measured something else — same-contig junction recall — and was accurate about that.

What the cited evidence does support: callers recover same-contig breakend junctions at a measurable rate. Job 20745149 puts Manta at 0.738 over the 480 inversion-oriented records it can represent as breakends (the aggregate 0.380 is floored by our own SVTYPE choice, since direct junctions are correctly called <DEL>/<DUP:TANDEM> and cannot match a BND record). Delly emits no BND records at all — it types the same junctions <INV> — which is itself evidence that BND is a representation convention rather than an event class.

Re-stating this criterion honestly requires either generating real translocations or narrowing the criterion to what is actually produced. That decision is deferred to the PCAWG measurement pass.

Known limitations carried into v1.12.0

• Breakpoint double-counting — resolved in v1.14.1 (#236). The regular per-contig pass used to also generate reads across BND / INV breakpoint positions (from the unbroken reference), leaving those loci with regular + junction reads — ~2× coverage for homozygous variants. Fixed: the regular pass now suppresses the broken-allele fraction of junction-crossing read-pairs (collect_suppressible_junctions / suppress_junction_double_count; DEL/BND/INV/CNV-loss, with DUP/CNV-gain intentionally excluded).
• Default model has no BND / INV. The bundled gnomAD-SV-derived default_sv_model() filters BND / INV out at training time. Resolved in v1.12.1 (#217) via tools/inject_cancer_sv_model.py: the bundled COSMIC tumor model (tools/cosmic_v104_pancancer_model.json.gz) now carries a literature-derived sv_model covering all six SV types. The germline-side default_sv_model() was already populated at heuristic rates. A data refit from real cancer-SV calls shipped in v1.14.0 (#218) — the bundled tumor model’s sv_model is now counted from the PCAWG consensus SV/CNV callsets (+ gnomAD-SV for INS length); see the “Cancer-specific SvModel defaults” note below.
• FORMAT/AD/DP/AF asymmetric across SV types. Symbolic SVs (BND, INV, DEL, DUP, CNV) emit .,.:.:. placeholders since span-based depth doesn’t fit a point-based counter. INS is now literal and gets real counts. (Span-based depth for symbolic SVs is tracked in #411.)
• No mobile-element library for de novo INS. Insertions are random novel sequence, not drawn from L1 / Alu / SVA / HERV consensus. Library-driven path is a #190 v2 extension; viral integration is a separate design problem.

Stage 3 — advanced cancer patterns (defer until stage 2 lands)

9. #191 — chromothripsis / chromoplexy / BFB design pass.
10. #192 — ecDNA / double-minute design pass.

Stage 3 starts with design notes, not implementation. Either of these could prove out-of-scope for eidolon and be better served by a sibling tool that consumes eidolon’s output.

Open questions

These are documented here so future work has a record of what’s been decided vs what’s still up for grabs:

• Subclonal heterogeneity. The v1 orchestration modelled two populations (a single global purity), collapsing every somatic variant to ~one effective VAF. Shipped via #405 (native sampler). A tumor is now a mixture of subclones at distinct cancer-cell fractions (CCF): each somatic variant’s observed alt fraction composes as purity × dosage × CCF. The architecture is authored inline (subclones:), fitted from real data (subclones_file: — PyClone-VI / PCAWG cluster tables), or replayed from a real somatic VCF (somatic_vcf:, honored at its observed VAF). Ground truth is emitted per somatic record as INFO/EIDOLON_CCF (intended CCF) and INFO/EIDOLON_VAF (intended observed VAF). Real-data validated on Delta two ways: generative (soybean scaffold, 200×) reproduces a synthetic 3-cluster spectrum unbiased (mean err −0.003, MAE 0.026); reproductive replay of a real tumor (SEQC2 HCC1395 somatic SNVs, GRCh38 chr1, 200×) reproduces its full empirical VAF distribution (2,948 sites, all deciles) at Pearson r 0.99, unbiased (−0.004), MAE 0.024 — see report §3.12 and scripts/delta/{run_subclonal_vaf_validation,run_hcc1395_reproductive}.sh.

• Mutation signatures (COSMIC SBS). Cancer mutation patterns are dominated by tissue- and exposure-specific signatures (UV → SBS7 in melanoma, tobacco → SBS4 in lung, MMR deficiency → SBS6/15 in colorectal, etc.). eidolon’s trinucleotide-frequency model captures these per-corpus, and context-weighted SNP placement (#372, v1.20.0) now reproduces context-specific signatures — validated on HCC1395 (SBS-96 cosine 0.72 → 0.99) and via signature extraction on chr1–3 (SBS1 + APOBEC recovered; #320 closed). What remains optional is explicit signature authoring (load a named COSMIC SBS profile and mix signatures at chosen proportions), a separate enhancement.

• Subclonal structural variants. Shipped via #537. De-novo SVs now receive the configured CCF, reciprocal translocation BND mates share one CCF, and INFO/EIDOLON_CCF / INFO/EIDOLON_VAF are emitted in the somatic truth. CCF scales supported SV depth and junction evidence as dosage × CCF; the no-subclone path is unchanged. Delta multi-contig validation exercised CCF 0.2/0.5/1.0 and retained 199 chimeric alignments. The remaining SV caveats, updated 2026-08-25: #498 (input BND inserted sequence), #500 (symbolic INS), #590 (literal deletions past ~a read length), and #589 (the golden BAM does not represent an insertion longer than a read in any CIGAR). #516 is fixed — large literal insertions are realized and the interim read-length cap is gone. #499’s over-delivery is largely gone as measured against a control run, but its mechanism is not fully explained (#584).

• Variant allele fraction (VAF) annotation in the truth VCF. Shipped via #176. Golden VCFs now carry FORMAT/GT:AD:DP:AF on every literal SNP/indel record, populated from a per-variant allelic-depth counter incremented by the gen-reads fragment loop at the existing per-read coin-flip site (fastq_tools.rs::generate_read). Symbolic SVs emit . placeholders for AD/DP/AF since their depth semantics are span-based, not point-based — that’s tracked in #411. Realistic FILTER values (LowQual etc.) were left out of v1; everything is still PASS.

• Cancer-specific SvModel defaults. First shipped in v1.12.1 (#217) as a literature-derived sv_model injection (tools/inject_cancer_sv_model.py, now deprecated). Replaced in v1.14.0 (#218) with a data-derived refit: the bundled tumor model’s sv_model is counted from the PCAWG consensus SV/CNV callsets (Li et al. 2020, n=2,748 donors) with gnomAD-SV supplying INS length, via tools/fetch_pcawg_sv_corpus.sh → tools/build_pcawg_sv_vcf.py → gen-mut-model → tools/normalize_pcawg_sv_model.py. The refit corrected the heuristic type mix (notably INV 0.08→0.134, BND 0.35→0.232) to match the published PCAWG ranking (DEL>DUP>BND>INV); the fitted per-base rate 3.46e-8 confirmed the prior 3.8e-8 estimate. INS and focal-CNV rates remain literature estimates (PCAWG MEI calls are controlled-access; focal-CNA segment counts are segmentation, not discrete events). The germline-side default_sv_model() stays gnomAD-derived. See CHANGELOG v1.14.0.

• Per-tissue tumor models (#202). Shipped — both halves now tissue-specific. Three complete tissue models are bundled alongside the pan-cancer one: tools/cosmic_per_tissue_{BRCA,skin,lung}.json.gz. Each pairs a per-tissue COSMIC SNP/indel spectrum with a per-tissue sv_model.

  SV half (first shipped in #237). The sv_model is fitted from that tissue’s PCAWG donors. tools/build_pcawg_sv_vcf.py gained --sample-sheet + --projects (filtering donors by aliquot_id → dcc_project_code via pcawg_sample_sheet.tsv), and tools/build_per_tissue_sv_models.sh drives build→fit→normalize per tissue, emitting tools/cosmic_pancancer_sv_{BRCA,skin,lung}.json.gz (pan-cancer SNP/indel base + per-tissue SV). Tissue → PCAWG project codes (SV-donor counts): BRCA → BRCA (211); skin → SKCM-US,MELA-AU (106, melanoma); lung → LUAD-US,LUSC-US (85). PCAWG’s codes diverge from TCGA’s and SKCM/LUAD alone are underpowered (~37 each), so skin/lung group related projects. The fits show expected signal — BRCA DUP-dominant (Dup 0.317, the tandem-duplicator phenotype) at ~2× SV burden; skin BND-enriched (0.306); lung DEL-dominant (0.349).

  SNP/indel half (#202 completion). Each tissue’s SNP/indel spectrum is trained from its COSMIC GenomeScreensMutant subset. Because the COSMIC VCF is tissue-aggregated, the split uses the sample-level TSV export: tools/fetch_cosmic_per_tissue_corpus.sh reads the Classification file (COSMIC_PHENOTYPE_ID → PRIMARY_SITE), collects every COSV mutation ID seen in that tissue’s phenotypes from the mutant TSV, and filters the GenomeScreensMutant VCF (which carries VCF-anchored alleles, so indels need no re-anchoring) down to those IDs. tools/build_per_tissue_models.sh orchestrates the full chain (fetch → gen-mut-model train → tools/graft_sv_model.py grafts the per-tissue sv_model on). The fits reproduce known tissue biology: skin is the most SNV-dominated (melanoma’s UV-driven C>T; indels only ~1.6% of records), while breast carries the most indels (~5.2%); lung sits between (~3.6%).

  Pick one at run time: gen-cancer-reads with tumor_model: tools/cosmic_per_tissue_BRCA.json.gz.

  Caveats: the bundled mutation_rate is COSMIC-corpus-aggregated (fraction of bp mutated across that tissue’s catalog), so it overstates single-tumor burden — set per-tumor burden via tumor_mutation_rate (defaults to 1e-5); what’s genuinely per-tissue is the spectrum (trinucleotide context, indel ratio, SV mix). INS length is shared from the gnomAD-derived base (germline, tissue-agnostic); INS and focal-CNV rates remain pan-cancer literature values, so a tissue’s CNV fraction shifts mainly because its SV burden (the denominator) differs. The cosmic_pancancer_sv_* files (pan-cancer SNP/indel + per-tissue SV) are retained as the sv_model donor for the graft step.

References

• Umbrella issue: #129 — Cancer modeling
• Cancer MVP sub-tickets: #183, #184, #185, #186
• SV gap tickets: #187, #188, #189, #190, #191, #192
• NEAT 2.1 prior art: ~/code/neat2/models/genReadsTumorTutorial.zip — 30-line simulate.sh reference implementation; ~/code/neat2/utilities/genMutModel.py — original tumor-model trainer that parsed ICGC simple-somatic-mutation TSVs.
• gnomAD-SV-derived germline defaults that the cancer side will complement: eidolon-core/src/models/sv_model_defaults.rs (refit in v1.10.1).

Generating a Mutation Model

eidolon now includes the ability to read in real data and learn the parameters to reproduce it in a simulation run. So far, this is limited to the mutation model, using the eidolon gen-mut-model subcommand.

$ eidolon gen-mut-model -c gen_mut_model_config.yml

The inputs are a single-sample VCF and the reference FASTA the VCF was called against. eidolon computes statistics for indels and SNPs and builds the trinucleotide model of SNP generation, as in the original eidolon. The output is a gzipped JSON model file that can be passed directly to gen-reads via its mutation_model config key.

Copy template_config/gen_mut_model_template.yml to a directory of your choosing and fill in the fields:

# required: path to the reference FASTA used to call the VCF
reference: /path/to/reference.fa

# required: path to a single-sample VCF file
# The VCF must include FORMAT and SAMPLE columns, and GT must be present in FORMAT.
# QUAL=. is accepted and treated as 0.
vcf_file: /path/to/variants.vcf

# required: output path for the generated model; should end in .json.gz
output_file: /path/to/output_model.json.gz

# optional: restrict model learning to regions in this BED file
# set to . (dot) to use the entire reference (default)
bed_file: .

# optional: set to true to overwrite an existing output file (default: false)
overwrite_output: false

VCF requirements

• Single sample only — multi-sample VCFs are not yet supported (tracked in #412)
• Each variant record must have GT in the FORMAT column; eidolon hard-errors if GT is missing
• QUAL=. is accepted and treated as quality score 0
• Insertion and deletion lengths are each fitted from the VCF. If one type is absent (an insertions-only VCF, say), that type’s length distribution is the default one and a warning says so; the type that was observed keeps its fitted lengths. With no indels at all, both use the defaults.
• Each SNP’s and indel’s REF must match the reference at its position. At the first record that disagrees, gen-mut-model stops with an error naming that record and both bases. A mismatch usually means the VCF was called against a different build, or lifted over badly. If the VCF does belong to this reference, drop the mismatching records deliberately and fit on the result:

  bcftools norm -f reference.fa --check-ref x in.vcf.gz -Oz -o checked.vcf.gz
  

• A SNP at a contig edge has no trinucleotide context. It is left out of the model entirely (not counted toward the mutation rate or the SNP/indel split), with a warning giving the count. If every SNP is at an edge, gen-mut-model errors out rather than writing a model with no SNP context data.
• A bed_file must name at least one contig in the reference. If none match, the error lists both sets of names; a chr1 versus 1 mismatch is the usual cause.

Caveats: Only one sample can be read at this point (#412). Currently, high-mutation regions and common variants features from Python NEAT are not yet implemented (#413).

Try it out and let us know if you run into any issues!

We will continue to improve this process in the future. If you have suggestions for parsing names to facilitate this, let us know in the Issues section!

Generating a Sequencing Error Model

eidolon can also learn a sequencing error model from real FASTQ data using the eidolon gen-seq-error-model subcommand. This reads per-base quality scores from the FASTQ to build a Markov quality score model and computes an average base error rate. Optionally, you can supply a BAM file or a custom TSV to set the SNP transition matrix (which base errors are most likely to become which other bases).

$ eidolon gen-seq-error-model -c gen_seq_error_model_config.yml

The output is a gzipped JSON model file that can be passed directly to gen-reads via its seq_error_model config key.

Copy template_config/gen_seq_error_model_template.yml to a directory of your choosing and fill in the fields:

# required: path to input FASTQ file (.fastq or .fastq.gz)
fastq_file: /path/to/reads.fastq.gz

# required: output path for the generated model; should end in .json.gz
output_file: /path/to/seq_error_model.json.gz

# optional: set to true to overwrite an existing output file (default: false)
overwrite_output: false

# optional: number of reads to fit from, sampled uniformly across the file; 0 = every read (default: 0)
max_reads: 0

# optional: quality score ASCII offset; 33 for Illumina 1.8+/Sanger, 64 for older Illumina (default: 33)
qual_offset: 33

# optional: largest read length to model in bp; 0 = no limit (default: 1000)
# One transition matrix is stored per read position (~69 KiB each), so this bounds memory:
# 1000 bp is about 67 MiB. A longer read is an error, not a silent truncation. Raise it for
# genuinely longer reads; note that long-read error models are not supported yet (#319).
max_model_read_length: 1000

# optional: list of Q-score bins to quantize the learned model to (e.g. NovaSeq 6000 uses
# [2, 12, 23, 37]). When set, each observed Q-score is snapped to its nearest bin before
# the model is built, and gen-reads will emit only these values. Omit for a continuous model.
# Bins must be < 94 and cannot include 31 (encodes to '@' under Phred+33).
# binned_quality_bins: [2, 12, 23, 37]

# optional: aligned BAM file to infer the SNP transition matrix from read-vs-reference mismatches.
# The BAM must have MD tags. If your aligner did not add them, run:
#   samtools calmd -b aligned.bam reference.fa > aligned_with_md.bam
# Ignored if transition_matrix_file is also set.
# If the BAM yields no mismatches at all, eidolon errors out rather than quietly using the
# default matrix. To use the default deliberately, omit this key.
bam_file: /path/to/aligned.bam

# optional, requires bam_file: the sample's own variants. Mismatches at the positions its
# records cover are left out of the matrix, since they are variants, not sequencing errors.
known_variants_vcf: /path/to/sample_variants.vcf.gz

# optional: how bam_file fits the matrix. `overlap` (default) counts only disagreements
# between overlapping mates, which read the same molecule, so variants and library damage
# cancel. `mismatch` counts every read-vs-reference mismatch.
bam_method: overlap
bam_min_mapq: 20

# optional: custom 4x4 SNP transition matrix TSV (rows/columns: A C G T).
# A single header line is allowed. Diagonal values are ignored. Each row needs four
# finite, non-negative values with some weight off the diagonal; a malformed row stops
# the run and names its line.
# Takes precedence over bam_file.
transition_matrix_file: /path/to/matrix.tsv

SNP transition matrix priority

1. transition_matrix_file (explicit TSV) — highest priority
2. bam_file (inferred from MD-tagged BAM mismatches)
3. Built-in default matrix (inherited from Python NEAT) — used when neither is provided

Supplying bam_file is a request to fit the matrix from data, so if the BAM yields no read-vs-reference mismatches, eidolon errors out instead of falling back to the default. A model that looks trained and is actually the default is indistinguishable from a trained one downstream. Omitting bam_file is how you ask for the default matrix.

BAM MD tag requirement

The BAM path requires MD tags to identify reference bases at mismatch positions. MD is an optional SAM tag, so check your aligner: BWA-MEM writes it, minimap2 needs --MD, and STAR needs --outSAMattributes MD. If your BAM lacks them, generate them with:

samtools calmd -b aligned.bam reference.fa > aligned_with_md.bam
samtools index aligned_with_md.bam

Without MD tags the inference finds nothing, which is an error rather than a silent fallback.

TSV format

The transition matrix TSV has 4 data rows (one per reference base: A, C, G, T) and 4 whitespace-separated columns (one per read base: A, C, G, T). The first row is skipped if it is non-numeric (treated as a header). Diagonal values are zeroed automatically. Rows are re-normalized to sum to 1.

A    C    G    T
0.0  0.5  0.3  0.2
0.5  0.0  0.3  0.2
0.4  0.3  0.0  0.3
0.3  0.3  0.4  0.0

Binned quality scores

Modern Illumina platforms (NovaSeq 6000, NextSeq 1000/2000, the simplified MiSeq reporting modes) no longer emit the full continuous Q0–Q40 range. Instead, they quantize each per-base quality into a small set of discrete bins — for example, NovaSeq 6000 emits only {2, 12, 23, 37}. Set binned_quality_bins in the config to mimic this behaviour:

binned_quality_bins: [2, 12, 23, 37]

When this field is set, gen-seq-error-model snaps each observed Q-score from the input FASTQ to its nearest bin (ties round down) before learning seed, transition, and global-frequency counts. The resulting model file is flagged binned_scores: true and lists the bin values as its quality_score_options. gen-reads then samples only those bin values when emitting reads — no extra flag needed on the read-generation side.

Validation rules:

• Bins must be non-empty integers in [0, 94) (entries are sorted and deduped automatically).
• Value 31 is rejected — under Phred+33 it encodes to @, which would corrupt FASTQ output. If you need a bin near Q31, pick 30 or 32.
• Bins with no observed counts in the training FASTQ are still kept in quality_score_options; their transition rows fall back to a uniform distribution, and a warning is logged listing the empty bins.

Q31 (@) and the first quality character

Phred+33 encodes Q31 as @ — the same byte that begins a FASTQ record header. eidolon never writes it as the first character of a quality line, because a reader that scans for @ to find record boundaries would mis-frame the file. (A reader that consumes fixed four-line records, including eidolon’s own, is unaffected.)

What that means depends on the model:

• Binned models reject 31 outright, at config time — see the validation rules above.
• Continuous models may legitimately learn Q31, since real data contains it, and will emit it anywhere in a read except position 1. There, eidolon substitutes the nearest other quality score the model actually contains, with ties going to the lower score. A model learned over the usual Q0–Q40 range therefore substitutes Q30. The rest of the read is untouched.

This is a deliberate, bounded departure from the fitted model: the distribution at position 1 is shifted off Q31 by design, and nowhere else is affected. A model whose only learned score is Q31 has nothing to substitute and is reported as an error rather than emitting a quality line that cannot be read.

Some common platform bin sets (consult your sequencer’s documentation for the authoritative list):

Generating a GC Bias Model

eidolon can learn a GC bias model from a reference FASTA and an aligned BAM file using the eidolon gen-gc-bias-model subcommand. It walks the BAM once to accumulate per-base reference coverage, tiles the reference in fixed-size windows, computes the GC% of each window, looks up the mean coverage over that window, and builds a 101-bin weight table (one bin per integer GC percentage, 0–100%). Each bin’s weight is its median window coverage divided by the median over all fitted windows. The median is used because a few pileup windows (collapsed repeats) can move a bin’s mean arbitrarily far, and a simulator working from a single-copy reference cannot reproduce them. Every fit also writes <output_file>.bins.tsv, listing each bin’s window count, median, mean, weight, and whether the weight was fitted or interpolated. The resulting model can be passed to gen-reads to make fragment start positions favour regions whose GC content matches the coverage bias observed in your data.

$ eidolon gen-gc-bias-model -c gen_gc_bias_model_config.yml

The output is a gzipped JSON model file that can be passed directly to gen-reads via its gc_bias_model config key.

Copy template_config/gen_gc_bias_model.yml to a directory of your choosing and fill in the fields:

# required: reference FASTA used to compute GC content
# Both plain and gzipped (.fa.gz) references are accepted.
reference: /path/to/reference.fa

# required: aligned BAM. Coverage is accumulated in process per reference base.
bam_file: /path/to/aligned.bam

# optional: minimum mapping quality applied while accumulating coverage.
# Default 0 matches `samtools depth` defaults.
min_mapq: 0

# optional: BED file restricting model inference to target regions
# Use "." (dot) to use the entire reference (default)
bed_file: .

# required: output path for the generated model; should end in .json.gz
output_file: /path/to/gc_bias_model.json.gz

# optional: set to true to overwrite an existing output file (default: false)
overwrite_output: false

# Window size in base pairs for GC% and coverage averaging.
# Should match the typical read length (or fragment length for long reads).
# Default: 100
window_size: 100

# Step between successive windows. Must not exceed window_size.
# Use the same value as window_size for non-overlapping windows (recommended).
# Values smaller than window_size produce overlapping windows.
# Default: window_size
window_stride: 100

# Bins with fewer windows than this are not fitted from their own windows: they take the
# weight interpolated between the nearest fitted bins, or the nearest fitted bin's weight
# past either end. Increase this to require more evidence before trusting a bin.
# Default: 10
min_windows_per_bin: 10

BAM filtering: unmapped, secondary, and supplementary records are skipped automatically. Reads with mapping quality <= min_mapq are dropped before contributing to coverage. The default min_mapq: 0 reproduces samtools depth defaults; set it to 20 if your downstream gen-reads runs assume mapq-filtered coverage.

Large genomes: The tool walks the BAM once and allocates a per-contig depth array (u32 per reference base) only for contigs that receive at least one record. Peak memory for hg38 at 30× is on the order of all-contigs-summed depth arrays — see the HPC section below for sizing.

Long reads: Set window_size to approximately your typical read length (e.g. 5000–50000 for ONT or PacBio). Larger windows mean fewer windows per contig and faster model building. Non-overlapping windows (window_stride: window_size) are recommended so each observation is independent.

If the model has no effect in gen-reads: If every GC% bin in your reference has fewer observations than min_windows_per_bin, all weights will be neutral (1.0) and no bias will be applied. Soft-masked (lowercase) reference bases are counted as the bases they are. A BAM that covers only part of the reference, such as an exome or a panel, is refused when the median window has no coverage: set bed_file to the regions it covers. This is logged as a warning. To diagnose, lower min_windows_per_bin or increase the region covered (use a larger reference or remove the BED restriction).

Using the model in gen-reads

# in gen_reads_config.yml
gc_bias_model: /path/to/gc_bias_model.json.gz

# optional: inflate fragment count to compensate for low-weight regions
# true (default): total coverage stays close to the requested depth
# false: coverage in low-GC regions will be lower than requested
gc_bias_normalize_coverage: true

Generating a Fragment Length Model

eidolon can learn a fragment length model from real paired-end alignment data using the eidolon gen-frag-length-model subcommand. It reads a BAM or SAM file, collects the template lengths (TLEN) of confidently-mapped concordant read pairs, trims extreme outliers, and writes a FragmentLengthModel that gen-reads can use.

Since v3.3.0 the model keeps the observed shape by default. Real size-selected libraries are right-skewed — a steep left edge from size selection, and a long tail of larger fragments that escaped it. A normal distribution is symmetric by construction and cannot hold that tail, and the tail is exactly what a paired-end SV caller thresholds against when it decides an insert is “larger than expected”. The builder therefore fits a discrete (empirical) distribution over the observed lengths. Sparse histograms are smoothed so the support has no holes in it, and a model too sparse to have a shape is refused rather than written. Set distribution: normal for the pre-v3.3.0 two-parameter fit, which is still the right choice for small or targeted BAMs (exome, amplicon) where there is not enough data to estimate a shape.

The shipped default — used when a config supplies neither fragment_model nor fragment_mean — is fitted from GIAB HG002 NIST_Illumina_2x250bps (HiSeq 2500, GRCh38, 61.6M read pairs), the same library the shipped quality model comes from: 982 bins over 3–984 bp, mean 408.5, sd 93.2, skew +0.207. Fragment length is set by library chemistry, so treat it as one real library rather than a universal shape: if yours differs, build your own with the command above. Full provenance is in eidolon-core/src/models/model_data/README.md.

$ eidolon gen-frag-length-model -c gen_frag_length_model_config.yml

The output is a gzipped JSON model file that can be passed directly to gen-reads via its fragment_model config key.

Copy template_config/gen_frag_length_model_template.yml to a directory of your choosing and fill in the fields:

# required: path to input BAM or SAM file
# The file must contain paired-end aligned reads. BAM is strongly recommended.
input_file: /path/to/aligned.bam

# required: output path for the generated model; should end in .json.gz
output_file: /path/to/fragment_length_model.json.gz

# optional: `discrete` (default) keeps the measured shape, including the right tail.
# `normal` collapses it to mean + standard deviation, which is the pre-v3.3.0 behaviour
# and is the right choice for sparse input (exome, amplicon, a small targeted BAM).
distribution: discrete

# optional: floor on the TOTAL number of observations, below which the model is refused.
# Default is 2 to handle smaller datasets. Set to 0 to disable filtering entirely.
# NOTE: before v3.3.0 this also DELETED any individual length seen fewer than min_reads
# times, which punched holes in the distribution — hardest in the tail, where counts are
# lowest. Sparse lengths are now handled by smoothing instead.
min_reads: 2

# optional: set to true to overwrite an existing output file (default: false)
overwrite_output: false

Read filtering rules

Only reads that satisfy all of the following are used:

• Paired-end and first in pair
• Mapped (not unmapped, secondary, or supplementary)
• Mate mapped to the same reference sequence
• Mapping quality > 10

Outlier filtering

Fragment lengths that exceed median + 10 × MAD (median absolute deviation) are removed before fitting. This mirrors the filtering in Python NEAT’s fragment length modeler. If no lengths survive the filter, lower min_reads or set it to 0.

Building multiple BAM-derived models in one pass

If you need both a fragment length model and a GC bias model from the same BAM, eidolon gen-bam-models walks the BAM once and feeds both observers from the single pass — half the I/O of running the two per-tool commands back-to-back.

$ eidolon gen-bam-models -c gen_bam_models_config.yml

Copy template_config/gen_bam_models.yml and fill in the sections for whichever models you want:

bam_file: /path/to/aligned.bam
min_mapq: 0   # walker-level filter; FragLengthObserver enforces MAPQ > 10 internally

frag_length:
  output_file: /path/to/frag_length_model.json.gz
  overwrite_output: false
  min_reads: 2

gc_bias:
  reference: /path/to/reference.fa
  output_file: /path/to/gc_bias_model.json.gz
  overwrite_output: false
  bed_file: .
  window_size: 100
  window_stride: 100
  min_windows_per_bin: 10

Both sections are optional but at least one must be present. The output models are interchangeable with what gen-frag-length-model and gen-gc-bias-model would produce on their own, so the resulting files plug into gen-reads exactly the same way.

One filter for the shared walk. The top-level min_mapq is the only MAPQ knob exposed to the unified walker, and it gates both observers. The standalone tools differ: gen-frag-length-model hard-codes its own MAPQ > 10 cutoff (via BamWalkFilter::for_frag_length()), while gen-gc-bias-model honors the min_mapq you pass it. If you want the unified runner’s frag-length output to match the standalone command byte-for-byte, set min_mapq: 10; if you want the GC bias output to match a standalone run with min_mapq: 20, set min_mapq: 20. When the two policies need to differ, run gen-bam-models twice — one config per output — and the per-observer BAM iteration savings still beat running the standalone commands.

Model-builder resource baseline & output fidelity

The model builders (gen-seq-error-model, gen-frag-length-model, gen-gc-bias-model, gen-mut-model, gen-bam-models) are what a first-contact user runs on their own data. eidolon/tests/model_parity.rs pins the fit algorithms on tiny fixtures; this document records how they behave on real, full-size input — the resource envelope and the proof that a fitted model actually shapes gen-reads output.

Harness: scripts/delta/model_builders.sbatch, fed by scripts/delta/stage_soy.sh.

Resource envelope (Delta) — soybean ~978 Mb and human GRCh38 ~3.1 Gb

soy uses two staging modes: chr scopes to the largest contig (fast correctness check); FULL_GENOME=1 stages the whole ~200-contig assembly. human is the full GRCh38 with the GIAB HG002 v4.2.1 truth VCF (stage_hg002.sh). Numbers are wall-clock and peak RSS from /usr/bin/time -v.

Source runs: soy chr = job 19887967, soy full = 19888625, human = 19899988. All overall: PASS.

Memory scales with reference size — plan accordingly

The heavy builders (gc_bias, bam_models) hold the reference plus per-window state, so peak RSS tracks genome size, not read count. Measured:

• ~4 GB peak at ~1 Gb (soybean full genome).
• 11.4 GB peak at ~3.1 Gb (human GRCh38) — the soy→human extrapolation (predicted ~11–13 GB) held; gc_bias landed exactly, bam_models mildly sublinear.

Trivial on an HPC node, but a user building bam_models/gc_bias against a human reference on a 16 GB laptop is near the edge — budget ~12 GB. seq_error is the exception — it streams and stays flat (18–35 MB) regardless of genome or input size.

Nothing approached the harness’s 128 GB / 8 h request at either scale; no OOM, no runtime cliff, clean on both a ~200-contig assembly and full GRCh38.

A real bug the human run caught

The first human run (job 19899126) FAILed mut_model: gen-mut-model rejected the GIAB truth VCF with MalformedVcf("FORMAT list and sample list different lengths"). The VCF is spec-compliant — trailing per-sample FORMAT fields may be dropped (all but GT) — so eidolon’s reader was too strict. Fixed in #364 (read_open_vcf + extract_gt_str pad dropped trailing fields; over-long samples still rejected); the re-run (19899988) passes all five. This is exactly the first-contact-on-real-data failure the harness exists to surface — the H1N1 fixture never exercised a dropped-FORMAT record.

Output fidelity — do the fitted numbers reach the reads?

A model that parses is not the same as a model that shapes the output. Verified links (all are real, runnable tests, not assertions on paper):

The model_builders.sbatch round-trip is only a usability smoke test (reads exist and are correctly-lengthed); the fidelity tests above prove the build→simulate path is numerically faithful through the model-file load path — for all four builders.

Because H1N1’s ~40% GC can’t populate distinct GC bins, the gc_bias test builds a synthetic 2-contig reference (a 20%-GC and an 80%-GC contig) so the weight table is deterministic; the others reuse the bundled H1N1 fixture.

PCAWG SV measurement pass (P2) — results

Status: MEASURED. Every number below was computed from the PCAWG open-tier consensus_sv corpus on 2026-08-02. Reproduction commands are in the appendix.

Why this was run

The SV model’s BND component rested on an unchecked premise: the BND share (23.2%) is the PCAWG per-donor mean of svclass == "TRA" rows, while the generator can only emit same-contig junctions. This pass measures what the source data says, before any code is written, so P3’s scope follows from evidence. See docs/cancer_simulator.md (retraction) and docs/claude_engineering_audit.md §6.

Corpus and coverage

Coverage assertions passed: every row landed in exactly one cell of every cross-tab, and cross-tab totals equal the row count.

M1 — TRA is inter-chromosomal. Exactly.

svclass       n    intra-chrom  inter-chrom   % inter
DEL       86,139       86,139            0      0.0%
DUP       76,214       76,214            0      0.0%
h2hINV    39,495       39,495            0      0.0%
t2tINV    38,851       38,851            0      0.0%
TRA       68,547            0       68,547    100.0%
TOTAL    309,246                                  coverage OK

68,547 of 68,547 TRA rows are inter-chromosomal. Zero exceptions. The four span classes are 100% intra-chromosomal, also with zero exceptions.

This is the load-bearing result. The 23.2% BND share is a translocation rate, and the generator emits none of it — it produces same-contig junctions instead.

There is no such thing as an intra-chromosomal TRA in this corpus. The event class eidolon generates for its BND budget does not exist in the data the budget came from.

M2 — Strand spectrum

svclass       n        ++       +-       -+       --
DEL       86,139     0.0%   100.0%     0.0%     0.0%
DUP       76,214     0.0%     0.0%   100.0%     0.0%
h2hINV    39,495   100.0%     0.0%     0.0%     0.0%
t2tINV    38,851     0.0%     0.0%     0.0%   100.0%
TRA       68,547    24.6%    25.4%    25.3%    24.7%

Known-answer control: passed. Each span class is definitionally one strand orientation, at 100%. That confirms the format reading; had DEL been mixed, nothing else here would be trustworthy.

Two conclusions, and they point in opposite directions:

For inter-chromosomal junctions, uniform geometry is empirically correct. TRA orientation is 24.6 / 25.4 / 25.3 / 24.7 — indistinguishable from uniform. default_bnd_geometry_weights()’s uniform 25%, written as a deliberate non-claim, turns out to be right — for translocations.

For intra-chromosomal junctions, uniform is emphatically wrong, because the four orientations are not one class at all. They are the four separate classes, with distinct rates:

So a same-contig junction with uniformly-drawn orientation is not a BND — it is a DEL, a DUP, or an INV, drawn at the wrong rates, under a SVTYPE=BND label. eidolon already generates DEL/DUP/INV separately, so the intra-contig BND budget is duplicating categories it is also emitting properly elsewhere.

M3 — Translocation partner structure

Since TRA is 100% inter-chromosomal, there is no empirical intra-chromosomal partner distance to fit. The current model’s uniform-over-contig mate draw (sv_model.rs:782-785) models a quantity that does not appear in the source at all.

The real question is which chromosomes pair. They do not pair uniformly:

68,547 TRA junctions across 276 distinct chromosome pairs
(uniform would be ~0.36% per pair)

  chr1  - chr12   2.46%   <- ~7x uniform
  chr12 - chr6    1.73%
  chr12 - chr5    1.32%
  chr1  - chr2    1.21%
  chr2  - chr3    1.11%

endpoint share:  chr1 8.17%, chr12 8.14%, chr6 6.46% ... chr22 1.96%, chrY 0.12%

chr1’s 8.17% tracks its share of genome length; chr12’s 8.14% does not — chr12 is ~4.3% of the genome, so it is roughly 2× enriched. Partner choice is therefore not explained by length alone, and a length-weighted draw would be an improvement over uniform but still wrong.

M4 — Spans vs the fitted log-normal

type       n    source        median      p90      p95       p99       max
DEL    86,139   EMPIRICAL       92kb   16.6Mb   34.4Mb    85.0Mb   236.4Mb
                fitted LN      115kb    8.4Mb   28.3Mb   276.2Mb  unbounded
DUP    76,214   EMPIRICAL      190kb   16.2Mb   35.2Mb    88.6Mb   236.4Mb
                fitted LN      221kb    8.8Mb   25.1Mb   178.0Mb  unbounded
INV    78,346   EMPIRICAL      1.2Mb   36.8Mb   60.1Mb   116.2Mb   237.6Mb
                fitted LN      445kb   55.0Mb  215.3Mb  2786.2Mb  unbounded

(INV empirical here pools h2hINV+t2tINV; the fit uses h2hINV only. Medians are comparable, the count is not.)

This partially corrects a claim in docs/claude_engineering_audit.md §5.6. That section called the fitted tails “biologically implausible.” The tail shape is indeed wrong — the fitted p99 overstates the empirical p99 by 3.2× (DEL) and 24× (INV). But the premise that multi-megabase SVs are implausible was mine, not the data’s: empirical maxima reach ~236 Mb, i.e. whole-chromosome scale. Very large somatic SVs are real.

That inverts the practical concern. The engineering cap (max_length_fraction, default 0.25 of contig length — ~62 Mb on chr1) is not protecting against an unrealistic tail; it is truncating events that genuinely occur. The log-normal body is a fair approximation; the tail is wrong in both directions at once — too heavy in the fit, too tightly clipped in the sampler.

M5 — The shipped model reproduces exactly

Recomputing normalize_pcawg_sv_model.py:108-114 from the corpus:

        recomputed   shipped    delta
  Del       0.2919    0.2919  -0.0000
  Dup       0.2582    0.2582  -0.0000
  Inv       0.1338    0.1338  +0.0000
  Bnd       0.2323    0.2323  -0.0000
  Cnv       0.0559    0.0559  -0.0000
  Ins       0.0279    0.0279  -0.0000

All six match to four decimal places. This had never been verifiable — the .counts.json sidecar was not committed — and it is now confirmed.

Two details that had to be right for this to reproduce, both of which are correct:

• The denominator is all 2,748 donors, not the 2,605 with calls. Donors with zero SVs are real observations; excluding them would inflate every rate by 5.5%.
• t2tINV is deliberately skipped (build_pcawg_sv_vcf.py:143-145) because a balanced inversion appears as two BEDPE rows bracketing one footprint. This is a documented, considered decision, and the corpus independently supports it: h2hINV 39,495 vs t2tINV 38,851, near-equal, consistent with pairing.

The fitting arithmetic was never the defect. The numbers are right; what is wrong is the semantics of the event they are attached to.

M6 — Which BND representation are we modeling?

Not measured — a decision to record, prompted by the observation that BND is a notation for an unresolved adjacency rather than an event class. The corpus makes the choice concrete:

• PCAWG consensus uses svclass to classify every junction it can, reserving TRA for inter-chromosomal ones. Under this convention, “BND” ≡ translocation, and intra-chromosomal junctions are never BND.
• Manta (short read) buckets inversion-oriented adjacencies as BND and calls direct ones DEL/DUP:TANDEM. Delly types the same junctions <INV>. Both observed on identical data, job 20745149.

These are incompatible, and the model currently mixes them: it takes a rate from the PCAWG convention and emits records under something closer to the Manta one.

Proposal for P3: declare bnd_representation: pcawg_consensus explicitly in the model, meaning BND ≡ inter-chromosomal translocation, with intra-chromosomal junctions emitted as DEL/DUP/INV. One convention first; knobs later.

Decision-rule outcome

The rule was fixed in writing before the data was seen:

Consequences for P3:

1. Make sample_variants genome-aware so a BND mate can land on another contig. The read path already supports this (runner.rs:2360); only the sampler does not.
2. Keep uniform BND geometry — measured correct for translocations. The existing non-claim survives contact with the data.
3. Stop emitting intra-contig BND. Under the PCAWG convention it is not a category; those junctions are DEL/DUP/INV, already generated separately at their own rates.
4. Draw partner chromosomes from the empirical pair distribution, not uniformly and not by length alone.
5. Revisit the size cap, which truncates events the corpus shows are real; and the log-normal tail, which is too heavy by 3–24× at p99.

Appendix — reproduction

# fetch (7.7 MB; script verifies md5s)
tools/fetch_pcawg_sv_corpus.sh --out-dir /tmp/pcawg     # or fetch consensus_sv/ only

cd /tmp/pcawg/consensus_sv
mkdir -p icgc tcga
tar xzf final_consensus_sv_bedpe_passonly.icgc.public.tgz -C icgc
tar xzf final_consensus_sv_bedpe_passonly.tcga.public.tgz -C tcga

# flatten: donor, chrom1, start1, chrom2, start2, strand1, strand2, svclass
: > all_sv.tsv
for f in $(find icgc tcga -name '*.bedpe.gz' | sort); do
  zcat "$f" | awk -F'\t' -v d="$(basename "$f" | cut -d. -f1)" -v OFS='\t' \
    'NR==1{next} NF>=11 {print d,$1,$2,$4,$5,$9,$10,$11}' >> all_sv.tsv
done
wc -l all_sv.tsv                     # expect 309246

# M1 — intra vs inter by class
awk -F'\t' '{n[$8]++; if($2==$4) a[$8]++} END{for(c in n) printf "%-8s %7d %7d %6.1f%%\n", c,n[c],a[c]+0,100*(n[c]-(a[c]+0))/n[c]}' all_sv.tsv

# M2 — strand spectrum by class
awk -F'\t' '{n[$8]++; k[$8 FS $6 $7]++} END{for(x in k) print x, k[x]}' all_sv.tsv | sort

# M3 — translocation chromosome pairing
awk -F'\t' '$8=="TRA"{a=$2;b=$4; if(a>b){t=a;a=b;b=t} print a"-"b}' all_sv.tsv | sort | uniq -c | sort -rn | head

# M5 — per-donor means -> type probabilities (denominator = 2748, INV = h2hINV only)
awk -F'\t' '{n[$8]++} END{for(c in n) print c, n[c]}' all_sv.tsv

Filtering Your Data

To run filter-reads, you must first copy the filter_reads_template.yml file from eidolon/template_config/ to a directory of your choosing. Then you can edit the file in your favorite editor. The configuration has four fields:

bed_file: /path/to/my.bed # required

A filename (full path is best) for a bed file with regions to target for filtering. It should be tab-separated in standard bed format.

In configuration:

files_to_filter: [
  # required, in list form
]

Filenames (full path is best) for files to filter by the bed file. For example, if you ran paired-ended fastqs with the vcf, your list would look like:

files_to_filter: [
   /my/home/output/bacteria_r1.fastq.gz,
   /my/home/output/bacteria_r2.fastq.gz,
   /my/home/output/bacteria.vcf.gz,
]

Each file must end in .fastq, .fastq.gz, .vcf or .vcf.gz; any other extension is refused. Gzipped and unzipped files are both read.

filter_key: .

This is the key appended to the filename, before extensions. For example, if your filename is “miscanthus_r1.fastq.gz”, and you set the key as “_only_genes”, the output file would be “miscanthus_r1_only_genes.fastq.gz”. The default is “_filter”. Output is always gzip-compressed and written next to its input, so an unzipped input gains .gz: “miscanthus_r1.fastq” becomes “miscanthus_r1_only_genes.fastq.gz”.

overwrite_output: false # optional, default false

Set to true to allow filter-reads to overwrite existing output files. When false, the run will error if the output file already exists.

Once your files are entered and your config is saved, you can run eidolon:

$ eidolon filter-reads -c my_config.yml

Your output will give you the filenames and success status.

IMPORTANT: This feature only works on fastq files generated by eidolon. It is untested on generic VCF files, but should work. The reason it only works with eidolon-generated fastq files is that eidolon uses a naming scheme that the parser can use to identify where the read is located (enabling filtering), whereas an average fastq will not have that info.

Comparing a downstream caller’s VCF against the golden VCF

After running gen-reads, you can validate a downstream variant caller against the simulated truth using eidolon compare-vcfs. The subcommand classifies every variant into TP / FN / FP, catches denotation-different alternates that exact matching would mis-classify as both an FN and an FP (e.g., a left-aligned indel in the called VCF versus the same indel right-aligned in the golden), and attributes the surviving FNs to specific reasons drawn from the simulator’s configuration.

$ eidolon compare-vcfs -c compare_vcfs_config.yml

Copy template_config/compare_vcfs_template.yml and fill in:

golden_vcf: /path/to/golden.vcf.gz
called_vcf: /path/to/called.vcf.gz
reference:  /path/to/reference.fa
output_dir: /path/to/report_dir
overwrite_output: false

# Optional: restrict comparison to these regions.
target_bed: .
# Optional: comma-separated list. "." treats every called contig as simulated.
contigs_simulated: .
# Optional: BED used by the simulator. Drives FN attribution.
mutation_bed: .
# Optional: TSV remapping BED chrom names (e.g., 1<TAB>chr1).
chrom_aliases: .

# Filters
include_homs: false       # ALT == REF rows
include_filtered: false   # FILTER != PASS / "."

# Equivalence sweep (NEAT 2.1 algorithm, ported)
equivalence_window: 50    # ±N bp around each FP
fast: false               # skip the sweep entirely

# Outputs
write_fp_vcf: false       # also produce FP.vcf

The run produces four files under output_dir:

• comparison_summary.json — schema-versioned (currently 1.3.0), machine-readable. Includes TP/FN/FP totals + per-contig breakdown, precision / recall / F1, the FN attribution roll-up (outside_simulated_contigs, outside_mutation_bed, outside_target_bed, unknown), per-VCF skip counters (multiallelic, homozygous_ref, filtered, outside_target_bed, outside_simulated_contigs, symbolic), and any chrom-naming-mismatch warnings. Symbolic / structural ALTs (<DEL>, <DUP>, <CNV>, …) are byte-incomparable, so they’re counted into skipped.symbolic and excluded from TP/FN/FP classification.
• comparison_summary.txt — same content, human-readable.
• FN_with_reasons.vcf — every surviving FN, annotated with a EIDOLON_REASON INFO tag listing the attribution reasons.
• FP.vcf (optional) — every surviving FP, as-is. Off by default; enable with write_fp_vcf: true.

Equivalence detection. For each false-positive variant, compare-vcfs takes a ±equivalence_window bp window of the reference and applies both the FP set and the FN set within that window. If the resulting byte sequences are identical, the two sets are alternative spellings of the same edit and every consumed FN is promoted to TP. Set fast: true to skip this pass; the report’s totals.equivalents_promoted counts how many TPs were rescued by it.

FN attribution. Each surviving FN is tagged with at least one reason. If the FN’s contig wasn’t in contigs_simulated, that reason is reported alone (the BED checks are skipped because they presuppose the contig was simulated). Otherwise the configured BEDs are checked and unknown is the fallback. Aliases configured via chrom_aliases are applied to BED chrom names at load time, so you can compare against a reference that uses one convention (e.g., chr1) when your BED uses another (e.g., 1). When a BED has any chrom that doesn’t appear in the reference’s FASTA contig list (post-alias), compare-vcfs surfaces a warning in the report — useful for catching single-typo BEDs that would otherwise misattribute every variant on that contig.

Running on HPC

eidolon runs as a single process and fits naturally onto a single HPC compute node. No MPI, distributed computing, or special environment setup is required. The notes below cover resource budgets for whole-genome human-scale runs.

gen-reads

gen-reads is already multi-threaded via rayon (see Parallel Processing above). Set num_threads in your config to match your CPU allocation:

num_threads: 16

Memory scales with num_threads × largest_contig_size, not total genome size — each thread processes one contig at a time. For hg38, chr1 is ~249 Mbp, so 16 simultaneous threads need roughly 16 GB for sequence data plus overhead. Budget additional scratch space for per-contig temp files before the final assembly: roughly 3× the expected output size.

Example SLURM header:

#!/bin/bash
#SBATCH --ntasks=1
#SBATCH --cpus-per-task=16
#SBATCH --mem=48G
#SBATCH --time=2:00:00
#SBATCH --tmp=300G

gen-gc-bias-model

Single-threaded; walks the BAM once and accumulates per-base reference coverage as a Vec<u32> per observed contig. Peak memory is roughly 4 bytes × total reference bases that received at least one record. For a full human genome (~3.2 Gbp) this is ~13 GB. For a single chromosome or a targeted/exome BAM it is much smaller.

If memory is tight on whole-genome runs, split the BAM by chromosome and run gen-gc-bias-model once per chromosome on a representative subset that spans your GC range of interest, then average or combine the resulting models.

Recommended SLURM header (whole-genome BAM):

#SBATCH --ntasks=1
#SBATCH --cpus-per-task=1
#SBATCH --mem=16G
#SBATCH --time=1:00:00

gen-mut-model

Single-threaded; processes the VCF one chromosome at a time. A full human genome WGS VCF (~10 GB) typically completes in 30–60 minutes with peak RSS under 4 GB.

#SBATCH --ntasks=1
#SBATCH --cpus-per-task=1
#SBATCH --mem=8G
#SBATCH --time=1:30:00

gen-seq-error-model

Single-threaded; streams through the FASTQ. For a full genome FASTQ (~600 M reads), expect 30–60 minutes. max_reads caps how many records the fit uses, drawn uniformly from the whole file rather than taken from its start. A FASTQ is written in flowcell order, so its first records all come from one tile of one lane. A capped run still reads every record, once to count them and once to sample, so it shortens the fit but not the pass over the file:

max_reads: 5000000   # a uniform sample of ~5 M records; 0 uses every record

Measured on Delta with eidolon v4.0.0, one CPU, on a 29 GB gzipped HG002 FASTQ: a 5 M-record cap took 11 minutes and 53 MB; fitting every record took 14 minutes. Most of the time is the pass over the file, which is why the cap saves little.

#SBATCH --ntasks=1
#SBATCH --cpus-per-task=1
#SBATCH --mem=4G
#SBATCH --time=1:30:00   # capped or not, every record is read

gen-frag-length-model

Single-threaded; streams TLEN fields from the BAM. A full human genome BAM at 30× typically completes in under 15 minutes with peak RSS under 2 GB.

#SBATCH --ntasks=1
#SBATCH --cpus-per-task=1
#SBATCH --mem=4G
#SBATCH --time=0:30:00

Environment notes

• The eidolon binary has no runtime dependencies beyond a standard C library (glibc), which is present on all Linux HPC systems.
• No module loads or conda environments are required — copy the release binary to your scratch or project directory and run it directly.
• If you compile from source on the cluster, ensure cmake is available (module load cmake on most systems) before running cargo build --release.
• You do not need a Rust module from your administrators. If the site has no rust module, rustup installs a toolchain into $HOME/.cargo as an ordinary user, and it persists across jobs:

  curl --proto '=https' --tlsv1.2 -sSf https://sh.rustup.rs | sh -s -- -y --default-toolchain stable
  source "$HOME/.cargo/env"
  

• If you do build on the cluster, send the artifacts to scratch: export CARGO_TARGET_DIR="$SCRATCH/cargo-target/eidolon". The toolchain itself is modest, but a release build’s target/ directory runs to several GB of crate artifacts against a ~5 MB binary, and an HPC home quota is usually small enough that this is what fills it.

Running eidolon on an HPC cluster

This guide is for users who want to run eidolon at scale — from a bacterial genome on one node up to a whole human (or larger) genome across many nodes — on a SLURM cluster. It distills what we learned validating eidolon on NCSA Delta (see docs/access_report_draft.md §3.6 for the underlying measurements) into cluster-agnostic guidance.

The SLURM scripts referenced here live in scripts/delta/. They carry Delta specifics (module names, an --account, conda env names); adapt those three things to your site and the rest is portable.

1. The one-paragraph mental model

eidolon is single-thread efficient, with low, flat memory. In-process (rayon) threading does scale — but only up to the genome’s contig count (eidolon’s unit of parallelism), so it caps at ~6–9× and does nothing for a single-contig input. The way to use a big machine is therefore many single-threaded processes, each simulating a different region of the genome, packed onto nodes: this parallelizes at sub-contig granularity and, on current builds, packs near-linearly — many more independent shards than threading one process can ever reach. This is region sharding, and it is the HPC story.

Three practical consequences follow, and the rest of this guide is just their detail:

1. Run one thread per process (num_threads: 1) and pack many per node. (Threading a single process is a fine secondary lever for a multi-contig genome run as one job — ~6–9× — but sharding gives far more.)
2. Run on exclusive nodes. On a shared node a co-tenant’s bandwidth-hungry job can inflate per-window time many-fold with no warning; a whole node to yourself is the only way to get reproducible wall-clock (see §5).
3. Pack densely. Independent single-thread shards pack at ~99% efficiency to ~64 per node (one per physical core) and ~75% at 128 — so match K to core count, not the old “keep it low” rule (see §5).

2. Before anything at scale: the log-level footgun

Set --log-level info (the default as of v1.19.1) and never run a large job at debug or trace.

eidolon’s per-base debug!/trace! events fire on the order of coverage × read_length × reference_bp times. At trace, a single chromosome-22 run wrote a multi-gigabyte .neat.log and spent most of its wall-clock on log I/O. Before v1.19.1 the default was trace; if you are on an older build, pass --log-level info explicitly on every run. On any version, only turn up verbosity for a small debugging run, never a production one.

eidolon --log-level info gen-reads -c sim.yml    # good for production
# eidolon --log-level debug ...                  # ONLY on a tiny debugging input

3. Quickstart: a small genome on one node

For anything that fits comfortably in one node’s memory and finishes in a reasonable single-process wall-clock (bacteria, fungi, small invertebrates), you do not need sharding at all. One single-threaded process:

# sim.yml
reference: /path/to/genome.fa
read_len: 250
coverage: 30
ploidy: 2
paired_ended: true
fragment_mean: 350
fragment_st_dev: 50
produce_fastq: true
produce_vcf: true
produce_bam: false
overwrite_output: true
output_dir: /scratch/$USER/run1
output_filename: sample
rng_seed: my run 1
num_threads: 1

#SBATCH --partition=<your-cpu-partition>
#SBATCH --nodes=1
#SBATCH --exclusive          # or a generous --cpus-per-task + --mem if you can't get exclusive
#SBATCH --time=01:00:00
eidolon --log-level info gen-reads -c sim.yml

If you want more throughput on that node, do not raise num_threads — instead run several of these processes concurrently on different regions (see the sharding approach below), or several independent samples at once.

Memory: eidolon’s resident memory is dominated by the reference contigs it loads, not by coverage. Reference loading is target_bed-aware and loads only the contigs a run touches, so peak RSS per process ≈ the size of the largest contig that process simulates plus a small overhead. Budget accordingly: a 25 Mb window on human chr1 still loads all of chr1 (~250 MB); a whole small genome loads all of it.

4. Whole-genome scale: region-sharded array jobs

The genome is partitioned into contiguous anchor-windows. eidolon’s target_bed is a generation-time filter in absolute reference coordinates that anchors each fragment to exactly one window (reads may extend past a window edge; ownership is by anchor point), so the union of all shards reconstructs a whole-genome run with no double-counting and no gaps, and the per-shard outputs merge by simple concatenation. Each window is one single-threaded eidolon process; a SLURM array spreads them across nodes.

The tooling (in scripts/delta/):

Do a fail-fast smoke test first (see §7) so you never commit a big node request to a config that dies in the first shard.

Then the full run:

# 1. index + partition (25 Mb windows → ~124 shards for GRCh38 primary assembly)
samtools faidx REF.fa
N=$(bash scripts/delta/make_shard_beds.sh REF.fa.fai 25000000 $SCRATCH/shards)

# 2. choose packing density K and derive node-tasks (K=64 packs ~99%; see §5)
K=64; T=$(( (N + K - 1) / K ))

# 3. submit the exclusive-node array (%20 caps concurrent nodes — tune to availability)
SHARD_DIR=$SCRATCH/shards OUTROOT=$SCRATCH/wg_run REFERENCE=REF.fa SHARDS_PER_NODE=$K \
  sbatch --array=0-$((T-1))%20 scripts/delta/genome_array.sbatch

# 4. after the array finishes, merge (on a compute node — it's I/O-heavy)
SHARD_OUTROOT=$SCRATCH/wg_run sbatch scripts/delta/merge_shards.sh

For a simulation-only timing run (no downstream alignment/calling — the number users actually want to plan around) this is the whole job: steps 1–3 produce the FASTQ + truth VCF and the array’s time.txt files hold the per-window wall-clock.

5. Choosing packing density (K) and whether to request a reservation

Good news: on current builds this is barely a decision. Independent single-thread shards pack near-linearly — measured on a real 25 Mb human window (30×), one exclusive node runs:

Per-window time is essentially flat to 64 shards/node (one per physical core), dropping only when you oversubscribe hyperthreads at 128. Default to K=64. A whole human genome (306 × 25 Mb windows) is then just ~5 node-tasks, each finishing in ~2 min of compute.

Two knobs:

• K=128 packs the node fullest → fewest node-tasks (~3 for a human genome) at ~75% efficiency. Use it when whole nodes are scarce.
• Lower K only if a genome’s largest contig is big and RAM is tight — peak RSS per node ≈ K × (largest contig loaded), since each shard loads its window’s contig.

Do you need a reservation? Almost never. A whole genome needs only ~3–5 exclusive nodes at K=64–128, which schedules easily even on a busy cluster — you’re not trying to grab 150 nodes for a low-K run. Submit with a %-throttle matched to node availability and let it drain; with no hard deadline, a reservation buys nothing.

Still run exclusive. Dense packing only works because you hold the whole node. A shared per-core slice re-exposes eidolon to cross-tenant memory-bandwidth theft, which is unbounded and makes wall-clock unreproducible — that finding is real and independent of packing density.

Cost: a whole-genome 30× simulation is only ~1.5–5 core-hours of actual compute (soybean ~1.5, human ~4.9). Trivial against any allocation — the binding constraint is wall-clock / node availability, not credits.

6. Reproducibility scope

eidolon is deterministic and parallelism-invariant where it matters, with one documented boundary:

• Determinism: same seed + same reference + same params → byte-identical output.
• Thread-invariant: 1 thread vs N threads → byte-identical. Parallelism never perturbs results.
• Shard-order independent + disjoint: the sharded whole-genome path gives the same merged result regardless of shard completion order, with zero duplicate variants across windows. This is why sharding is safe.
• Contig-order sensitive (by design, monolithic runs only): a single monolithic (non-sharded) run consumes one global RNG stream in contig order, so reversing contig order in the FASTA changes the realization (even the variant count). Consequence: to reproduce a monolithic run exactly, use the byte-identical reference (same contig order), not merely the same genome. The sharded path is immune to this — each shard seeds independently from an absolute-coordinate window. (Tracked as #322: optional per-contig seeding to make even monolithic output contig-order-independent.)

Practical guidance: for reproducible large runs, prefer the sharded path and record your seed(s) + the exact reference file (or its checksum).

7. Fail-fast before you commit big node requests

A large array that dies on a bad config or a malformed-output bug wastes both your time and your queue priority. Two cheap guards:

1. Smoke-test one window first. scripts/delta/wg_smoke.sh runs a single small window as a normal (non-exclusive, few-minute) job and then validates the output the way a downstream aligner would — it checks that every FASTQ record has len(seq) == len(qual) and that reads and truth variants were actually produced. This is the exact class of defect (a malformed FASTQ that wc -l can’t see but bwa-mem2 silently truncates on) that cost us a full at-scale run during adapter validation. Run it whenever you change the build, the config, or the reference:

   REFERENCE=$SCRATCH/neat_data/genome.fa bash scripts/delta/wg_smoke.sh
   

2. Throttle the first real array (--array=...%N with a small N) so a systematic problem surfaces on the first few nodes, not all of them at once.

genome_array.sbatch is idempotent (a completed shard is skipped on retry) and uses --requeue, so a task killed by a sick node reruns elsewhere and only the unfinished windows in its batch redo — you never lose completed work.

8. Shared-filesystem hygiene (Lustre / GPFS)

Large simulation runs move a lot of FASTQ/BAM. On a shared parallel filesystem:

• Stage heavy I/O on node-local disk ($SLURM_TMPDIR / local NVMe) when your site provides it, and copy only the small final artifacts back. A dropped Lustre OST otherwise wedges jobs in uninterruptible I/O wait. (germline_e2e.sbatch demonstrates this pattern.)
• Keep manifests and bookkeeping on durable storage, never scratch — scratch can throw transient EIO or hit quota mid-run, and you don’t want a driver half-submit because a tiny write failed. (run_genome_reps.sh writes its manifest to $HOME.)
• Prune intermediates. Whole-genome 30× FASTQ is ~35–40 GB per run; delete BAMs and FASTQ after scoring if you only need the metrics (PRUNE=1 in the harnesses).
• Watch your scratch quota — on Delta it is a project quota, not per-user, so a teammate’s data counts against yours.

9. Troubleshooting

Measurements in this guide come from the eidolon validation campaign on NCSA Delta; see docs/access_report_draft.md (§3.6 for scaling, §3.10 for reproducibility, §3.11 for the malformed-FASTQ story) and docs/regression_protocol.md for how the performance and fidelity baselines are frozen and gated.

eidolon re-validation / regression protocol

When a feature changes eidolon (e.g. adapter modeling, #125), this protocol confirms the change does not:

• (M) blow up memory — peak RSS regression,
• (S) slow down — wall-clock regression,
• (D) emit errant data — fidelity or determinism regression.

It reuses the Delta validation harnesses already built, and gates a candidate build against a frozen baseline with statistically-grounded tolerances.

Why this is now possible: the replication study gives the noise floor

The Phase-2 replication runs measured the run-to-run standard deviation of every metric (see replication_audit.md). That sd is exactly what separates a real regression from normal scatter:

A gate fires when a candidate moves more than k·sd below baseline (default k = 2 ≈ 95 % confidence). Without the reps you could not set these — you’d be guessing whether a 0.97→0.94 dip is a bug or noise.

The suite: failure mode → metric → harness → gate

All harnesses already emit the needed numbers via the collect_* scripts and archive_run (which also records git SHA + core-hours for provenance).

Tiers — run the cheapest tier that covers the change

Tier 0 alone catches the big three failure modes grossly (crash, memory blowup, non-determinism, fidelity collapse). Tiers 1–2 add resolution and breadth.

Feature-touch matrix — which extra checks a change demands

Worked example — adapter modeling (#125)

Adapters add read-level sequence (3′ adapter readthrough, etc.). Two question sets:

Regression (must NOT change):

• Determinism still byte-identical (exact).
• Wall-clock / peak RSS within tolerance — adapter logic must not inflate either.
• Truth VCF unchanged — adapters are a read artifact, not variants (md5 of truth VCF vs no-adapter run must match).
• Post-trim fidelity == baseline: simulate with adapters → trim (cutadapt / Trim Galore) → align → call; SNP/indel recall/precision must match the no-adapter baseline within 2 sd. Adapters that survive trimming and corrupt calls would show up here.

New-feature checks (adapter-specific) — built as a Tier-2 seam: run_adapter_validation.sh (CONTROL=1) runs a four-arm matrix that isolates the adapter effect from the short-insert coverage effect (off rejects short inserts, so a naïve off-vs-on comparison is confounded):

The callability check rests on the adapter-only contrasts (on_trim − short_ctrl, on_raw − on_trim), not the confounded off-vs-on. collect_adapter_validation.sh with CANDIDATE_TSV=… emits adapter_* metrics that regression_gate.sh checks against the adapter_* rows in baseline_metrics.tsv — gating callability (recall stays high → catches a malformed-FASTQ/collapse regression), readthrough presence (on_raw soft-clip elevated), trim recovery, and the null adapter effect.

# after the 12 jobs finish:
CANDIDATE_TSV=$RESULTS_DIR/adapter_candidate.tsv \
  bash scripts/delta/collect_adapter_validation.sh $RESULTS_DIR/adapter_validation.manifest
bash scripts/delta/regression_gate.sh scripts/delta/baseline_metrics.tsv \
  $RESULTS_DIR/adapter_candidate.tsv 2

Result (#125): adapter readthrough is detectable and fully callable — the pure adapter effect is within replication noise; the only fidelity difference vs the long-insert baseline is the short-insert coverage effect, present with adapters entirely absent (short_ctrl). Lesson: the collapse bug that motivated this seam (a malformed FASTQ record) was invisible to zcat/wc read counts and only surfaced through a real aligner — so the gate keys on aligned/called metrics, not read counts.

Tier plan for adapters: Tier 0/1 (default path, adapters off) on every commit / PR; run the Tier-2 adapter seam above when the change touches the adapter or fragment-generation code, or pre-release.

Baseline management

• Freeze baseline metrics on develop into a committed baseline_metrics.tsv (one row per metric: value, sd, n, git SHA, date). Re-freeze on each release.
• A candidate run produces candidate_metrics.tsv in the same shape.
• regression_gate.sh diffs the two, applies the gates + tolerances above, prints a PASS/FAIL table, and exits non-zero on any FAIL (CI / --dependency-friendly).

Re-freezing PERF rows without disturbing fidelity (the v1.19.1 logging fix)

A change can be output-preserving but performance-changing — the v1.19.1 logging fix (#340) is the canonical case. Every Phase-1/2 run used the old default log level trace, whose per-base debug events fire ~coverage × read_len × ref_bp times and burned most of the wall-clock on log I/O; so the frozen *_wall_s / *_rss_mb rows were captured under that handicap and are now stale (the fixed build is faster). The fix does not change output bytes, so every fidelity / determinism / recall row stays valid — only the perf rows must be re-measured and re-frozen. Two steps:

1. Prove the change is output-preserving (before trusting the split above). Run baseline_capture.sbatch on the fixed build and compare its baseline.json fastq_md5 / vcf_records_md5 against a pre-fix run’s:

   bash scripts/delta/rebaseline_perf.sh --check-identity OLD.json NEW.json
   

   If this is not PASS, the change altered output and the fidelity rows must be re-validated too — it is not a pure perf/logging change.

2. Re-freeze only the perf rows from a fresh benchmark on the fixed build:

   MODE=submit TIER=0 LABEL=logfix bash scripts/delta/regression_suite.sh   # TIER=1 also re-freezes chr22_*
   # …when the perf job finishes…
   MODE=collect MANIFEST=$HOME/regression_logfix.manifest \
     bash scripts/delta/regression_suite.sh > $HOME/logfix.candidate.tsv
   bash scripts/delta/rebaseline_perf.sh $HOME/logfix.candidate.tsv          # rewrites *_wall_s / *_rss_mb only
   

   rebaseline_perf.sh touches only rows ending in _wall_s / _rss_mb (preserving their gate + tier), writes a .bak, and leaves every other row untouched — so a stray fidelity value in the candidate can never overwrite a fidelity baseline. Commit the result with the fixed build’s git SHA.

Automation (built)

Three files implement this, thin wrappers over the existing harnesses — no new simulation logic:

• scripts/delta/baseline_metrics.tsv — the frozen baseline, seeded from the Phase-1/2 validation + replication numbers (value, sd, gate, tier). Re-freeze on each release by running the collect step on develop.
• scripts/delta/regression_suite.sh — MODE=submit TIER=n submits the tier’s jobs (order-independence + germline + perf, plus chr22 germline/cancer at Tier 1) and writes a manifest; MODE=collect MANIFEST=… extracts metrics into a candidate TSV. (Tier 2 — variety, SV per-type reps, sweeps — is driven by run_input_variety.sh / run_cancer_sweeps.sh REPS=3 and appended.)
• scripts/delta/regression_gate.sh baseline candidate [max_tier] — applies the gates, prints a PASS/FAIL table, and exits non-zero on any FAIL (CI-friendly).

One command (fire-and-forget) on a feature branch, from repo root:

MODE=run TIER=1 LABEL=adapters bash scripts/delta/regression_suite.sh

This submits the tier’s jobs plus a dependency gate job (--dependency=afterany) that auto-collects and writes a PASS/FAIL verdict to $RESULTS_DIR/regression_<label>.verdict.txt when they finish — no babysitting. A crashed harness surfaces as a FAIL (non-numeric), not a silent pass.

Manual three-step (if you want to inspect between stages):

MODE=submit  TIER=1 LABEL=cand bash scripts/delta/regression_suite.sh
MODE=collect MANIFEST=<path>   bash scripts/delta/regression_suite.sh > candidate.tsv
bash scripts/delta/regression_gate.sh scripts/delta/baseline_metrics.tsv candidate.tsv 1

There is no automatic trigger (git hook / CI) — GitHub runners can’t reach Delta, so a run is kicked off deliberately (the convention: MODE=run TIER=1 on the feature branch before merging; TIER=2 for SV/perf-touching changes).

The gate logic and all collect parsers are unit-tested; the sbatch submit path runs on Delta. Capturing a fresh baseline is just the same collect step pointed at a develop run.

Versioning and the public API

As of v3.0.0, eidolon follows Semantic Versioning 2.0.0: a MAJOR bump means something you may depend on changed incompatibly, a MINOR bump adds functionality compatibly, and a PATCH bump is fixes only. Releases before v3.0.0 were versioned less strictly — notably v1.11.0 and v1.12.0 changed the FASTQ read-name format in minor bumps.

SemVer requires a project to say what its public API is. For eidolon:

Public API — a change here means a MAJOR bump

• Names and semantics of emitted VCF INFO tags, the VCF sample-column name, and the FILTER / FORMAT fields eidolon writes
• FASTQ/BAM read-name (QNAME) format — the EIDOLON_generated_ / EIDOLON_chimeric_ prefixes and the positional fields encoded after them
• CLI subcommand names, flags, and configuration-YAML keys
• Model-file compatibility, in the backward direction only — a model built by an older eidolon stays readable by a newer one. See Model files below for what that does and does not cover

Not public API — may change in a MINOR or PATCH release

• The eidolon-core Rust library surface; it exists to serve the binary
• Log output, progress reporting, and human-readable messages
• The exact simulated content for a given seed — reads are a random draw, so sampler and model changes legitimately alter output while preserving the format

Versioned separately from eidolon itself

• The compare-vcfs JSON report carries its own schema_version field (currently 1.3.0), bumped only when a backward-incompatible field change lands. Read that field rather than inferring the report’s shape from the eidolon version — the two move independently, so a minor eidolon release can carry a new schema, and a major one need not

If you parse eidolon’s output, the first list is what you are relying on, and a major version bump is your signal to check this file before upgrading.

Model files

“Compatible” covers one direction, and neither direction is currently checked automatically.

• Old model, newer eidolon — supported. Fields added since the model was written deserialize to their shipped defaults, so the model keeps its previous behavior for anything it does not carry, rather than failing to load.
• New model, older eidolon — NOT supported, and NOT detected. Unrecognized fields are ignored silently. An older binary will read a newer model, discard what it does not understand, and simulate using its own defaults with no warning that it did so. Do not move model files backwards across releases.
• Model files carry no version stamp. A model cannot tell you which eidolon wrote it, so neither direction above can be verified — by you or by eidolon. Tracked in #708.

The practical rule: rebuild models with the eidolon you intend to run, and treat a model file as belonging to the version that produced it.

SV support matrix — what eidolon can reproduce, and from where

Measured 2026-08-02 on the H1N1 fixture (8 contigs), 60× diploid, ploidy: 2, heterozygous unless noted. Every row is an observation, not a design intent. Cells marked ⚠ or ✗ are known gaps, not omissions.

Two separate questions, and they are not required to agree:

• Input — you supply the record via input_vcf. The contract is fidelity: what you put in comes out, in the truth VCF and in the reads.
• De novo — eidolon samples the record from a model. The contract is realism: the distribution matches the source corpus.

Input fidelity is the higher priority. A researcher modelling a specific variant needs to know it reached the output; a mismatched distribution is a modelling argument, but a variant that silently fails to appear is a broken tool.

2026-08-16 characterization addendum

The focused Gate 2b checks and the first Delta performance/BAM measurements are now complete. These are baselines for the current develop build, not pooled recall estimates.

Post-#537 subclonal SV smoke (develop 4a3ea4a)

Delta jobs 21281985 and 21282371 exercised the merged subclone-SV implementation on multi-contig S. cerevisiae. The higher-rate run produced 12 SV records (7 at CCF 0.2, 1 at CCF 0.5, and 4 clonal at CCF 1.0); all 12 carried EIDOLON_CCF and EIDOLON_VAF. Three reciprocal translocation pairs shared CCFs of 0.2, 1.0, and 0.2. After bwa-mem2 alignment, 199 chimeric alignments were retained, 100% of 2,426,389 reads mapped, 69 were supplementary, and 99.99% were properly paired. This validates CCF stamping, reciprocal-BND pairing, and read-level junction survival; it is not a caller recall estimate.

The full-GRCh38 post-merge stress run (job 21282504, 30×, purity 0.6, nominal SV_RATE_SCALE=1.0) produced 140 tagged SV records: 31 at CCF 0.2, 43 at CCF 0.5, and 66 at CCF 1.0. Its merged BAM contained 588,995,182 reads, all mapped, with 2,837 chimeric alignments, 2,216 supplementary records, 99.99% proper pairing, and 71 GB compressed size. This confirms the subclonal SV path at primary-assembly scale; it is still a simulation/read-level result rather than a pooled caller-recall estimate.

Local performance and BAM baseline

The FASTQ-only benchmark (Delta job 21188169, three replicates per row) remained consistent with the earlier chr22 measurement after the BAM changes:

The yeast thread sweep was 22.85 s (1 thread), 11.15 s (4), 7.18 s (16), and 7.68 s (64). This benchmark does not produce BAMs. A separate chr22 BAM timing smoke produced 7,780,196 records, including 60 chimeric records, in 3:05.39 at 230.5 MB peak RSS; samtools quickcheck passed and the BAM was 843 MB. The result is a BAM-path baseline, not a before/after comparison against an older binary.

A full-primary-assembly BAM-only run (Delta job 21262840, GRCh38 at 30×, 128 requested threads) completed successfully in 23:56.16. It wrote 583,088,206 records to a 62 GB BAM; samtools quickcheck passed. Process-level /usr/bin/time reported 13.6 GiB peak RSS and 924% CPU utilization (about 9.2 cores), showing that this path does not scale to all 128 requested threads. Slurm reported a much larger 111 GiB MaxRSS; that accounting value is retained as a cluster-side note, while the process-level measurement is the comparable RSS figure. This run produced no input SVs, so it characterizes full-genome BAM throughput and resource use rather than BND recovery.

First multi-contig GRCh38 SV smoke

Delta job 21170801 (30×, purity 0.6, SV_RATE_SCALE=1.0, Manta only) produced 137 somatic SV records: 29 DEL, 32 DUP, 15 INV, 2 INS, 52 BND, and 7 CNV. All 52 BND records parsed, and the scorer controls passed (truth-vs-truth recall 1.000; shifted-truth decoy recall 0.000). There were no intra-contig BND spans because this smoke’s BND set was inter-contig, so BNDspan scoring was correctly skipped.

These are PASS-only caller measurements from one replicate, not simulator correctness verdicts. Most importantly, the independent INS probe found 2/2 planted insertions in the reads even though Manta emitted no PASS INS calls. Five de novo insertions longer than the 151 bp read length were dropped as expected under #516. The run retained 30 GB normal and 70 GB tumor BAMs with PRUNE_BAM=0 for follow-up inspection.

Follow-up characterization is complete. The initial smoke is retained above as a historical single-run baseline; the follow-up observations below provide the caller comparison and the nominal-rate Manta repeatability check. These are still separate from the larger pooled campaign reported in the whole-genome tier below.

Follow-up multi-contig GRCh38 observations

Delta job 21226562 reran the same 30×/0.6/SV_RATE_SCALE=1.0 design with both Manta and Delly on one simulated tumor/normal pair. The truth again contained 137 somatic records (including 52 BNDs), all 52 BNDs parsed with zero unpaired or mispaired junctions, and all scorer controls passed. Manta recovered 42/52 BNDs (recall 0.808, precision 0.457); Delly recovered 20/52 (recall 0.385, precision 1.000). Overall scoreable recall/precision were 0.795/0.747 for Manta and 0.821/0.660 for Delly. The independent insertion probe found both planted insertions in the reads (2/2); the zero manta_INS score is therefore a PASS-filter limitation, not evidence that the simulator omitted the insertions.

Delta job 21249416 supplied the third nominal-rate Manta observation. It reproduced the same truth counts and clean BND geometry (52/52 parsed; 0 unpaired or mispaired), with scorer self-tests and decoys passing. Manta scored 67 TP, 11 FN, and 20 FP overall (recall 0.859, precision 0.770), and 42 TP, 10 FN, and 50 FP for BND (recall 0.808, precision 0.457). Across the three Manta observations, BND recall was approximately 0.77–0.81 while aggregate recall remained about 0.86. The consistent truth-side checks and the caller-to-caller difference in BND recovery support a caller/representation limitation rather than a simulator defect.

Input VCF → output

Untested cells: none remain

Every row above is measured. Read-level evidence is depth against a separate no-variant control run, or BAM CIGAR I/D counts against that control’s sequencing-error background (15 I, 20 D in the probed window), or chimeric-read counts.

The “~8% high” figure is size-dependent, and the size it was measured at is not realistic. #499 measured a 1.2 kb event on the 13.5 kb H1N1 fixture, a window where depth noise alone is ~13% sigma. Re-measured on a 10 kb event over a 1 Mb reference, depth reads 0.98 +/- 1.9% — accurate. The three rows above therefore report a bias that does not survive at realistic event size, and are kept only because they are what H1N1 actually shows.

The residual mechanism is real and tracked as #584: junction reads are emitted on top of coverage-multiplied reads rather than as part of them (+0.52 at multiplier 1.0 on a 300 bp DUP). Their count does not vary with event length, so the excess as a fraction of the event scales roughly as 1/length — negligible at 10 kb, dominant at 1 kb. Anything calibrating against eidolon depth needs this number at its own event size, not this table’s.

The confirmed failures

[#498] BND with inserted sequence — FIXED. VCF 4.2 allows a breakend ALT to carry novel bases inserted at the junction, common in real rearrangements since NHEJ frequently leaves sequence at the breakpoint. parse_bnd_alt accepted it while get_bnd_pieces rebuilt the read from reference pieces only, so the insert never reached the output while the truth VCF kept it — a benchmark asserting an insertion the reads did not contain.

parse_bnd_alt now returns the inserted bases and get_stitched_sequence splices them between the two reference pieces. Two details are load-bearing and each has its own test:

• Which end of the replacement sequence holds the reference base flips with the form — t[p[ and t]p] start with it, [p[t and ]p]t end with it. Getting it backwards splices the anchor into the read and drops one inserted base.
• The insert comes OUT of the fragment budget (L1 + insert + L2 == frag_len), it does not extend it. The inserted bases consume READ bases but no REFERENCE bases, so a correct pair spans less reference than the library mean — a reduced insert size is one of the signatures a caller uses to detect an insertion at a junction. Letting the fragment grow keeps the span at the library mean, erases that signature, and puts the mate further into the second piece than the fragment distribution allows: systematically discordant pairs at every inserted junction, and an inferred event size that does not match what was planted. It does not change depth — R1 and R2 are read_len each however long the fragment is. Invisible to every read-content check, since the insert is present, adjacent and correctly oriented; caught only by mutating the arithmetic.

The insert is never reverse-complemented: it is written in the orientation of the record declaring it, while the rev flags describe only how each reference piece is read.

[#500] Symbolic <INS> — FIXED. SVTYPE=INS;SVLEN=60 used to be preserved in the truth VCF while the reads came out behaviourally identical to the no-variant control — same I count, same D count, same depth. A 60 bp insertion was requested and nothing was inserted, so a benchmark built from that truth scored a caller as having missed an insertion that was never in the data.

A symbolic <INS> from input_vcf is now realized: SVLEN novel bases are drawn by the same sample_novel_insertion_bases the de novo path uses — one model must not plant two different kinds of insertion depending on which door the record came through — and the record becomes a literal ALT, which is what routes the bases through the insertion machinery. The truth VCF then states what was actually planted rather than a size and a promise. #580 tracks making that sampler mobile-element-aware; both paths inherit it at once. A <INS> with no usable SVLEN has no length to synthesise and is dropped with a warning rather than kept as an unachievable truth record.

The de novo path did synthesise novel sequence for INS all along, so the capability existed; it was simply not reached from input_vcf. That last sentence was an assertion when first written and is now measured: jobs 20885875/20892075 planted 4 insertions of 296–2500 bp whose novel bases are present in the truth VCF, at a rate consistent with the model’s Ins = 0.0279 (4 observed against ~5.2 expected-visible). It was nearly retracted as false in 2026-08 on the strength of truth INS: 0 in the SV harness — which turned out to be measuring the harness’s own selector, not the generator. See §“De novo generation” below and docs/claude_engineering_audit.md §5.3.

[#516] Large insertions were only PARTIALLY realized — FIXED. Kept here because the defect’s shape is the most instructive one in this repo: the head was always present while the interior collapsed to zero, so any probe near the insertion’s start said it worked, which is where a casual check looks. The truth VCF meanwhile declared the full SVLEN.

Measured on the H1N1 fixture at read_len=100, fragment_mean=250, GT 0/1, counting reads containing a 30-mer from the head / middle / far end of the inserted sequence. Identical fixture and config before and after, so the columns are directly comparable:

The fix samples fragments in altered-haplotype coordinates, where an insertion has real width, so fragments begin inside it. eidolon/tests/long_insertion_rework.rs pins eight criteria on eidolon’s own FASTQ, and eidolon/tests/gate2_realigned_ins.rs pins the same insertion through bwa-mem2: a homozygous 600 bp insertion at 60x yields 5/5 interior probes against 0/5 in a same-seed no-variant control, 96 breakpoint clips against 0, and 51 reads unmapped-with-mapped-mate against 0.

At genome scale, de novo on a 5 Mb chr22 window (default germline model, sv_rate_scale=8): 427 of 440 insertions ≥200 bp have their middle sequence in the reads, with zero fallbacks. The remaining 13 are thinly covered rather than absent — mean interior depth 3.96 against 6.77 for the rest, still 68% of their length covered. Interior coverage of a heterozygous insertion runs at ~half the flanking reference (7.0 vs 13.0 reads per 30-mer), which is the correct ratio, not a shortfall.

What is NOT yet established: Delta job 21378484 (chr22, SV_RATE_SCALE=30, COSMIC pancancer somatic) verified 4 of 4 planted insertions at 5/5 interior probes including one of 2127 bp — real evidence, but a denominator of four. verify_planted_ins now refuses to let a run that small read as validation. A larger-denominator campaign is outstanding.

The historical mechanism, kept because it explains why the defect was invisible:

Fragments and read windows were chosen purely in reference offsets (cover_dataset, generate_fragments.rs:323), and an insertion has zero reference width, so no read window could ever begin inside one. Reads are assembled per-read by walking a reference slice and expanding variants inline, bounded by bases_written < read_length (fastq_tools.rs:429) with a break 'outer at :555 that discards the rest of the insert silently — no log, no counter. Hence:

novel bases visible in one read = min(L, read_length − anchor_offset − 1)

At most read_length − 1 inserted bases are ever realized, at any declared SVLEN. Probe visibility follows exactly: head needs anchor_offset ≤ 69 (no L term), middle needs anchor_offset ≤ 84 − L/2 (impossible for L ≥ 170), tail needs ≤ 99 − L (impossible for L ≥ 100) — which reproduces all fifteen measured cells.

The head count is size-independent: with the heterozygous coin removed (GT 1/1) it is 42/41/50/50/50 across the five sizes. An earlier version of this note claimed it “rises with size” and read meaning into 16 → 27; that was RNG drift, since a larger insert fires break 'outer sooner and consumes fewer sequencing-error draws.

Consequence, measured: SV validation campaign 20925151 planted 22 de novo somatic insertions of 61–2155 bp and Manta detected exactly one — the 61 bp event, and only at MinSomaticScore. That is not primarily a caller limitation. Manta documents a fully assembled ceiling of “approximately twice the read-pair fragment size” and states it also reports very large insertions from a breakend signature alone, as IMPRECISE <INS> with LEFT_SVINSSEQ/RIGHT_SVINSSEQ. Most of the planted set was inside that range, and Manta reported nothing within 2 kb of any of them — because the evidence for the declared length was not in the reads. INS recall figures from campaigns before the fix are invalid above ~a read length and must not be re-quoted; DEL/DUP/INV/BND/CNV from those runs are unaffected. Job 21378484 is the first post-fix counter-example: Manta’s single INS call was an exact position-and-size match to a planted 93 bp insertion, dropped only by --passonly as MinSomaticScore, while the three larger insertions it did not call were all confirmed present in the reads at 5/5 interior probes. Evidence present, caller declined — which is a different finding from evidence absent, and the read-level probe is what separates them.

The interim read_len - 1 cap on de novo insertions has been REMOVED. While #516 was open, gen-reads refused to plant a de novo insertion whose novel sequence exceeded what a read could carry, rather than emit a truth record the reads could not support. That truncated the model’s own distribution — the bundled default puts Ins at log-normal (5.7, 1.0), median ~299 bp, so at read_len=150 roughly three quarters of INS draws were refused.

With insertions sampled in altered-haplotype coordinates the cap is unnecessary, and its removal is pinned by eidolon/tests/multi_sv_integration.rs, which now asserts the opposite of what it once did: a de novo insertion MAY declare more novel bases than read_len - 1. Reinstating the cap, or unwiring the sampler, fails that assertion while every other assertion in the file stays green — which is the specific regression it exists to catch. De novo insertion sizes now follow the model rather than the read length; measured largest realized on a 5 Mb chr22 window is 9,093 bp.

Pinned by runner.rs’s unrealizable_insertions_are_dropped_at_exactly_the_boundary and the_insertion_cap_does_not_touch_other_variant_types (a must-not-fire: dropping large DELs/DUPs/INVs would be far worse than the bug being worked around), plus an end-to-end assertion in multi_sv_integration.rs that no de novo INS in the truth VCF exceeds read_len - 1. Removing the call site makes that fail with eleven oversized records.

This is an interim measure. It removes the false-truth problem, not the capability gap: eidolon still cannot simulate insertions longer than a read.

[#500] A single breakend (A.) — FIXED by rejecting it. VCF 4.2 allows an unresolved partner. parse_bnd_alt returns no mate, and accepting the record anyway was worse than ignoring it: depth fell 72.0 → 42.4 (−41%) with zero chimeric reads. The truth declared a breakend, the reads contained no junction, and a depth caller saw a partial deletion no record described.

It is now rejected during input_vcf filtering, with a warning naming the record. That keeps it out of the truth VCF and leaves coverage alone — the regression test asserts both, and the depth ratio it checks is the old failing threshold inverted (0.59 measured then, > 0.85 required now). Rejection rather than realization because eidolon has no way to build a junction read without a partner; revisit if junction generation ever represents one.

Fragment spanning a whole DUP — retracted, it was never broken. This cell previously read ❌, citing a three-piece stitch (left + block + right) that get_dup_pieces cannot express, and “4 of 30 reads match no haplotype”. All three parts are withdrawn:

• Measured false. short_dup_spanned_by_a_fragment_still_matches_the_derived_haplotype sweeps DUP blocks of 60/100/150/200/400/800 bp against a ~200 bp fragment — every one spannable by a fragment — and all match the derived haplotype at the ≥90% threshold used throughout that file.
• The 13% was the error model. matches_with_one_small_indel exists because a purely positional comparison frameshifts on sequencing-error indels; its own note says that without it “a correct implementation looks ~9% broken”. 4 of 30 = 13% sits inside that band, so the original measurement almost certainly predates that tolerance.
• The issue reference was wrong. #474 is closed and concerns INV/DUP anchor-base conventions, not stitching — it never covered this.

Residual, and why this is ✅ rather than settled: the test asserts the reads that exist are correct, not that none are missing. Fragments needing three pieces being silently dropped would dip coverage across a short DUP and leave the test green. A depth comparison against a no-variant control would close that.

<INV> — investigated, and the test was wrong, not the code

The first measurement of this cell read 0.63, which looked like a balanced inversion losing a third of its coverage. It was flagged rather than filed, and that was the right call: diagnosis showed the fixture was at fault.

Against a no-variant control run of the identical config:

300bp inversion, 250bp fragments (H1N1_NA:300-600)
  inside inversion      0.62
  flank left            0.81
  flank right           0.76
  far away              1.00     <- coverage is fine elsewhere, so the loss is real

1200bp inversion, 250bp fragments (H1N1_PB2:500-1700)
  deep interior         1.00     <- NO loss
  near left junction    0.82
  near right junction   0.74
  outside               0.99

The inverted sequence is covered at full depth. The depletion is junction-proximal, extending roughly one fragment length around each breakpoint. A 300 bp inversion is smaller than a 250 bp fragment, so nearly every fragment overlapping it touches a junction and the entire event sits inside the dip zone — which is what produced 0.63.

The remaining junction dip (0.74–0.82 over ~180 bp) is consistent with the mechanism: a fragment spanning a junction becomes a chimeric read placed at one locus, so its coverage contribution to the other side is not counted. Real aligners soft-clip and place such reads at one side too, so some dip is expected. How much dip is realistic has not been compared against real data — that is unquantified, not verified.

Practical consequence for anyone writing SV tests: an inversion shorter than the fragment length is not a clean fixture. Its coverage behaviour is dominated by junction effects, not by the inversion.

<CNV> — semantics confirmed, but a systematic ~8% bias found (#499)

The first reading of this cell (1.67 for CN=4) matched no sensible interpretation. It was measured against an in-run control span on a 300 bp event — the same flawed method that produced the bogus <INV> result. Re-measured properly (1200 bp event, deep interior only, against a separate no-variant run):

CN   multiplier  expected  observed  obs/exp
 0      0.00       0.00      0.00     exact
 1      0.50       0.50      0.54     1.08
 2      1.00       1.00      1.00     exact
 3      1.50       1.50      1.61     1.07
 4      2.00       2.00      2.17     1.09
 6      3.00       3.00      3.28     1.09

Semantics answered: CN is the total copy number, multiplier CN/ploidy, and genotype is ignored when CN is present. That is consistent and defensible.

But the delivery is ~8% high whenever the multiplier is neither 0 nor 1 — and it is not CNV-specific. The same bias appears in DEL and DUP, both genotypes:

DEL het  0.54 vs 0.50   (1.079)      DUP het  1.61 vs 1.50   (1.073)
DEL hom  0.00 vs 0.00   exact        DUP hom  2.17 vs 2.00   (1.084)

coverage_multiplier_for returns the correct values; the gap is between the requested multiplier and the depth actually delivered. Mechanism not diagnosed — see #499.

Update 2026-08-23 — mechanism diagnosed, and it is size-dependent. #499 is closed. Two effects contribute, both boundary effects, both shrinking as the event grows relative to the fragment length: chimeric junction reads landing outside the coverage-multiplied budget (#499’s own, excess ~1/event_length), and fragments legitimately extending across the event’s coverage-multiplier boundary into its neighbours (introduced with the fragment-placement rewrite, and the thing that removes an artificial dead zone at every such boundary). See the size guidance below.

How big does an SV need to be for its declared depth to be accurate?

Rule of thumb: depth multipliers land within ~1% for events ≳10× the fragment length. Below ~5×, expect several percent of under-delivery. On a small contig it is worse, because the flanking baseline is narrow too.

Measured, het <DUP> (declared 1.50×), read_len: 100, fragment_mean: 250, 3 seeds averaged except where noted:

Read this conservatively: the 1,200 bp (93%) against 1,500 bp (~100%) difference is within seed-to-seed noise at that scale, so the useful signal is the trend, not any single row. What the table does establish is that both ends are real — a fragment-length-scale event genuinely loses several percent, and a realistic event genuinely does not.

Practical guidance:

• Validating depth-multiplier accuracy (DEL/DUP/CNV dosage, copy-ratio callers): use events ≳10× fragment length on a reference with wide flanking sequence. The depth_modulation_is_accurate_at_realistic_scale test does exactly this and is the reference measurement.
• Small references (viral genomes, small scaffolds, the H1N1 fixture) cannot host an event much larger than the fragment length and leave wide flanks — H1N1 is 2,280 bp total. They remain excellent for fast mechanism checks and are used that way throughout this file, but they should not be used to validate depth precision.
• Small SVs on any reference carry the same few-percent under-delivery. If you are simulating sub-kilobase events and depend on their exact dosage, measure it rather than assuming the declared multiplier.

Whole-genome tier — de novo, scored by callers

Everything above is H1N1: 8 contigs of 1–2.3 kb, largest event ~1.2 kb. That fixture cannot say whether an SV behaves at the sizes anyone simulates. This section is the other tier — GRCh38 at 30×, purity 0.6, SV_RATE_SCALE=1.0, scored by Manta 1.6.0 and Delly against truvari.

It answers a different question. The tables above are input fidelity — what you supplied came out. This is de novo generation as an independent tool can recover it. Neither substitutes for the other, and the gap between them is the subject of docs/sv_polish_roadmap.md.

⚠ Four replicates, pooled — against a target of ~8. aggregate_sv_reps.sh pools summed TP/FN/FP (not a mean of per-replicate ratios, which would weight a replicate with 3 events the same as one with 40). At SV_RATE_SCALE=1.0 GRCh38 yields ~25 translocations per run, so ~8 replicates is the design target. Per-replicate truth counts vary widely — BND 50–76, INV 9–20, DUP 18–34 — so single-replicate figures were never going to be stable.

⚠ Every recall below is PASS-only. truvari is run with --passonly, so a truth event whose only matching call was non-PASS is scored FN. Measured instance: chr10:45850637 was called by Manta with identical POS/END/SVLEN and scored a DUP false negative because its FILTER was MinSomaticScore. These recalls are therefore underestimates by an unmeasured amount — see #541.

Array 21072620 (2026-08-13, seeds 9–12, GRCh38 30x, purity 0.6, SV_RATE_SCALE=1.0), pooled over 4 replicates:

The INS rows ran in 1 and 2 of 4 replicates respectively and are not comparable with the rest; the aggregator says so itself rather than pooling them silently.

DUP: both callers miss the identical set, and one of them is ours

manta_DUP and delly_DUP agree to the record — TP=91, FN=14 — differing only in FP. Two independent callers missing the same events points at a property of the events, not the callers. Examined for one replicate (FN = 4):

So it is not a size floor (1.7 Mb is trivially detectable by depth) and not assembly gaps (three are 0.0% N; SV placement already requires a ±200 bp N-free window at both breakpoints, sv_model.rs ALIGNABLE_WINDOW, #224). One of the four is a filter artifact (#541). The remaining three — including a megabase-scale duplication with nothing called anywhere near it — are unexplained, and diagnosing them needs reads, so PRUNE_BAM=0 on a future run.

Calibration controls passed: truth-vs-truth recall 1.000 with all 71 scoreable records scored, decoy (shifted truth) recall 0.000; same for BND across all 50. So the matching configuration is not loose enough to accept anything.

BND geometry is fixed, confirmed at genome scale

This is the first whole-genome confirmation of the v3.1.0 fix. BND recall=0.000 was the symptom of a truth VCF describing a rearrangement the reads did not carry; it now reads:

BND geometry: 50 record(s), 50 parsed into a form, 0 unparsable
  t[p[  direct/deletion-like    9      t]p]  head-to-head            18
  [p[t  tail-to-tail           14      ]p]t  direct/duplication-like  9
  reciprocity: 0 unpaired, 0 mispaired

All four bracket forms present, every record parsed, 0 unpaired and 0 mispaired — the #451-era failure (truth emitted unmatchable by construction) is gone, and Manta independently recovers 46 of the 50 records from the reads alone.

INV precision 0.500 is ours, not the callers’ — and Gate 2 now shows the reads are clean

Both callers report TP=9, FP=9 — identically. The pipeline predicts this in its own pre-flight: 16 of the truth’s junctions are inversion-oriented (t]p] or [p[t), and a caller’s breakends for those convert to <INV> via Manta’s convertInversion.py, landing as INV false positives. Two independent callers arriving at the same number is the evidence that it is a representation artifact. INV recall 1.000 is real; INV precision from this tier is not quotable.

Gate 2 (eidolon/tests/gate2_realigned_inv.rs) closes the remaining inference. Realigning a homozygous 1.2 kb inversion with bwa-mem2 gives 123 breakpoint-clipped reads against 0 in the control, 98 same-orientation (FF/RR) pairs against 0, 0 everted (RF) pairs — so nothing suggests a duplication — and interior depth of 1.026x the flank, i.e. balanced. The reads carry an unambiguous inversion signature, so the false positives are not attributable to them. That was previously an argument from two callers agreeing; it is now a measurement.

INS: one planted, one verified present, zero found

truth INS: 1 (65 bp) was the #516 cap working as designed, and the drop log proved it rather than implying it: 2 insertions dropped (chr2, chr19) + 1 kept = 3 draws, against ~3.5 expected from Ins = 0.0279 over ~130 SVs, of which ~0.9 were expected to survive the 150 bp cap.

That cap no longer exists (#516 fixed, 2026-08-24). Every draw is now realized, so a run’s INS count is the model’s own rate rather than the rate times a survival fraction — which is why the numbers in this section cannot be compared to a post-fix run.

INS read support: 1 of 1 — 7 reads carry a 30-mer from the middle of the inserted sequence. This is the first campaign to check that the reads contained what the truth declared before the BAMs were pruned, which is exactly how #516 survived three campaigns.

⚠ But the check itself is not trustworthy yet. Job 21076622, same array, reports INS read support: 0 of 1 — which turned out to measure nothing. Its probe file is well formed (chr3 92124128 56bp + a valid 30-mer), yet no per-insertion line was emitted on stdout or stderr: count_probe_hits produced no output, the reporting loop ran zero iterations, and the 0 is an untouched initialiser printed against an independently counted n_probes. Worse, unsupported stayed 0, so the failure gate never fired and the replicate archived as a clean result (#540).

The 1 of 1 above came through the same code path. It did print its detail line, so it is probably a true pass, but treated the #516 cap as unverified in production until #540 is fixed and both replicates are re-checked. Moot for the cap itself, which is gone; the #540 concern about a replicate archiving a clean result without printing its detail line stands.

Neither caller found it. Manta’s only INS call was a false positive on another chromosome (chr11:99348945, 58 bp, non-PASS); nothing was called within 2 kb of the planted event. Delly emitted no INS records at all. So the miss is verified rather than inferred — the reads demonstrably carried the insertion. With n=1 this does not constrain a recall rate (95% CI ≈ [0, 0.98]); why it was missed is unresolved, and separating “somatic-score threshold at ~7 supporting reads” from “the reads lack a signature the aligner can present” needs PRUNE_BAM=0 and a look at the CIGARs.

A 30× GRCh38 replicate yields ~3 INS draws and ~1 after the cap. INS cannot be validated at caller level from single replicates at any runtime — it needs pooled replicates or an INS-enriched model. Across the ~8 replicates the aggregator targets that is still only ~8 planted insertions, so an INS-enriched model is likely the only route to a usable denominator.

What this tier does NOT establish

• Input fidelity at scale. Every cell above this section is still H1N1-only.
• Dosage. Detection is not dosage; the ~8% over-delivery (#499) was measured on a 1.2 kb event and has never been checked at these sizes.
• Subclonal behaviour. Resolved for de-novo SVs in #537. Supported DEL/DUP/INV/INS/CNV records now carry EIDOLON_CCF/EIDOLON_VAF and use CCF-adjusted depth and junction evidence; reciprocal translocation BND mates share one CCF. Delta yeast validation exercised CCF 0.2/0.5/1.0 and retained 199 chimeric alignments. This section’s older nominal-rate campaign predates that fix; the full-GRCh38 post-fix run above is the current release-candidate stress test.
• Small events. The planted set skews very large (DUP 3.4 Mb, DEL 9.1 Mb, INV 2.0 Mb). Multi-megabase events are easy to detect by depth; these recall figures say little about the 1–50 kb range where callers actually struggle.

De novo generation

BND de novo means inter-chromosomal translocation only. PCAWG’s TRA class is 100% inter-chromosomal (68,547/68,547, docs/pcawg_sv_measurement.md M1), and a same-contig junction is a DEL/DUP/INV by its orientation — already sampled at their own rates.

De novo INS is a LITERAL record, and it must carry SVTYPE=INS. Unlike every other sampled type, an insertion is emitted as REF=<anchor>, ALT=<anchor><novel bases> rather than as symbolic <INS> — that routes the novel sequence through gen-reads’ literal-insertion machinery, so reads spanning the locus actually carry the bases, and it matches what Manta emits for a resolved insertion. END equals POS (an insertion consumes no reference span) and SVLEN is positive and equal to ALT.len() - REF.len().

The tag is not cosmetic. Without it the record is indistinguishable from a small-indel-model insertion of the same size, so nothing downstream can attribute it — a length threshold cannot separate them, because the Ins length fit puts ~17% of structural insertions below 100 bp while the indel tail reaches ~93 bp. Measured consequence: SV validation jobs 20885875 and 20892075 planted 4 insertions between them (296, 1057, 1492, 2500 bp) and scored none, both reporting truth INS: 0, because sv_pipeline.sbatch selects SVs on ALT[0]~"<" || INFO/SVTYPE!=".". Every SV validation run before this had the same hole. Pinned by sampled_ins_carry_svtype_and_svlen_so_they_are_separable_from_indels and, across the sampler/writer boundary, sampled_ins_svtype_survives_the_vcf_writer.

De novo SVs are opt-in: sv_rate_scale defaults to 0.0. For human germline they are not recommended at all — constitutional balanced translocations occur in roughly 1 in 500 live births, so a normal genome carries ~zero. The machinery stays for cancer (rates are real and verified) and is expected to matter for plants.

Where input and de novo deliberately differ

These are not inconsistencies to reconcile. Input is “reproduce what I gave you”; de novo is “sample what the corpus says”. A researcher wanting a specific rearrangement should use input_vcf — that path preserves IDs, POS and ALT verbatim and generates junction reads.

How this was measured

# per cell: one variant, input_vcf, sv_rate_scale 0.0, produce_bam
eidolon gen-reads -c cell.yml
samtools index o.bam
samtools depth -a -r CONTIG:START-END o.bam | awk '{s+=$3;n++} END{print s/n}'   # vs a control span
zcat o_r1.fastq.gz | paste - - - - | awk '$1 ~ /EIDOLON_chimeric/' | grep -c PROBE

One methodological trap, hit while producing this table. Do not grep a FASTQ for a sequence probe without separating the sequence from the quality string: paste - - - - puts them on one line, and Phred+33 quality characters overlap the DNA alphabet (G is Q38, so a run of 12 Gs is common in quality). A homopolymer probe reported 71 false matches that way. Filter to the sequence field, and prefer a non-homopolymer probe. This is the same class of error as the FASTA reader accepting a FASTQ because > is Q29.

Current Benchmarks

We ran some benchmarks on eidolon on my home desktop, a very basic Linux desktop, using data pulled from public resources on the internet (mostly ncbi). The “yeast” is brewer’s yeast. These are the results:

[16:26:39] Build complete.
eidolon Benchmark Report
Run:        20260510_162639
Machine:    pop-os  |  8 logical CPUs
Coverage:   10x    |  Read length: 151 bp
Binary:     /home/joshfactorial/code/eidolon/target/release/eidolon
──────────────────────────────────────────────────────────────────────
SECTION 1: Single-ended read generation
Genome          Size(MB)   Wall time  Peak RSS(MB)      CPU%
──────────────────────────────────────────────────────────────────────
[16:26:40] Single-ended: ecoli  (4.6 MB)
ecoli                4.6     0:14.83         416.9       99%
[16:26:54]   Done: wall=0:14.83  peak_rss=416.9 MB  cpu=99%
[16:26:54] Single-ended: pneumonia  (2.1 MB)
pneumonia            2.1     0:05.39         197.0       99%
[16:27:00]   Done: wall=0:05.39  peak_rss=197.0 MB  cpu=99%
[16:27:00] Single-ended: yeast  (12 MB)
yeast                 12     0:10.94         981.7      352%
[16:27:11]   Done: wall=0:10.94  peak_rss=981.7 MB  cpu=352%
[16:27:11] Single-ended: c_elegans  (97 MB)
c_elegans             97     7:35.23        8591.2      387%
[16:34:46]   Done: wall=7:35.23  peak_rss=8591.2 MB  cpu=387%
──────────────────────────────────────────────────────────────────────
Section 1 — Passed: 4 / 4  |  Failed: 0
──────────────────────────────────────────────────────────────────────
SECTION 2: Paired vs single comparison  (genome: yeast, 10x)
  Paired-end uses fragment_mean=400, fragment_st_dev=50
  Same reference, same coverage, same seed — only the read-generation mode differs.
Mode               Wall time  Peak RSS(MB)      CPU%
──────────────────────────────────────────────────────────────────────
[16:34:46] Paired comparison — single-ended pass
single-ended         0:10.48         970.4      365%
[16:34:57]   Single done: wall=0:10.48  peak_rss=970.4 MB  cpu=365%
[16:34:57] Paired comparison — paired-ended pass
paired-ended         0:13.30        1126.8      354%
[16:35:10]   Paired done: wall=0:13.30  peak_rss=1126.8 MB  cpu=354%
──────────────────────────────────────────────────────────────────────
  Paired/single wall-time ratio: 1.27x  (10.48s → 13.30s)
  Paired/single peak-RSS ratio:  1.16x  (970.4 MB → 1126.8 MB)
──────────────────────────────────────────────────────────────────────

Upgrading from 3.x

v4.0.0 has three changes that can stop a run which completed under 3.x. Each fails before any output is written and says what to change.

Unknown config keys

Every top-level key must be one the subcommand reads; so must the keys inside adapters: (gen-reads) and frag_length: / gc_bias: (gen-bam-models). A key whose value is . is checked too.

A rejected key never did anything. In 3.x it was ignored, so the run used that setting’s default without warning: tumor_mutation_model:, meant as tumor_model:, simulated the tumor with the germline model. If a 3.x config now fails, the 3.x results it produced were made without that setting. Re-run anything that matters with the corrected key.

For a NEAT config, gen-reads names the eidolon equivalent where one exists (include_vcf → input_vcf, error_model → sequence_error_model, …) and otherwise says the option is not supported. Keys are not translated automatically.

gen-mut-model and REF mismatches

Each SNP’s and indel’s REF must match the reference at its position. The fit stops at the first record that does not, naming it and both bases. If the VCF does belong to this reference and only a few records disagree, drop them explicitly:

bcftools norm -f reference.fa --check-ref x in.vcf.gz -Oz -o checked.vcf.gz

If many records disagree, the VCF was most likely called against a different build. Fit against the matching reference instead.

Two related changes alter a model fitted from the same VCF:

• A SNP at a contig edge (no base on one side) is left out of the mutation rate. 3.x counted it.
• A bed_file that names no contig in the reference stops with both lists of contig names, usually a chr1 versus 1 mismatch, instead of “Unknown error”.

Output that changes without any action

These are fixes, not compatibility breaks, but they change output for the inputs they affect. Re-baseline anything measured on such inputs under 3.x:

• Soft-masked (lowercase) references now simulate exactly like uppercase ones.
• input_vcf runs that combine a literal deletion of 50 bp or more with a <DUP>, <DEL> or <CNV> no longer over-produce reads past the deletion.
• gen-cancer-reads with the bundled COSMIC models places somatic SNVs slightly differently; the models were rebuilt.
• gen-seq-error-model with max_reads samples the whole FASTQ, so a capped fit gives a different model and takes longer.
• filter-reads output is always gzip-compressed, so a plain x.fastq produces x_filter.fastq.gz.
• Golden BAM records gain an NM tag.

See CHANGELOG.md for the full list.

Upgrading from 2.0.0

v3.0.0 renamed the tokens eidolon emits. There is no behavior change — the same reads and variants are produced — but anything that parses eidolon output needs attention. In-tree consumers were all migrated; your own scripts were not.

VCFs — convert them

tools/migrate_legacy_tokens.sh old_truth.vcf.gz new_truth.vcf.gz
tools/migrate_legacy_tokens.sh --check old_truth.vcf.gz   # report only; exit 10 = legacy

Idempotent: safe to run on a file that’s already current, or on a non-eidolon VCF, so you can call it unconditionally in a pipeline. It rewrites the header declarations and the record fields together, so the result is a valid VCF.

Note you cannot work around this with a dual-name bcftools filter — bcftools rejects an undefined tag in a -i expression at parse time, so even INFO/EIDOLON_ORIGIN="somatic" || INFO/NEAT_ORIGIN="somatic" fails outright on a file that carries only one of the two. Convert the file instead.

Read names — NOT converted, and this is the part that can bite you

Read names are not rewritten, in either FASTQ or BAM. Rewriting files that are routinely hundreds of GB to change ~19 bytes per record isn’t a reasonable ask, so instead eidolon filter-reads accepts both prefixes natively — legacy FASTQs keep working with no conversion step, and it warns once when it sees an old prefix.

Nothing protects your own scripts. A grep '^@RNEAT_generated_', awk field split, or regex over read names will match zero records and exit 0 against v3.0.0 output — a silent empty result, not an error. Audit for the old prefix before upgrading:

grep -rn 'RNEAT_generated_\|RNEAT_chimeric_' your_scripts/

Best fix: make your parser prefix-agnostic

Only the prefix changed. Everything after _generated_ — <contig>_<start>_<end>_<uniq>/<mate> — is byte-for-byte identical between versions (verified across every read of a golden BAM re-prefixed both ways). So rather than rewriting files, strip the prefix and your script works against either version, and against any future rename:

# version-proof: drop whatever prefix is present, keep the encoded fields
samtools view x.bam | cut -f1 | sed -E 's/^[A-Z]+_generated_//'
zcat x_r1.fastq.gz  | awk 'NR%4==1' | sed -E 's/^@[A-Z]+_generated_//'

Same trick for chimeric reads: sed -E 's/^@?[A-Z]+_chimeric_//'.

If you are joining an old artifact against a new one on read name, normalize both sides this way. One unrelated gotcha for QNAME joins: eidolon’s golden BAM keeps the /1 /2 mate suffix, while bwa strips it from aligned BAMs — so strip that too (sed 's|/[12]$||') when joining a golden BAM to an aligner’s output. That mismatch predates this release.

If you must rewrite the files instead

# FASTQ (expensive — a full re-compress; prefer the prefix-agnostic parse above)
zcat old_r1.fastq.gz | sed 's/^@RNEAT_generated_/@EIDOLON_generated_/' | gzip > new_r1.fastq.gz

# BAM QNAMEs (QNAME is field 1, and no header line contains the prefix, so ^ is safe)
samtools view -h old.bam | sed 's/^RNEAT_generated_/EIDOLON_generated_/' | samtools view -b -o new.bam

BAM QNAMEs are covered by no in-tree tooling — filter-reads reads FASTQ, and the converter is bcftools-based so it cannot touch a BAM. Nothing eidolon ships parses BAM read names either, so this only affects your own scripts.

Project documents not in this site

docs/ also holds working documents that are not part of the user guide — design scopes, campaign plans, and the engineering audit. They are listed here so nothing in the repository is unreachable from the site, but they are not rendered as chapters: several are point-in-time records rather than current instructions, and folding them into the nav would present them as guidance.

Model provenance — what each shipped default is, where it came from, and whether it has been measured — lives in eidolon-core/src/models/model_data/README.md.