Skip to content

Commit ca769db

Browse files
authored
feat: Add random seed support for deterministic read simulation (#33)
Closes #30 Implements random seed support for both Illumina and ONT read simulators, enabling reproducible simulation results for scientific research and benchmarking. ## Features Delivered - Random seed parameter for Illumina pipeline (w-Wessim2 fragment simulation) - Random seed parameter for ONT pipeline (NanoSim native --seed support) - CLI options: --seed for both `reads illumina` and `reads ont` commands - Configuration file support in read_simulation.seed and nanosim_params.seed - Comprehensive logging of seed values for audit trails ## Quality Metrics - 582 tests pass (added 8 new tests) - 78% overall code coverage maintained - click_main.py: 51% → 90% coverage improvement - Zero regressions - All linting passes: ruff, mypy, bandit ## Implementation Approach - Followed DRY, KISS, SOLID principles - Backward compatible: seed=None maintains random behavior - Optional parameters throughout the stack - Comprehensive test coverage including determinism verification See IMPLEMENTATION_COMPLETE_ISSUE_30.md for full details.
2 parents fe17687 + 3557bb9 commit ca769db

15 files changed

Lines changed: 3077 additions & 2 deletions

CLAUDE.md

Lines changed: 250 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,250 @@
1+
# CLAUDE.md
2+
3+
This file provides guidance to Claude Code (claude.ai/code) when working with code in this repository.
4+
5+
## Overview
6+
7+
**MucOneUp** is a Python tool for simulating MUC1 VNTR (Variable Number Tandem Repeat) diploid references with customizable mutations and read simulation. The tool generates realistic genomic sequences with configurable repeat structures, applies targeted mutations, and can simulate both Illumina and Oxford Nanopore sequencing reads.
8+
9+
## Common Commands
10+
11+
### Installation
12+
```bash
13+
pip install .
14+
```
15+
16+
### Running Tests
17+
```bash
18+
# Run all tests
19+
python -m pytest tests/
20+
21+
# Run specific test file
22+
python -m pytest tests/test_simulate.py
23+
24+
# Run with verbose output (pytest)
25+
python -m pytest -v tests/
26+
27+
# Run with MucOneUp verbose flag
28+
muconeup --verbose --config config.json simulate --help
29+
muconeup -v --config config.json simulate --help # Short form
30+
```
31+
32+
### Basic Simulation
33+
```bash
34+
# Generate diploid haplotypes with random VNTR lengths
35+
muconeup --config config.json simulate --out-base output_name --out-dir output/
36+
37+
# Generate with fixed VNTR lengths
38+
muconeup --config config.json simulate --out-base output_name --fixed-lengths 60
39+
40+
# Generate series across length range (with progress bar)
41+
muconeup --config config.json simulate --out-base output_name --fixed-lengths 20-40 --simulate-series 5
42+
43+
# With verbose output
44+
muconeup --verbose --config config.json simulate --out-base output_name --fixed-lengths 20-40
45+
46+
# Apply mutations to specific positions
47+
muconeup --config config.json simulate --out-base output_name --mutation-name dupC --mutation-targets 1,25 2,30
48+
49+
# Dual simulation (normal + mutated)
50+
muconeup --config config.json simulate --out-base output_name --mutation-name normal,dupC
51+
```
52+
53+
### Progress Indicators
54+
When using `--simulate-series` with multiple iterations, a progress bar will appear:
55+
```
56+
Simulating 21 iterations [################-----------] 61% 00:02:15
57+
```
58+
This helps track long-running simulations.
59+
60+
### Read Simulation
61+
```bash
62+
# Simulate Illumina reads
63+
muconeup --config config.json --out-base output_name --simulate-reads
64+
65+
# With ORF prediction
66+
muconeup --config config.json --out-base output_name --output-orfs --orf-min-aa 100
67+
```
68+
69+
### Deterministic Read Simulation
70+
71+
Generate reproducible reads using the `--seed` parameter. This ensures identical read output across multiple runs, enabling reproducible research, fair benchmarking, and consistent debugging.
72+
73+
```bash
74+
# Illumina reads with seed (reproducible)
75+
muconeup --config config.json reads illumina sample.fa --seed 42
76+
77+
# ONT reads with seed (reproducible)
78+
muconeup --config config.json reads ont sample.fa --seed 42
79+
80+
# Full pipeline with seed (haplotypes + reads)
81+
muconeup --config config.json simulate --seed 42 --out-base sample
82+
muconeup --config config.json reads illumina sample.001.simulated.fa --seed 42
83+
84+
# Verify reproducibility
85+
muconeup --config config.json reads illumina test.fa --seed 42 --out-base run1
86+
muconeup --config config.json reads illumina test.fa --seed 42 --out-base run2
87+
diff run1_R1.fastq.gz run2_R1.fastq.gz # Should be identical
88+
```
89+
90+
**Important**: Same seed guarantees identical output ONLY when:
91+
- Using identical input files
92+
- Running on same platform/architecture
93+
- Using same tool versions (NanoSim, Python, reseq)
94+
95+
**Configuration Example**:
96+
```json
97+
{
98+
"read_simulation": {
99+
"simulator": "illumina",
100+
"seed": 42
101+
},
102+
"nanosim_params": {
103+
"training_data_path": "/path/to/model",
104+
"coverage": 30,
105+
"seed": 42
106+
}
107+
}
108+
```
109+
110+
### SNP Integration
111+
```bash
112+
# Generate random SNPs
113+
muconeup --config config.json --out-base output_name --random-snps --random-snp-density 1.0
114+
115+
# Apply SNPs from file
116+
muconeup --config config.json --out-base output_name --snp-file snps.tsv
117+
```
118+
119+
### Setting Up Conda Environments
120+
```bash
121+
# For Illumina read simulation (w-Wessim2)
122+
mamba env create -f conda/env_wessim.yaml
123+
124+
# For ONT read simulation (NanoSim)
125+
mamba env create -f conda/env_nanosim.yml
126+
```
127+
128+
## Architecture
129+
130+
### Core Simulation Pipeline
131+
132+
1. **Configuration Loading** (`config.py`): Validates JSON configuration against schema containing repeat definitions, probability transitions, mutation definitions, and tool paths.
133+
134+
2. **VNTR Generation** (`simulate.py`):
135+
- Builds haplotype chains by sampling repeats according to probability distributions
136+
- Forces canonical terminal block (6/6p → 7 → 8 → 9)
137+
- Supports both random length sampling and fixed-length generation
138+
- Can generate from predefined structure files
139+
140+
3. **Mutation Application** (`mutate.py`):
141+
- Applies insertions, deletions, replacements, or delete_insert operations
142+
- Supports strict mode to prevent auto-conversion of non-allowed repeats
143+
- Tracks mutated positions with "m" suffix in structure files
144+
- Records mutated VNTR unit sequences separately
145+
146+
4. **SNP Integration** (`snp_integrator.py`):
147+
- Parses TSV files with haplotype-specific SNPs (1-based haplotype, 0-based position)
148+
- Generates random SNPs with configurable density
149+
- Validates reference bases before applying (skippable in dual mutation mode)
150+
- Tracks successfully applied SNPs for reporting
151+
152+
5. **Sequence Assembly**: Concatenates left constant → repeat chain → right constant, supporting both hg19 and hg38 assemblies
153+
154+
### Read Simulation Pipelines
155+
156+
#### Illumina Pipeline (`read_simulator/pipeline.py`)
157+
Uses a port of w-Wessim2 with these steps:
158+
1. Replace Ns using reseq
159+
2. Generate systematic errors with reseq illuminaPE
160+
3. Convert to 2bit format with faToTwoBit
161+
4. Extract subset reference from sample BAM
162+
5. Align with pblat
163+
6. Simulate fragments (ported w-Wessim2 logic in `fragment_simulation.py`)
164+
7. Create paired reads with reseq seqToIllumina
165+
8. Split interleaved FASTQ
166+
9. Align to human reference with BWA MEM
167+
10. Optional coverage-based downsampling
168+
169+
#### ONT Pipeline (`read_simulator/ont_pipeline.py`)
170+
Uses NanoSim for Oxford Nanopore reads:
171+
1. Run NanoSim simulation with pre-trained models
172+
2. Align reads with minimap2
173+
3. Create indexed BAM output
174+
175+
### Key Modules
176+
177+
- **`cli.py`**: Argument parsing, simulation orchestration, dual simulation mode, series generation
178+
- **`distribution.py`**: Samples target VNTR length from normal/uniform distributions
179+
- **`probabilities.py`**: Weighted random selection for repeat transitions
180+
- **`fasta_writer.py`**: FASTA output with per-haplotype mutation annotations
181+
- **`io.py`**: Structure file parsing with embedded mutation information
182+
- **`translate.py`**: ORF prediction using orfipy
183+
- **`toxic_protein_detector.py`**: Scans ORFs for toxic features based on repeat structure and amino acid composition
184+
- **`simulation_statistics.py`**: Generates comprehensive JSON reports with runtime, haplotype metrics, and mutation details
185+
- **`analysis/vntr_statistics.py`**: Analyzes VNTR structures from CSV/TSV files, computes statistics (min/max/mean/median repeats), and builds transition probability matrices
186+
187+
### Configuration Structure
188+
189+
The `config.json` file contains:
190+
- **repeats**: Dictionary mapping repeat symbols (1, 2, X, A, B, etc.) to DNA sequences
191+
- **constants**: Left/right flanking sequences for hg19 and hg38, plus VNTR region coordinates
192+
- **probabilities**: State transition probabilities for repeat chain generation
193+
- **length_model**: Distribution parameters (normal/uniform) with min/max/mean/median repeats
194+
- **mutations**: Named mutation definitions with:
195+
- `allowed_repeats`: Valid repeat symbols for this mutation
196+
- `strict_mode`: Boolean to enforce allowed_repeats (prevents auto-conversion)
197+
- `changes`: List of operations (insert/delete/replace/delete_insert) with positions and sequences
198+
- **tools**: Command paths for external tools (reseq, bwa, samtools, etc.)
199+
- **read_simulation**: Parameters for Illumina pipeline (fragment size, coverage, threads, etc.)
200+
- **nanosim_params**: Parameters for ONT pipeline (training model path, coverage, read lengths)
201+
202+
### Mutation Strict Mode
203+
204+
By default, if a mutation target has a repeat not in `allowed_repeats`, it's auto-converted with a warning. Setting `"strict_mode": true` in a mutation definition causes an error instead. Random targets always respect `allowed_repeats` in both modes.
205+
206+
### Structure Files
207+
208+
Structure files can contain mutation information in header comments:
209+
```
210+
# Mutation Applied: dupC (Targets: [(1, 25)])
211+
haplotype_1 1-2-3-4-5-C-X-B-Xm-X-A-6-7-8-9
212+
haplotype_2 1-2-3-4-5-C-X-A-B-X-6p-7-8-9
213+
```
214+
215+
The "m" suffix marks mutated positions. Use `--input-structure` to generate from predefined chains.
216+
217+
### Dual Simulation Mode
218+
219+
When `--mutation-name normal,mutationName` is provided, two complete simulation runs occur:
220+
- `*.normal.fa` and `*.normal.simulation_stats.json`
221+
- `*.mut.fa` and `*.mut.simulation_stats.json`
222+
223+
In dual mode, SNP integration uses `skip_reference_check=True` for mutated sequences.
224+
225+
### Reference Assembly Support
226+
227+
The tool supports both hg19 and hg38 assemblies. Set `"reference_assembly": "hg38"` or `"hg19"` in config. Constants and VNTR regions are assembly-specific.
228+
229+
## Development Notes
230+
231+
- **Entry point**: `muconeup` command maps to `muc_one_up.cli:main` (defined in setup.cfg)
232+
- **Version management**: Single source in `muc_one_up/version.py`, imported by `__init__.py` and `cli.py`
233+
- **File naming**: Simulation outputs use numbered iterations (`.001`, `.002`, etc.) for series mode
234+
- **External tools**: All pipeline wrappers in `read_simulator/wrappers/` have timeout handling
235+
- **Intermediate files**: Pipeline uses underscore prefix (`_`) for temporary files
236+
- **Logging**: Configurable via `--log-level` (DEBUG/INFO/WARNING/ERROR/CRITICAL/NONE)
237+
- **Example data**: Located in `data/examples/`, includes `vntr_database.tsv` with real-world VNTR structures from published research
238+
239+
## Testing Strategy
240+
241+
Tests cover:
242+
- Configuration validation and loading (`test_config.py`)
243+
- Probability-based repeat selection (`test_probabilities.py`)
244+
- Haplotype simulation and chain building (`test_simulate.py`)
245+
- Mutation application and validation (`test_mutate.py`)
246+
- ORF prediction and translation (`test_translate.py`)
247+
- VNTR structure analysis (`test_vntr_statistics.py`)
248+
- CLI integration for vntr-stats command (`test_click_vntr_stats.py`)
249+
250+
When adding features, ensure corresponding tests validate both success and error cases.

0 commit comments

Comments
 (0)