diff --git a/.github/workflows/release.yml b/.github/workflows/release.yml index 8b62a94..3f8d118 100644 --- a/.github/workflows/release.yml +++ b/.github/workflows/release.yml @@ -108,6 +108,13 @@ jobs: steps: - uses: actions/download-artifact@v4 + - uses: actions/checkout@v4 + with: + path: source + + - name: Package example data + run: tar -czf madcats-demo.tar.gz --exclude=results -C source/examples madcats-demo + - name: Create Release uses: softprops/action-gh-release@v2 with: @@ -133,4 +140,5 @@ jobs: madcats-linux-x86_64/*.tar.gz madcats-macos-x86_64/*.tar.gz madcats-macos-aarch64/*.tar.gz + madcats-demo.tar.gz dist/* diff --git a/.gitignore b/.gitignore index 7c67c2a..9c780c5 100644 --- a/.gitignore +++ b/.gitignore @@ -184,3 +184,4 @@ docs/book/ # macOS .DS_Store +examples/madcats-demo/results/ diff --git a/README.md b/README.md index 90a700b..40d762a 100644 --- a/README.md +++ b/README.md @@ -1,4 +1,4 @@ -# MADCATS +# MADCATS MADCATS logo [![License: GPL v3](https://img.shields.io/badge/License-GPLv3-blue.svg)](https://www.gnu.org/licenses/gpl-3.0) [![PyPI](https://img.shields.io/pypi/v/madcats.svg)](https://pypi.org/project/madcats/) @@ -81,6 +81,8 @@ pip install . # or: python3 madcats.py --help (no install needed) ## Quick Start +To try MADCATS on a small simulated project first, see [Try the Example](https://madcats.readthedocs.io/en/latest/getting-started/example.html). + ```bash # 1) Reconstruct each marker's Ref/Alt tags from the DArTag design files madcats design \ @@ -131,7 +133,7 @@ Full documentation is available at [madcats.readthedocs.io](https://madcats.read - [Evidence Model](https://madcats.readthedocs.io/en/latest/reference/evidence-model.html): how alleles are counted, discovered and denoised The design rationale and version-by-version measurements behind every default are recorded in the -[development reference](docs/DEVELOPMENT.md). +[development record](https://github.com/Breeding-Insight/MADCATS/blob/main/docs/DEVELOPMENT.md). ## Commands diff --git a/assets/madcats-logo-full.png b/assets/madcats-logo-full.png new file mode 100644 index 0000000..41107b4 Binary files /dev/null and b/assets/madcats-logo-full.png differ diff --git a/docs/DEVELOPMENT.md b/docs/DEVELOPMENT.md index d06f957..9519c80 100644 --- a/docs/DEVELOPMENT.md +++ b/docs/DEVELOPMENT.md @@ -15,8 +15,7 @@ and VCF-ready output. `madcats.py` uses only Python 3's standard library, exposes all analysis thresholds on the command line, and writes a QC and provenance JSON beside every result. (The name also reflects the tool's origins as a FASTQ-to-MADC workflow, where MADC denotes the Missing -Allele Discovery Counts format; `madc_pipeline.py` remains as a -compatibility shim for existing scripts.) +Allele Discovery Counts format.) **A prior MADC is not an input to the workflow.** Generation uses the FASTQs plus the assay design: either the DArTag primer/probe files or a diff --git a/docs/src/README.md b/docs/src/README.md index 7771a9a..8417450 100644 --- a/docs/src/README.md +++ b/docs/src/README.md @@ -1,3 +1,5 @@ +MADCATS logo + # Welcome to MADCATS! [![License: GPL v3](https://img.shields.io/badge/License-GPLv3-blue.svg)](https://www.gnu.org/licenses/gpl-3.0) @@ -69,7 +71,8 @@ New to MADCATS? Start here: 1. [Installation](getting-started/installation.md): install with pip or download a binary 2. [Quick Start](getting-started/quickstart.md): your first MADC in two commands -3. [Choosing a Workflow](getting-started/workflows.md): which commands your platform and inputs need +3. [Try the Example](getting-started/example.md): run a small simulated project and check it against the truth +4. [Choosing a Workflow](getting-started/workflows.md): which commands your platform and inputs need ## Scope diff --git a/docs/src/SUMMARY.md b/docs/src/SUMMARY.md index aecdc6c..49674f2 100644 --- a/docs/src/SUMMARY.md +++ b/docs/src/SUMMARY.md @@ -6,6 +6,7 @@ - [Installation](getting-started/installation.md) - [Quick Start](getting-started/quickstart.md) +- [Try the Example](getting-started/example.md) - [Choosing a Workflow](getting-started/workflows.md) - [Input Formats](getting-started/input-formats.md) @@ -39,7 +40,6 @@ - [Architecture](development/architecture.md) - [Contributing](development/contributing.md) - [Changelog](development/changelog.md) -- [Development Reference](development/reference.md) ----------- diff --git a/docs/src/citation.md b/docs/src/citation.md index 1161bde..b06f7ea 100644 --- a/docs/src/citation.md +++ b/docs/src/citation.md @@ -47,5 +47,5 @@ polyploids from messy sequencing data. Genetics 210(3):789-807. In your methods section: > Microhaplotype read counts were generated from the demultiplexed FASTQ files with MADCATS -> v0.12.1 (Sandercock, 2026), using panel tags reconstructed from the vendor design files and +> v0.13.0 (Sandercock, 2026), using panel tags reconstructed from the vendor design files and > default evidence and denoising settings. diff --git a/docs/src/commands/generate.md b/docs/src/commands/generate.md index fa279ad..7aef3b1 100644 --- a/docs/src/commands/generate.md +++ b/docs/src/commands/generate.md @@ -1,6 +1,6 @@ # madcats generate -Generate a MADC-compatible microhaplotype count table from per-sample FASTQs and a panel. +Generate a MADC microhaplotype count table from per-sample FASTQs and a panel. ## Synopsis @@ -55,7 +55,7 @@ The [Evidence Model](../reference/evidence-model.md) explains stages 2–6 in de | `--sample-lookup ` | none | Two-column CSV/TSV mapping FASTQ-derived ids to sample names. Every run sample must be mapped | | `--tag-length ` | 81 | Tag length | -## Read Assignment (stage 2; these settings key the cache) +## Read Assignment (stage 2) | Option | Default | Description | |--------|---------|-------------| @@ -64,7 +64,7 @@ The [Evidence Model](../reference/evidence-model.md) explains stages 2–6 in de | `--primer-mismatches {0,1}` | 1 | Mismatches allowed while assigning a primer/anchor. With `0`, reads from Match haplotypes whose distinguishing SNP falls inside the anchor are lost | | `--gapped-rescue {on,off}` | `on` | Give reads that fail substitution-only identity one gapped alignment, so undesigned indels are counted | | `--max-indel-length ` | 10 | Widest gap (alignment band) a gapped alignment may use | -| `--partial-read-mode {off,assign}` | `assign` | Add a read shorter than the tag to the one full-length allele (already passing the evidence rules) whose start it matches, only when it is longer than every prefix shared between that locus's alleles. `off` drops such reads | +| `--partial-read-mode {on,off}` | `on` | Add a read shorter than the tag to the one full-length allele (already passing the evidence rules) whose start it matches, only when it is longer than every prefix shared between that locus's alleles. `off` drops such reads | | `--min-partial-length ` | tag length − 9 | Minimum length of a credited partial read | | `--cache-dir ` | none | Reusable per-sample count cache | | `--jobs ` | up to 4 | Samples processed in parallel (also parallelizes denoising) | @@ -85,8 +85,8 @@ The [Evidence Model](../reference/evidence-model.md) explains stages 2–6 in de | Option | Default | Description | |--------|---------|-------------| -| `--denoise-distance {0,1,2}` | 1 | Per-sample abundance-skew denoising distance | -| `--global-denoise-distance {0,1,2}` | 1 | Edit distance for the pooled same-allele merge; `0` disables | +| `--denoise-distance {0,1,2}` | 1 | Per-sample denoising: the most differing bases allowed between a rare sequence and its parent; `0` disables | +| `--global-denoise-distance {0,1,2}` | 1 | Most differing bases allowed in the pooled same-allele merge; `0` disables | | `--denoise-alpha ` | 2.0 | UNOISE-style abundance-skew exponent | ## Output (stages 6–7) diff --git a/docs/src/development/architecture.md b/docs/src/development/architecture.md index 030aca0..8fb2faf 100644 --- a/docs/src/development/architecture.md +++ b/docs/src/development/architecture.md @@ -9,10 +9,11 @@ madcats.py`), as the installed `madcats` command, or frozen into a standalone bi | Path | Contents | |------|----------| | `madcats.py` | The whole tool: library functions and the command line | -| `madc_pipeline.py` | Compatibility shim re-exporting `madcats` under the tool's former name | | `tests/` | Unit and end-to-end tests (`unittest`) | | `scripts/` | Diagnostic and conversion scripts (e.g. `madc_to_reference_tags.py`, `validate_inferred_design.py`) | | `docs/` | This documentation (mdBook), plus `DEVELOPMENT.md`, the long-form development record | +| `examples/madcats-demo/` | Simulated example project and the script that generates it | +| `assets/` | Full-size logo (2028 × 2028 PNG); the docs use a 600 px copy in `docs/src/images/` | | `packaging/` | PyInstaller entry point for the standalone binaries | | `pyproject.toml` | Package metadata for PyPI | | `.github/workflows/` | Tests on every push; PyPI and binary releases on version tags | diff --git a/docs/src/development/changelog.md b/docs/src/development/changelog.md index 9076f8c..3393bb2 100644 --- a/docs/src/development/changelog.md +++ b/docs/src/development/changelog.md @@ -1,17 +1,28 @@ # Changelog Notable changes by version. Measurements behind each change are in the -[Development Reference](reference.md). +[development record](https://github.com/Breeding-Insight/MADCATS/blob/main/docs/DEVELOPMENT.md) on GitHub. + +## 0.13.0 + +- Removed the `madc_pipeline` compatibility module. Import `madcats` instead; the bundled + scripts already do. +- `--partial-read-mode` takes `on`/`off` like the other switches; the old value `assign` still + works as `on`, with a warning. +- Clearer command-line help: every command has a description, every option has help text, and + the program is named `madcats` in usage lines. +- The development record (`docs/DEVELOPMENT.md`) is no longer published on the docs site; it + stays in the repository. +- A small simulated example project, with its true haplotypes and genotypes, is included in + `examples/madcats-demo` and attached to each release as `madcats-demo.tar.gz`. +- The test suite checks that the command reference pages match each command's `--help`. -## Unreleased +## 0.12.1 - Packaged for PyPI: `pip install madcats` installs a `madcats` command. - Standalone binaries for Linux and macOS are built for each release. -- New documentation site; the former README is now the development reference. +- New documentation site; the former README is now the development record, `docs/DEVELOPMENT.md`. - Licensed under GPL-3.0-or-later. - -## 0.12.1 - - `generate`'s finalisation (per-sample denoising and the pooled pass) runs across the `--jobs` workers. The selection is identical to the sequential passes. On a 44-run rhAmpSeq family the pooled pass alone had taken 1,123 s in one process. @@ -28,7 +39,7 @@ Notable changes by version. Measurements behind each change are in the ## 0.11.0 -- Partial reads (`--partial-read-mode assign`, default): reads shorter than the tag are +- Partial reads (`--partial-read-mode assign`, now `on`; default): reads shorter than the tag are added to an allele that already passed the evidence rules and whose start they match, only when this cannot change an allele ratio. `--min-partial-length` defaults to tag length − 9. The experimental `extend` mode was removed. @@ -97,5 +108,5 @@ Notable changes by version. Measurements behind each change are in the - `infer-design` for panels without a primer file. - 0.3.x–0.4.x: quality-gated discovery and the exploratory audits that motivated 0.5.0. -- The tool was renamed from `madc_pipeline` to MADCATS. `madc_pipeline.py` remains as a - compatibility shim. +- The tool was renamed from `madc_pipeline` to MADCATS, with `madc_pipeline.py` kept as a + compatibility shim until 0.13.0. diff --git a/docs/src/development/contributing.md b/docs/src/development/contributing.md index eb3792f..6ae3ea3 100644 --- a/docs/src/development/contributing.md +++ b/docs/src/development/contributing.md @@ -28,10 +28,18 @@ indel support, normalization and window freezing, output framing, paired-end mer window-scan, gapped alignment, denoising and partial reads. Tests run on every push and pull request across Python 3.9–3.13 (`.github/workflows/tests.yml`). +Two tests guard the documentation. `tests/test_docs.py` checks that each page in +`docs/src/commands/` lists exactly the options that command's `--help` shows, with matching +choices and defaults, so update the page whenever you add or change an option. +`tests/test_example.py` runs the example in `examples/madcats-demo` and checks the output against +its simulated truth. If a change alters the example's numbers, update +[Try the Example](../getting-started/example.md) to match. To change the example data itself, +edit and rerun `examples/madcats-demo/make_demo_data.py`. + Changes that touch counting should also be checked against real data. The standard is byte-for-byte: a refactor regenerates existing tables identically from cache, and a behaviour change is measured against a comparator before it becomes a default. Record those measurements in -`docs/DEVELOPMENT.md`. +[`docs/DEVELOPMENT.md`](https://github.com/Breeding-Insight/MADCATS/blob/main/docs/DEVELOPMENT.md). ## Code Style diff --git a/docs/src/development/reference.md b/docs/src/development/reference.md deleted file mode 100644 index fb8ad9e..0000000 --- a/docs/src/development/reference.md +++ /dev/null @@ -1 +0,0 @@ -{{#include ../../DEVELOPMENT.md}} diff --git a/docs/src/getting-started/example.md b/docs/src/getting-started/example.md new file mode 100644 index 0000000..596b903 --- /dev/null +++ b/docs/src/getting-started/example.md @@ -0,0 +1,110 @@ +# Try the Example + +MADCATS comes with a small simulated DArTag-style project that runs in seconds. The true +haplotypes and genotypes behind its reads are included, so you can check every row of the output +against what was simulated. + +## Get the Data + +Download it from the latest release (MADCATS 0.13.0 and later): + +```bash +curl -LO https://github.com/Breeding-Insight/MADCATS/releases/latest/download/madcats-demo.tar.gz +tar -xzf madcats-demo.tar.gz +cd madcats-demo +``` + +The same files are in the repository under `examples/madcats-demo`. + +| File | Contents | +|------|----------| +| `primers.csv` | DArT-style primer CSV for the 8 markers | +| `probe_design.tsv` | Probe design with one bracketed SNP per marker | +| `fastq/` | One FASTQ per sample, `S01.fastq.gz` … `S12.fastq.gz` | +| `sample_names.csv` | Renames the samples from `S01` to `Line_01`, and so on | +| `truth_haplotypes.csv` | Every simulated haplotype's 81-bp tag sequence | +| `truth_genotypes.csv` | Each sample's simulated diploid genotype at each marker | +| `make_demo_data.py` | The script that made all of the above; rerunning it reproduces the files exactly | + +## What Was Simulated + +Twelve diploid samples, eight markers with 81-bp tags, and about 550 reads per sample. Beyond each +marker's designed Ref and Alt alleles: + +- **Chr02** and **Chr04** carry an extra microhaplotype: a second SNP on the Ref or Alt background. +- **Chr06** carries a haplotype with an undesigned 1-bp insertion. +- **Chr08** has no Alt allele in any sample. +- **Chr03** and **Chr07** have their primer on the reverse strand of the probe-design target. +- Every read has a 0.2% chance per base of a sequencing error, about 8% of reads are shorter than + the tag, and each sample has a few off-target reads. + +## Run It + +Reconstruct the panel from the design files, then count: + +```bash +madcats design \ + --primers primers.csv \ + --probe-design probe_design.tsv \ + --output results/reference_tags.csv + +madcats generate \ + --fastq-dir fastq \ + --reference-tags results/reference_tags.csv \ + --sample-lookup sample_names.csv \ + --output results/demo_MADC.csv +``` + +`generate` writes the MADC table, `results/demo_MADC.csv.qc.json` and +`results/demo_MADC.csv.variants.csv`. + +## Check the Results + +The table has 19 rows: a Ref and an Alt row for each of the 8 markers, plus the 3 microhaplotypes +MADCATS discovered. With MADCATS 0.13.0 the discovered rows are: + +| AlleleID | Variants | MarkerType | Reads (all samples) | +|----------|----------|------------|--------------------:| +| `Chr02_001257730\|RefMatch` | `67:A>T` | `snp` | 150 | +| `Chr04_002318006\|AltMatch` | `25:G>T;57:A>T` | `snp` | 122 | +| `Chr06_001644492\|RefMatch` | `29:G>GC;clip=1` | `indel` | 128 | + +Things to look for: + +- **Every simulated haplotype is a row, and nothing else is.** Compare the `AlleleSequence` column + with `truth_haplotypes.csv`. The sequencing errors were merged back into their true haplotypes by + denoising, so none of them became rows. +- **The insertion was recovered.** Its reads failed the substitution-only identity check and were + placed by gapped alignment. `Variants` reports the insertion left-aligned (`29:G>GC`), and + `clip=1` notes that it pushed the last Ref base out of the 81-bp window. +- **Chr08 still has an Alt row.** Canonical Ref and Alt rows are always written. Its single read + comes from a sequencing error that happened to hit the SNP base. +- **Counts follow the genotypes.** In each `Line_` column, a heterozygous sample splits its reads + between two haplotypes and a homozygous sample puts nearly all of them on one. Compare with + `truth_genotypes.csv`. +- **Every read is accounted for.** The `read_accounting` section of the QC JSON shows 6,989 reads in + total, 129 assigned by gapped alignment, 466 partial reads added to counts, 851 merged away by + per-sample denoising, and 64 in error sequences that failed the evidence rules. Each sample's + entry in `sample_qc` also counts its off-target reads (`reads_no_primer`). + +## Try Changing Something + +Each of these reruns `generate` with one setting changed: + +| Change | What happens | +|--------|--------------| +| `--reproducible-min-samples 5 --depth-min-reads 0` | The 3 discovered rows disappear: each extra haplotype is carried by only 4 samples, and the depth route is off | +| `--gapped-rescue off` | The Chr06 insertion row disappears; its reads are no longer assigned | +| `--partial-read-mode off` | Same 19 rows, but reads in rows fall from 6,661 to 6,195 | + +To compare two runs, use `validate`: + +```bash +madcats validate --generated results/other_MADC.csv --expected results/demo_MADC.csv +``` + +## Next Steps + +- [Choosing a Workflow](workflows.md): which commands your own data needs +- [Evidence Model](../reference/evidence-model.md): how reads are assigned, discovered and denoised +- [Output Formats](../reference/output-formats.md): every column of the MADC table diff --git a/docs/src/getting-started/installation.md b/docs/src/getting-started/installation.md index a98d88b..ba0b267 100644 --- a/docs/src/getting-started/installation.md +++ b/docs/src/getting-started/installation.md @@ -111,7 +111,7 @@ madcats --version ``` ``` -madcats 0.12.1 +madcats 0.13.0 ``` To run the test suite from a source checkout: @@ -123,4 +123,5 @@ python3 -m unittest discover -s tests -v ## Next Steps - [Quick Start](quickstart.md): generate your first MADC +- [Try the Example](example.md): run a small simulated project in seconds - [Choosing a Workflow](workflows.md): which commands your data needs diff --git a/docs/src/getting-started/quickstart.md b/docs/src/getting-started/quickstart.md index 5f49898..893a9f7 100644 --- a/docs/src/getting-started/quickstart.md +++ b/docs/src/getting-started/quickstart.md @@ -89,6 +89,7 @@ See [Output Formats](../reference/output-formats.md) for details. ## Next Steps +- [Try the Example](example.md): run these steps on a small simulated project and check the output - [Choosing a Workflow](workflows.md): routes for other platforms and inputs - [madcats generate](../commands/generate.md): every option, grouped by pipeline stage - [Evidence Model](../reference/evidence-model.md): when a new sequence is reported as a row diff --git a/docs/src/images/madcats-logo.png b/docs/src/images/madcats-logo.png new file mode 100644 index 0000000..abd9af2 Binary files /dev/null and b/docs/src/images/madcats-logo.png differ diff --git a/docs/src/reference/evidence-model.md b/docs/src/reference/evidence-model.md index 8b252b3..c5c8eed 100644 --- a/docs/src/reference/evidence-model.md +++ b/docs/src/reference/evidence-model.md @@ -2,7 +2,7 @@ This page explains how `madcats generate` decides which reads count, which sequences become rows, and how sequencing errors are removed. The design rationale and the measurements behind each default -are in the [Development Reference](../development/reference.md). +are in the [development record](https://github.com/Breeding-Insight/MADCATS/blob/main/docs/DEVELOPMENT.md) on GitHub. ## Read Assignment @@ -61,8 +61,8 @@ depth route. Canonical sequences and sequences that pass the evidence rules are protected: denoising never merges them into anything else. -1. Per-sample denoising (UNOISE-style abundance skew). Within each sample, a rare sequence within - `--denoise-distance` (1) edits of a more abundant parent is merged into the parent when the +1. Per-sample denoising (UNOISE-style abundance skew). Within each sample, a rare sequence that differs at + no more than `--denoise-distance` (1) bases from a more abundant parent is merged into the parent when the parent is at least `2^(alpha·d + 1)` times more abundant (`--denoise-alpha` 2.0). The closest parent is chosen first, then the most abundant. 2. Pooled same-allele merge (`--global-denoise-distance`, 1). A sequence that failed the evidence diff --git a/docs/src/reference/validation.md b/docs/src/reference/validation.md index 92a85f1..f69b891 100644 --- a/docs/src/reference/validation.md +++ b/docs/src/reference/validation.md @@ -2,7 +2,7 @@ MADCATS has been validated on five platforms, each against the best independent comparison available for that dataset. This page summarizes the results. The full record, including every -version-to-version comparison, is in the [Development Reference](../development/reference.md). +version-to-version comparison, is in the [development record](https://github.com/Breeding-Insight/MADCATS/blob/main/docs/DEVELOPMENT.md) on GitHub. ## Performance at a Glance diff --git a/docs/theme/favicon.png b/docs/theme/favicon.png new file mode 100644 index 0000000..118a342 Binary files /dev/null and b/docs/theme/favicon.png differ diff --git a/examples/madcats-demo/README.md b/examples/madcats-demo/README.md new file mode 100644 index 0000000..14b0b0e --- /dev/null +++ b/examples/madcats-demo/README.md @@ -0,0 +1,17 @@ +# MADCATS example: a simulated DArTag-style project + +Eight markers, twelve diploid samples, and the true haplotypes and genotypes behind every read. +The data are simulated by `make_demo_data.py`; rerunning it reproduces these files exactly. + +```bash +madcats design --primers primers.csv --probe-design probe_design.tsv \ + --output results/reference_tags.csv + +madcats generate --fastq-dir fastq --reference-tags results/reference_tags.csv \ + --sample-lookup sample_names.csv --output results/demo_MADC.csv +``` + +Then compare `results/demo_MADC.csv` with `truth_haplotypes.csv` and `truth_genotypes.csv`. + +The walkthrough, including what was simulated and what to look for in the output, is at +https://madcats.readthedocs.io/en/latest/getting-started/example.html diff --git a/examples/madcats-demo/fastq/S01.fastq.gz b/examples/madcats-demo/fastq/S01.fastq.gz new file mode 100644 index 0000000..2f18872 Binary files /dev/null and b/examples/madcats-demo/fastq/S01.fastq.gz differ diff --git a/examples/madcats-demo/fastq/S02.fastq.gz b/examples/madcats-demo/fastq/S02.fastq.gz new file mode 100644 index 0000000..d0eea5a Binary files /dev/null and b/examples/madcats-demo/fastq/S02.fastq.gz differ diff --git a/examples/madcats-demo/fastq/S03.fastq.gz b/examples/madcats-demo/fastq/S03.fastq.gz new file mode 100644 index 0000000..bacdddd Binary files /dev/null and b/examples/madcats-demo/fastq/S03.fastq.gz differ diff --git a/examples/madcats-demo/fastq/S04.fastq.gz b/examples/madcats-demo/fastq/S04.fastq.gz new file mode 100644 index 0000000..9ceb1ad Binary files /dev/null and b/examples/madcats-demo/fastq/S04.fastq.gz differ diff --git a/examples/madcats-demo/fastq/S05.fastq.gz b/examples/madcats-demo/fastq/S05.fastq.gz new file mode 100644 index 0000000..3597cca Binary files /dev/null and b/examples/madcats-demo/fastq/S05.fastq.gz differ diff --git a/examples/madcats-demo/fastq/S06.fastq.gz b/examples/madcats-demo/fastq/S06.fastq.gz new file mode 100644 index 0000000..e98eacf Binary files /dev/null and b/examples/madcats-demo/fastq/S06.fastq.gz differ diff --git a/examples/madcats-demo/fastq/S07.fastq.gz b/examples/madcats-demo/fastq/S07.fastq.gz new file mode 100644 index 0000000..4f4f0bd Binary files /dev/null and b/examples/madcats-demo/fastq/S07.fastq.gz differ diff --git a/examples/madcats-demo/fastq/S08.fastq.gz b/examples/madcats-demo/fastq/S08.fastq.gz new file mode 100644 index 0000000..0706369 Binary files /dev/null and b/examples/madcats-demo/fastq/S08.fastq.gz differ diff --git a/examples/madcats-demo/fastq/S09.fastq.gz b/examples/madcats-demo/fastq/S09.fastq.gz new file mode 100644 index 0000000..db37114 Binary files /dev/null and b/examples/madcats-demo/fastq/S09.fastq.gz differ diff --git a/examples/madcats-demo/fastq/S10.fastq.gz b/examples/madcats-demo/fastq/S10.fastq.gz new file mode 100644 index 0000000..4437c7c Binary files /dev/null and b/examples/madcats-demo/fastq/S10.fastq.gz differ diff --git a/examples/madcats-demo/fastq/S11.fastq.gz b/examples/madcats-demo/fastq/S11.fastq.gz new file mode 100644 index 0000000..7b427af Binary files /dev/null and b/examples/madcats-demo/fastq/S11.fastq.gz differ diff --git a/examples/madcats-demo/fastq/S12.fastq.gz b/examples/madcats-demo/fastq/S12.fastq.gz new file mode 100644 index 0000000..3377174 Binary files /dev/null and b/examples/madcats-demo/fastq/S12.fastq.gz differ diff --git a/examples/madcats-demo/make_demo_data.py b/examples/madcats-demo/make_demo_data.py new file mode 100644 index 0000000..7e5d530 --- /dev/null +++ b/examples/madcats-demo/make_demo_data.py @@ -0,0 +1,207 @@ +"""Write the MADCATS example dataset: a small simulated DArTag-style project. + +Everything in this folder except README.md is produced by this script, so the +data can be regenerated exactly and every planted feature is visible here: + + python make_demo_data.py + +The panel has 8 markers with 81-bp tags and the project has 12 diploid +samples. On top of each marker's designed Ref and Alt alleles the reads carry: + +- Chr02 and Chr04: an additional microhaplotype (a second SNP on the Ref or + Alt background) that `generate` should discover; +- Chr06: an undesigned 1-bp insertion haplotype, recovered by gapped alignment; +- Chr08: no sample carries the Alt allele, so its Alt row is still written but + holds only a stray sequencing-error read; +- Chr03 and Chr07: primers on the reverse strand of the probe-design target; +- sequencing errors (0.2% per base), reads shorter than the tag, and + off-target reads, so denoising, partial-read crediting and the QC read + accounting all have something to do. + +truth_haplotypes.csv and truth_genotypes.csv record what was simulated. +""" +import csv +import gzip +import random +from pathlib import Path + +HERE = Path(__file__).resolve().parent +SEED = 2026 +TAG_LENGTH = 81 +PRIMER_LENGTH = 20 +GENOME_LENGTH = 140 +PRIMER_OFFSET = 10 # bases of target sequence before the left primer +AMPLICON_LENGTH = 110 # bases sequenced from the start of the primer +N_SAMPLES = 12 +READS_PER_COPY = (25, 45) # reads per haplotype copy per sample +ERROR_RATE = 0.002 +PARTIAL_FRACTION = 0.08 +OFF_TARGET_READS = (10, 20) + +COMPLEMENT = str.maketrans("ACGT", "TGCA") + +# (marker, primer on the reverse strand of the target, variant position in the tag) +MARKERS = [ + ("Chr01_000412883", False, 34), + ("Chr02_001257730", False, 48), + ("Chr03_000086214", True, 41), + ("Chr04_002318006", False, 57), + ("Chr05_000970341", False, 29), + ("Chr06_001644492", False, 52), + ("Chr07_000203517", True, 63), + ("Chr08_001105968", False, 45), +] + + +def reverse_complement(sequence): + return sequence.translate(COMPLEMENT)[::-1] + + +def other_base(rng, base): + return rng.choice([candidate for candidate in "ACGT" if candidate != base]) + + +def substitute(sequence, position, base): + return sequence[:position] + base + sequence[position + 1:] + + +def build_marker(rng, index, name, reverse, variant_position): + """Return the design rows and the simulated haplotypes (full amplicons) of one marker.""" + genome = "".join(rng.choice("ACGT") for _ in range(GENOME_LENGTH)) + snp = PRIMER_OFFSET + variant_position - 1 + ref_base = genome[snp] + alt_base = other_base(rng, ref_base) + left_primer = genome[PRIMER_OFFSET:PRIMER_OFFSET + PRIMER_LENGTH] + right_primer = reverse_complement(genome[-PRIMER_OFFSET - PRIMER_LENGTH:-PRIMER_OFFSET]) + + if reverse: + target = reverse_complement(genome) + bracket = GENOME_LENGTH - 1 - snp + bracketed = f"[{reverse_complement(ref_base)}/{reverse_complement(alt_base)}]" + else: + target = genome + bracket = snp + bracketed = f"[{ref_base}/{alt_base}]" + target_sequence = target[:bracket] + bracketed + target[bracket + 1:] + + amplicon = genome[PRIMER_OFFSET:PRIMER_OFFSET + AMPLICON_LENGTH] + tag_snp = variant_position - 1 + haplotypes = { + "Ref": (amplicon, "designed Ref allele"), + "Alt": (substitute(amplicon, tag_snp, alt_base), "designed Alt allele"), + } + if index == 1: + second = 66 + haplotypes["RefMatch"] = ( + substitute(amplicon, second, other_base(rng, amplicon[second])), + f"Ref allele plus a second SNP at tag position {second + 1}", + ) + if index == 3: + second = 24 + alt = haplotypes["Alt"][0] + haplotypes["AltMatch"] = ( + substitute(alt, second, other_base(rng, alt[second])), + f"Alt allele plus a second SNP at tag position {second + 1}", + ) + if index == 5: + insertion = 30 + haplotypes["Insertion"] = ( + amplicon[:insertion] + other_base(rng, amplicon[insertion]) + amplicon[insertion:-1], + f"Ref allele with an undesigned 1-bp insertion near tag position {insertion}", + ) + + primer_row = [name, left_primer, right_primer, ref_base, alt_base, variant_position] + design_row = [name, target_sequence] + return primer_row, design_row, haplotypes + + +def genotypes_for(index, haplotypes): + """Fixed diploid genotypes: every allele is carried by at least three samples.""" + dosage_cycle = [0, 1, 2, 1, 0, 1, 2, 2, 1, 0, 1, 1] + calls = [] + for sample in range(N_SAMPLES): + alt_copies = dosage_cycle[(sample + index) % N_SAMPLES] + if index == 7: + alt_copies = 0 + pair = ["Alt"] * alt_copies + ["Ref"] * (2 - alt_copies) + extra = [name for name in haplotypes if name not in ("Ref", "Alt")] + if extra and sample % 3 == 0: + pair[sample % 2] = extra[0] + calls.append(sorted(pair)) + return calls + + +def add_errors(rng, sequence): + return "".join(other_base(rng, base) if rng.random() < ERROR_RATE else base for base in sequence) + + +def sample_reads(rng, markers, sample): + reads = [] + for (_, haplotypes, calls) in markers: + for copy in calls[sample]: + amplicon = haplotypes[copy][0] + for _ in range(rng.randint(*READS_PER_COPY)): + read = add_errors(rng, amplicon) + if rng.random() < PARTIAL_FRACTION: + read = read[:rng.randint(TAG_LENGTH - 9, TAG_LENGTH - 1)] + reads.append(read) + for _ in range(rng.randint(*OFF_TARGET_READS)): + reads.append("".join(rng.choice("ACGT") for _ in range(AMPLICON_LENGTH))) + rng.shuffle(reads) + return reads + + +def write_fastq(path, sample_name, reads): + # mtime=0 and no stored filename keep regenerated files byte-identical + with open(path, "wb") as raw, gzip.GzipFile(filename="", mode="wb", fileobj=raw, mtime=0) as handle: + for number, read in enumerate(reads, start=1): + handle.write(f"@{sample_name}:{number}\n{read}\n+\n{'I' * len(read)}\n".encode()) + + +def main(): + rng = random.Random(SEED) + markers = [] + primer_rows, design_rows = [], [] + for index, (name, reverse, variant_position) in enumerate(MARKERS): + primer_row, design_row, haplotypes = build_marker(rng, index, name, reverse, variant_position) + primer_rows.append(primer_row) + design_rows.append(design_row) + markers.append((name, haplotypes, genotypes_for(index, haplotypes))) + + with (HERE / "primers.csv").open("w", newline="") as handle: + writer = csv.writer(handle, lineterminator="\n") + writer.writerow(["MarkerName", "LeftPrimerSequence", "RightPrimerSequence", + "VariantRef", "VariantAlt", "VariantPositionInTagSequence"]) + writer.writerows(primer_rows) + with (HERE / "probe_design.tsv").open("w", newline="") as handle: + writer = csv.writer(handle, delimiter="\t", lineterminator="\n") + writer.writerow(["MarkerName", "TargetSequence"]) + writer.writerows(design_rows) + + samples = [f"S{number:02d}" for number in range(1, N_SAMPLES + 1)] + with (HERE / "sample_names.csv").open("w", newline="") as handle: + writer = csv.writer(handle, lineterminator="\n") + writer.writerow(["fastq_id", "sample_name"]) + writer.writerows([sample, f"Line_{sample[1:]}"] for sample in samples) + + with (HERE / "truth_haplotypes.csv").open("w", newline="") as handle: + writer = csv.writer(handle, lineterminator="\n") + writer.writerow(["MarkerName", "Haplotype", "TagSequence", "Description"]) + for name, haplotypes, _ in markers: + for haplotype, (amplicon, description) in haplotypes.items(): + writer.writerow([name, haplotype, amplicon[:TAG_LENGTH], description]) + with (HERE / "truth_genotypes.csv").open("w", newline="") as handle: + writer = csv.writer(handle, lineterminator="\n") + writer.writerow(["Sample", "MarkerName", "Genotype"]) + for sample_index, sample in enumerate(samples): + for name, _, calls in markers: + writer.writerow([f"Line_{sample[1:]}", name, "/".join(calls[sample_index])]) + + fastq_dir = HERE / "fastq" + fastq_dir.mkdir(exist_ok=True) + for sample_index, sample in enumerate(samples): + write_fastq(fastq_dir / f"{sample}.fastq.gz", sample, sample_reads(rng, markers, sample_index)) + + +if __name__ == "__main__": + main() diff --git a/examples/madcats-demo/primers.csv b/examples/madcats-demo/primers.csv new file mode 100644 index 0000000..46db86b --- /dev/null +++ b/examples/madcats-demo/primers.csv @@ -0,0 +1,9 @@ +MarkerName,LeftPrimerSequence,RightPrimerSequence,VariantRef,VariantAlt,VariantPositionInTagSequence +Chr01_000412883,AGATATGCTGGGTAGAGGTC,AGGAACGGATCGAGGACTGC,G,C,34 +Chr02_001257730,GTCATACCAATCTACCCCCT,GCTCACTCCCTATCCACCTA,G,C,48 +Chr03_000086214,ATGCCAAGTTGGCACAGGGA,ATTTCACGTAACTCCATTAC,T,C,41 +Chr04_002318006,TCTGGAGGTGCGCGCAGGTT,GATGTACGTGCCCACTACTA,A,T,57 +Chr05_000970341,CGGCTCCACCTTAGCGATTC,CGCTTTTCGTATTGCGTAAA,G,C,29 +Chr06_001644492,CTATCGCGGTAGCAGCACTT,TAGCTTCGCGATGCTTCCAA,C,G,52 +Chr07_000203517,GTATGTTGGTTTGGCAACCG,ATTGATAAGGCTTGCAGTTA,T,C,63 +Chr08_001105968,CGGTACGATAGGTGCCTGAG,GTCTGAAAATTAGCGCATGA,G,A,45 diff --git a/examples/madcats-demo/probe_design.tsv b/examples/madcats-demo/probe_design.tsv new file mode 100644 index 0000000..9021377 --- /dev/null +++ b/examples/madcats-demo/probe_design.tsv @@ -0,0 +1,9 @@ +MarkerName TargetSequence +Chr01_000412883 AGACTTTCAAAGATATGCTGGGTAGAGGTCGAGGTTATTATTT[G/C]TTACCAATTCTCATTGTGTTTCGGAACTTGCGTTTTAGGTATGTCTTAGTGACTCTAAATACCAAGGCAGTCCTCGATCCGTTCCTAATAAGGAAT +Chr02_001257730 GTGATTCCCTGTCATACCAATCTACCCCCTGTTATGCGCGTTTGTCGTTAGACCAAT[G/C]TCAGCGCAGCGGCAGATCAAGCAGGAGGCGGAATGTAAACAGAAGGTATGCTTAGGTGGATAGGGAGTGAGCAACAAACGGA +Chr03_000086214 AACCACGGCCATTTCACGTAACTCCATTACTGGGGCCTCTTTGGCGTGCAATGAGCCCCGGTCCCAGTTGAATCGTTGCCCTGTACTAG[A/G]GGCAAACCGCCGCAGGTAGTTCCCTGTGCCAACTTGGCATGGGAGAAACG +Chr04_002318006 CCTCGGTTCCTCTGGAGGTGCGCGCAGGTTTAGTGATCTGGATCAGGCGTTTGAACAGGACTGGAC[A/T]ACGCTCCGATCAAGTACCTGGGGTGTGGATCATGGTCGGTGCATAGTAGTGGGCACGTACATCCTCCGTCGGT +Chr05_000970341 CCCCCAAGGCCGGCTCCACCTTAGCGATTCGCATAAAG[G/C]GACCATTCCCCTCGCTTTGCGTTACCAAGCACCTCTACGGGTATTGGTCACCGTTAGGATTGAGGAAAATATTTACGCAATACGAAAAGCGTCGTCTCGTA +Chr06_001644492 CCTCATTTATCTATCGCGGTAGCAGCACTTAATTACCAGCGGAGCACATAGTACTAGACAG[C/G]AGATGACCTAGCTTACAATTATCCCCCGTGCTAAGACCTCGCTACATATTGGAAGCATCGCGAAGCTACTCTAACAGC +Chr07_000203517 TTCAATAATAATTGATAAGGCTTGCAGTTAAGCTTTCCGATCATAGTTTTCGAGAAAGGGTTCAGCA[A/G]ATGACCCGACCTCTAAAGAGCGTCGGACCACGACGAGTGCGACGGTTGCCAAACCAACATACGTTATAGCAT +Chr08_001105968 CTAGCTAGCACGGTACGATAGGTGCCTGAGTGTGCTCGGGGACACGGAGTTTTA[G/A]TTGTGAGTGGTAGTTAAGAAGTGTGGACCAACCAGTAAGCATCAAAATACTAGTGTCATGCGCTAATTTTCAGACTCGTCCGCGA diff --git a/examples/madcats-demo/sample_names.csv b/examples/madcats-demo/sample_names.csv new file mode 100644 index 0000000..7ac76f3 --- /dev/null +++ b/examples/madcats-demo/sample_names.csv @@ -0,0 +1,13 @@ +fastq_id,sample_name +S01,Line_01 +S02,Line_02 +S03,Line_03 +S04,Line_04 +S05,Line_05 +S06,Line_06 +S07,Line_07 +S08,Line_08 +S09,Line_09 +S10,Line_10 +S11,Line_11 +S12,Line_12 diff --git a/examples/madcats-demo/truth_genotypes.csv b/examples/madcats-demo/truth_genotypes.csv new file mode 100644 index 0000000..00a51f7 --- /dev/null +++ b/examples/madcats-demo/truth_genotypes.csv @@ -0,0 +1,97 @@ +Sample,MarkerName,Genotype +Line_01,Chr01_000412883,Ref/Ref +Line_01,Chr02_001257730,Ref/RefMatch +Line_01,Chr03_000086214,Alt/Alt +Line_01,Chr04_002318006,AltMatch/Ref +Line_01,Chr05_000970341,Ref/Ref +Line_01,Chr06_001644492,Insertion/Ref +Line_01,Chr07_000203517,Alt/Alt +Line_01,Chr08_001105968,Ref/Ref +Line_02,Chr01_000412883,Alt/Ref +Line_02,Chr02_001257730,Alt/Alt +Line_02,Chr03_000086214,Alt/Ref +Line_02,Chr04_002318006,Ref/Ref +Line_02,Chr05_000970341,Alt/Ref +Line_02,Chr06_001644492,Alt/Alt +Line_02,Chr07_000203517,Alt/Alt +Line_02,Chr08_001105968,Ref/Ref +Line_03,Chr01_000412883,Alt/Alt +Line_03,Chr02_001257730,Alt/Ref +Line_03,Chr03_000086214,Ref/Ref +Line_03,Chr04_002318006,Alt/Ref +Line_03,Chr05_000970341,Alt/Alt +Line_03,Chr06_001644492,Alt/Alt +Line_03,Chr07_000203517,Alt/Ref +Line_03,Chr08_001105968,Ref/Ref +Line_04,Chr01_000412883,Alt/Ref +Line_04,Chr02_001257730,Ref/RefMatch +Line_04,Chr03_000086214,Alt/Ref +Line_04,Chr04_002318006,Alt/AltMatch +Line_04,Chr05_000970341,Alt/Alt +Line_04,Chr06_001644492,Alt/Insertion +Line_04,Chr07_000203517,Ref/Ref +Line_04,Chr08_001105968,Ref/Ref +Line_05,Chr01_000412883,Ref/Ref +Line_05,Chr02_001257730,Alt/Ref +Line_05,Chr03_000086214,Alt/Alt +Line_05,Chr04_002318006,Alt/Alt +Line_05,Chr05_000970341,Alt/Ref +Line_05,Chr06_001644492,Ref/Ref +Line_05,Chr07_000203517,Alt/Ref +Line_05,Chr08_001105968,Ref/Ref +Line_06,Chr01_000412883,Alt/Ref +Line_06,Chr02_001257730,Alt/Alt +Line_06,Chr03_000086214,Alt/Alt +Line_06,Chr04_002318006,Alt/Ref +Line_06,Chr05_000970341,Ref/Ref +Line_06,Chr06_001644492,Alt/Ref +Line_06,Chr07_000203517,Alt/Ref +Line_06,Chr08_001105968,Ref/Ref +Line_07,Chr01_000412883,Alt/Alt +Line_07,Chr02_001257730,Alt/RefMatch +Line_07,Chr03_000086214,Alt/Ref +Line_07,Chr04_002318006,AltMatch/Ref +Line_07,Chr05_000970341,Alt/Ref +Line_07,Chr06_001644492,Insertion/Ref +Line_07,Chr07_000203517,Ref/Ref +Line_07,Chr08_001105968,Ref/Ref +Line_08,Chr01_000412883,Alt/Alt +Line_08,Chr02_001257730,Alt/Ref +Line_08,Chr03_000086214,Ref/Ref +Line_08,Chr04_002318006,Alt/Ref +Line_08,Chr05_000970341,Alt/Ref +Line_08,Chr06_001644492,Ref/Ref +Line_08,Chr07_000203517,Alt/Ref +Line_08,Chr08_001105968,Ref/Ref +Line_09,Chr01_000412883,Alt/Ref +Line_09,Chr02_001257730,Ref/Ref +Line_09,Chr03_000086214,Alt/Ref +Line_09,Chr04_002318006,Alt/Ref +Line_09,Chr05_000970341,Ref/Ref +Line_09,Chr06_001644492,Alt/Ref +Line_09,Chr07_000203517,Alt/Alt +Line_09,Chr08_001105968,Ref/Ref +Line_10,Chr01_000412883,Ref/Ref +Line_10,Chr02_001257730,Alt/RefMatch +Line_10,Chr03_000086214,Alt/Ref +Line_10,Chr04_002318006,AltMatch/Ref +Line_10,Chr05_000970341,Alt/Ref +Line_10,Chr06_001644492,Alt/Insertion +Line_10,Chr07_000203517,Alt/Ref +Line_10,Chr08_001105968,Ref/Ref +Line_11,Chr01_000412883,Alt/Ref +Line_11,Chr02_001257730,Alt/Ref +Line_11,Chr03_000086214,Ref/Ref +Line_11,Chr04_002318006,Alt/Ref +Line_11,Chr05_000970341,Alt/Alt +Line_11,Chr06_001644492,Alt/Ref +Line_11,Chr07_000203517,Ref/Ref +Line_11,Chr08_001105968,Ref/Ref +Line_12,Chr01_000412883,Alt/Ref +Line_12,Chr02_001257730,Ref/Ref +Line_12,Chr03_000086214,Alt/Ref +Line_12,Chr04_002318006,Alt/Alt +Line_12,Chr05_000970341,Alt/Ref +Line_12,Chr06_001644492,Ref/Ref +Line_12,Chr07_000203517,Alt/Ref +Line_12,Chr08_001105968,Ref/Ref diff --git a/examples/madcats-demo/truth_haplotypes.csv b/examples/madcats-demo/truth_haplotypes.csv new file mode 100644 index 0000000..4911dda --- /dev/null +++ b/examples/madcats-demo/truth_haplotypes.csv @@ -0,0 +1,20 @@ +MarkerName,Haplotype,TagSequence,Description +Chr01_000412883,Ref,AGATATGCTGGGTAGAGGTCGAGGTTATTATTTGTTACCAATTCTCATTGTGTTTCGGAACTTGCGTTTTAGGTATGTCTT,designed Ref allele +Chr01_000412883,Alt,AGATATGCTGGGTAGAGGTCGAGGTTATTATTTCTTACCAATTCTCATTGTGTTTCGGAACTTGCGTTTTAGGTATGTCTT,designed Alt allele +Chr02_001257730,Ref,GTCATACCAATCTACCCCCTGTTATGCGCGTTTGTCGTTAGACCAATGTCAGCGCAGCGGCAGATCAAGCAGGAGGCGGAA,designed Ref allele +Chr02_001257730,Alt,GTCATACCAATCTACCCCCTGTTATGCGCGTTTGTCGTTAGACCAATCTCAGCGCAGCGGCAGATCAAGCAGGAGGCGGAA,designed Alt allele +Chr02_001257730,RefMatch,GTCATACCAATCTACCCCCTGTTATGCGCGTTTGTCGTTAGACCAATGTCAGCGCAGCGGCAGATCTAGCAGGAGGCGGAA,Ref allele plus a second SNP at tag position 67 +Chr03_000086214,Ref,ATGCCAAGTTGGCACAGGGAACTACCTGCGGCGGTTTGCCTCTAGTACAGGGCAACGATTCAACTGGGACCGGGGCTCATT,designed Ref allele +Chr03_000086214,Alt,ATGCCAAGTTGGCACAGGGAACTACCTGCGGCGGTTTGCCCCTAGTACAGGGCAACGATTCAACTGGGACCGGGGCTCATT,designed Alt allele +Chr04_002318006,Ref,TCTGGAGGTGCGCGCAGGTTTAGTGATCTGGATCAGGCGTTTGAACAGGACTGGACAACGCTCCGATCAAGTACCTGGGGT,designed Ref allele +Chr04_002318006,Alt,TCTGGAGGTGCGCGCAGGTTTAGTGATCTGGATCAGGCGTTTGAACAGGACTGGACTACGCTCCGATCAAGTACCTGGGGT,designed Alt allele +Chr04_002318006,AltMatch,TCTGGAGGTGCGCGCAGGTTTAGTTATCTGGATCAGGCGTTTGAACAGGACTGGACTACGCTCCGATCAAGTACCTGGGGT,Alt allele plus a second SNP at tag position 25 +Chr05_000970341,Ref,CGGCTCCACCTTAGCGATTCGCATAAAGGGACCATTCCCCTCGCTTTGCGTTACCAAGCACCTCTACGGGTATTGGTCACC,designed Ref allele +Chr05_000970341,Alt,CGGCTCCACCTTAGCGATTCGCATAAAGCGACCATTCCCCTCGCTTTGCGTTACCAAGCACCTCTACGGGTATTGGTCACC,designed Alt allele +Chr06_001644492,Ref,CTATCGCGGTAGCAGCACTTAATTACCAGCGGAGCACATAGTACTAGACAGCAGATGACCTAGCTTACAATTATCCCCCGT,designed Ref allele +Chr06_001644492,Alt,CTATCGCGGTAGCAGCACTTAATTACCAGCGGAGCACATAGTACTAGACAGGAGATGACCTAGCTTACAATTATCCCCCGT,designed Alt allele +Chr06_001644492,Insertion,CTATCGCGGTAGCAGCACTTAATTACCAGCCGGAGCACATAGTACTAGACAGCAGATGACCTAGCTTACAATTATCCCCCG,Ref allele with an undesigned 1-bp insertion near tag position 30 +Chr07_000203517,Ref,GTATGTTGGTTTGGCAACCGTCGCACTCGTCGTGGTCCGACGCTCTTTAGAGGTCGGGTCATTTGCTGAACCCTTTCTCGA,designed Ref allele +Chr07_000203517,Alt,GTATGTTGGTTTGGCAACCGTCGCACTCGTCGTGGTCCGACGCTCTTTAGAGGTCGGGTCATCTGCTGAACCCTTTCTCGA,designed Alt allele +Chr08_001105968,Ref,CGGTACGATAGGTGCCTGAGTGTGCTCGGGGACACGGAGTTTTAGTTGTGAGTGGTAGTTAAGAAGTGTGGACCAACCAGT,designed Ref allele +Chr08_001105968,Alt,CGGTACGATAGGTGCCTGAGTGTGCTCGGGGACACGGAGTTTTAATTGTGAGTGGTAGTTAAGAAGTGTGGACCAACCAGT,designed Alt allele diff --git a/madc_pipeline.py b/madc_pipeline.py deleted file mode 100644 index 9adfdde..0000000 --- a/madc_pipeline.py +++ /dev/null @@ -1,11 +0,0 @@ -"""Compatibility shim: the tool is now MADCATS (madcats.py). - -Everything that imported or invoked madc_pipeline continues to work; -new code should import madcats directly. -""" -from madcats import * # noqa: F401,F403 -from madcats import main, __version__ # noqa: F401 - -if __name__ == "__main__": - import sys - sys.exit(main()) diff --git a/madcats.py b/madcats.py index aa92293..532075c 100644 --- a/madcats.py +++ b/madcats.py @@ -152,7 +152,7 @@ # ============================================================================ -__version__ = "0.12.1" +__version__ = "0.13.0" # Increment only when the serialized per-sample cache or raw-read processing # changes. Output-row filtering does not invalidate reusable FASTQ counts. @@ -5048,6 +5048,14 @@ def proper_fraction(value: str) -> float: return parsed +def partial_read_mode(value: str) -> str: + """``--partial-read-mode``: on/off, still accepting the 0.12 spelling ``assign``.""" + if value == "assign": + print("warning: --partial-read-mode assign is now spelled 'on'", file=sys.stderr) + return "on" + return value + + def ratio_at_least_one(value: str) -> float: parsed = float(value) if parsed < 1: @@ -5071,10 +5079,24 @@ def _add_fastq_arguments(parser: argparse.ArgumentParser, directory_help: str = def _add_design_command(subparsers) -> None: - design = subparsers.add_parser("design", help="reconstruct Ref/Alt tags from DArT design tables") - design.add_argument("--primers", required=True, help="primer-sequence CSV") - design.add_argument("--probe-design", required=True, help="probe/marker design TSV or CSV") - design.add_argument("--tag-length", type=positive_int, default=TAG_LENGTH_DEFAULT) + design = subparsers.add_parser( + "design", + help="reconstruct each marker's Ref/Alt tags from DArT primer and probe-design files", + description=( + "Reconstruct each marker's canonical Ref and Alt tags from a DArT primer-sequence CSV and " + "probe-design file, and write them as a reference-tag CSV for generate. Uses no sequencing data." + ), + ) + design.add_argument("--primers", required=True, help="DArT primer-sequence CSV") + design.add_argument( + "--probe-design", required=True, help="probe-design TSV or CSV with one bracketed variant per marker, e.g. [C/T]" + ) + design.add_argument( + "--tag-length", + type=positive_int, + default=TAG_LENGTH_DEFAULT, + help=f"tag length to build, e.g. 54, 81 or 109 (default: {TAG_LENGTH_DEFAULT})", + ) design.add_argument("--output", required=True, help="reference-tag CSV to write") design.set_defaults(func=run_design) @@ -5082,9 +5104,15 @@ def _add_design_command(subparsers) -> None: def _add_infer_design_command(subparsers) -> None: infer = subparsers.add_parser( "infer-design", - help="infer a reference-tag panel from the probe design plus FASTQ evidence when the primer file is unavailable", + help="infer the reference-tag panel from the probe design and the reads when the primer file is unavailable", + description=( + "Infer each marker's Ref and Alt tags from the probe design and the reads, for panels whose primer " + "file is not available. Reads vote for each locus's tag window: exact prefix matches decide first, " + "and a fuzzy pass runs only for loci with no exact votes. Loci that cannot be resolved are left out " + "and listed, with the reason, in the QC report." + ), ) - infer.add_argument("--probe-design", required=True, help="probe/marker design TSV or CSV") + infer.add_argument("--probe-design", required=True, help="probe-design TSV or CSV with bracketed variants") _add_fastq_arguments(infer) infer.add_argument( "--leading-anchor", @@ -5101,7 +5129,7 @@ def _add_infer_design_command(subparsers) -> None: "--min-identity", type=probability, default=0.80, - help="identity envelope for the fuzzy pass, matching generate (default: 0.80)", + help="minimum identity for fuzzy-pass votes, as in generate (default: 0.80)", ) infer.add_argument( "--dominance-ratio", @@ -5113,7 +5141,7 @@ def _add_infer_design_command(subparsers) -> None: "--fuzzy-min-reads", type=nonnegative_int, default=30, - help="identity-envelope votes required to anchor a locus with no exact support; 0 disables the fuzzy pass (default: 30)", + help="fuzzy-pass votes needed to resolve a locus with no exact votes; 0 disables the fuzzy pass (default: 30)", ) infer.add_argument( "--fuzzy-min-samples", @@ -5131,7 +5159,7 @@ def _add_infer_design_command(subparsers) -> None: "--anchor-max-length", type=positive_int, default=30, - help="longest pseudo-primer anchor; also capped below each variant (default: 30)", + help="longest pseudo-primer anchor; an anchor always ends before the variant (default: 30)", ) _add_jobs_argument(infer) infer.add_argument("--output", required=True, help="reference-tag CSV to write") @@ -5142,7 +5170,14 @@ def _add_infer_design_command(subparsers) -> None: def _add_merge_pairs_command(subparsers) -> None: merge = subparsers.add_parser( "merge-pairs", - help="stitch overlapping paired-end reads into single fragment reads (R2 reverse-complemented; read-through adapter tails removed by the overlap geometry)", + help="merge overlapping paired-end reads into single fragment reads", + description=( + "Merge each overlapping R1/R2 pair into one fragment read, before normalize-reads or generate. R2 is " + "reverse-complemented and the best-scoring overlap is used (overlap length with a penalty for " + "each mismatch). Adapter read-through falls outside the " + "overlap and never reaches the merged read, so no adapter trimming is needed. Use this only when " + "mates overlap." + ), ) merge.add_argument("--fastq-dir", required=True, help="folder holding paired FASTQ(.gz) files") merge.add_argument("--r1-suffix", default="_R1.fastq.gz", help="filename suffix of the R1 mate (default: _R1.fastq.gz)") @@ -5158,8 +5193,8 @@ def _add_merge_pairs_command(subparsers) -> None: merge.add_argument( "--unmerged", choices=("r1", "pool", "drop"), default="r1", help=( - "unmergeable pairs: keep the R1 (default), keep both mates (pool — for dual-anchored " - "designs whose mates cover opposite locus halves), or drop" + "pairs that do not merge: keep the R1 (default), keep both mates (pool; for dual-anchored " + "designs whose mates cover opposite halves of the locus), or drop the pair" ), ) _add_jobs_argument(merge) @@ -5170,12 +5205,14 @@ def _add_merge_pairs_command(subparsers) -> None: def _add_window_scan_command(subparsers) -> None: scan = subparsers.add_parser( "window-scan", - help=( - "sweep candidate tag lengths on a sample subset and report the selection curves " - "(retention, cleanliness, information) with a recommended length" + help="try several tag lengths on a subset of samples and recommend one", + description=( + "Run normalize-reads and generate at each candidate tag length on a subset of samples, measure " + "read retention, cleanliness and haplotype information, and recommend a length. Use it before " + "the first normalize-reads run of a new panel." ), ) - scan.add_argument("--probe-design", required=True, help="DArT probe/marker design TSV or CSV") + scan.add_argument("--probe-design", required=True, help="probe-design TSV or CSV with bracketed variants") _add_fastq_arguments(scan) scan.add_argument( "--lengths", required=True, @@ -5197,8 +5234,8 @@ def _add_window_scan_command(subparsers) -> None: "--scan-draws", type=positive_int, default=1, help=( "independent sample draws to sweep (seeds --scan-seed, +1, +2, ...; default: 1); with more than one, " - "each draw's recommendation and their agreement are reported, so a recommendation that depends " - "on the subset is visible without external truth" + "each draw's recommendation and their agreement are reported, so you can see whether the " + "recommendation depends on which samples were scanned" ), ) scan.add_argument( @@ -5213,12 +5250,21 @@ def _add_window_scan_command(subparsers) -> None: "informative eligible length is recommended (default: 0.10)" ), ) - scan.add_argument("--anchor-length", type=positive_int, default=12) - scan.add_argument("--kmer-size", type=positive_int, default=21) - scan.add_argument("--str-penalty", type=float, default=0.002) + scan.add_argument( + "--anchor-length", type=positive_int, default=12, + help="as in normalize-reads; used for the coverage pass and every candidate run (default: 12)", + ) + scan.add_argument( + "--kmer-size", type=positive_int, default=21, + help="as in normalize-reads; used for the coverage pass and every candidate run (default: 21)", + ) + scan.add_argument( + "--str-penalty", type=float, default=0.002, + help="as in normalize-reads; used for the coverage pass and every candidate run (default: 0.002)", + ) scan.add_argument( "--group-routing", choices=("nearest", "covered"), default="nearest", - help="passed through to the candidate normalize runs (use covered for merged tiled input)", + help="passed to the candidate normalize-reads runs; use covered for merged, tiled input (default: nearest)", ) _add_jobs_argument(scan) scan.add_argument("--out-dir", required=True, help="folder for per-length runs and the scan report") @@ -5228,15 +5274,21 @@ def _add_window_scan_command(subparsers) -> None: def _add_normalize_reads_command(subparsers) -> None: normalize = subparsers.add_parser( "normalize-reads", - help="normalize bidirectional, variable-start amplicon reads (e.g. AgriSeq) into pipeline-ready pseudo-tag FASTQs plus a window panel", + help="orient and window reads that do not start at a primer (e.g. AgriSeq), ready for generate", + description=( + "For platforms such as AgriSeq (Ion Torrent) whose reads do not start at a fixed primer. Each read " + "is located against the probe-design targets, oriented forward, and cut to a fixed-length window " + "containing the variant. Writes one FASTQ per sample plus a window reference-tags CSV, ready for " + "generate. Pass the first run's CSV to later runs with --reuse-tags." + ), ) - normalize.add_argument("--probe-design", required=True, help="probe/marker design TSV or CSV with bracketed variants") + normalize.add_argument("--probe-design", required=True, help="probe-design TSV or CSV with bracketed variants") _add_fastq_arguments(normalize, directory_help="folder containing one raw FASTQ(.gz) per sample") normalize.add_argument( "--tag-length", type=positive_int, default=41, - help="excised window length; shorter windows suffer fewer indel losses on indel-prone chemistries (default: 41)", + help="window length to cut; shorter windows lose fewer reads to indels on indel-prone chemistries (default: 41)", ) normalize.add_argument( "--window-placement", @@ -5270,11 +5322,11 @@ def _add_normalize_reads_command(subparsers) -> None: choices=("nearest", "covered"), default="nearest", help=( - "how a read covering several member windows of one LocusGroup is routed: 'nearest' feeds " - "exactly one tile (window center nearest the read center; default), 'covered' emits the " - "excision at every member window the read fully covers — use with merged full-span fragments, " - "where nearest-routing is winner-take-all and starves the other tiles (each window still " - "sees each fragment at most once in both modes)" + "how a read covering several windows of one LocusGroup is routed: 'nearest' (default) sends it " + "to the one tile whose window center is nearest the read's center; 'covered' cuts a window at " + "every tile the read fully covers. Use covered with merged full-span fragments: they all share " + "the amplicon's center, so with nearest one tile would get almost every read. In both modes " + "each window counts each fragment at most once" ), ) normalize.add_argument( @@ -5283,7 +5335,8 @@ def _add_normalize_reads_command(subparsers) -> None: default=0.002, help=( "auto placement: coverage penalty per tandem-repeat bp (unit sizes 1-3) beyond 4 in a candidate " - "window; steers windows off repeat runs that cause allele-biased read dropout (0 disables; default: 0.002)" + "window; steers windows away from repeats where reads carrying one allele are lost more often " + "(0 disables; default: 0.002)" ), ) normalize.add_argument( @@ -5301,14 +5354,27 @@ def _add_normalize_reads_command(subparsers) -> None: def _add_generate_command(subparsers) -> None: """``generate``'s options, grouped by the pipeline stage each one feeds.""" - generate = subparsers.add_parser("generate", help="generate a MADC-compatible count table") + generate = subparsers.add_parser( + "generate", + help="count alleles in every sample and write the MADC table", + description=( + "Assign reads to markers, count exact tag sequences per sample, report new microhaplotypes that " + "pass the evidence rules, remove sequencing errors by denoising, and write a raw MADC table plus " + "a QC/provenance JSON. Run one generate over the whole cohort: the reproducibility route counts " + "supporting samples within the run. The raw MADC goes through HapApp before BIGr and other " + "downstream tools." + ), + ) panel = generate.add_argument_group("panel and samples (stage 1)") _add_fastq_arguments(panel) - panel.add_argument("--primers", help="DArT primer-sequence CSV") - panel.add_argument("--probe-design", help="DArT probe/marker design TSV or CSV") + panel.add_argument("--primers", help="DArT primer-sequence CSV (use with --probe-design)") + panel.add_argument("--probe-design", help="DArT probe-design TSV or CSV (use with --primers)") panel.add_argument( "--reference-tags", - help="platform-neutral reference-tag CSV made by the design command; replaces --primers/--probe-design", + help=( + "reference-tag CSV from design, infer-design or normalize-reads, or curated by hand; " + "use instead of --primers/--probe-design" + ), ) panel.add_argument( "--target-metadata", @@ -5328,14 +5394,22 @@ def _add_generate_command(subparsers) -> None: "FASTQ ids, matching the raw-MADC layout" ), ) - panel.add_argument("--tag-length", type=positive_int, default=TAG_LENGTH_DEFAULT) + panel.add_argument( + "--tag-length", + type=positive_int, + default=TAG_LENGTH_DEFAULT, + help=f"tag length; must match the panel, e.g. 41 after normalize-reads (default: {TAG_LENGTH_DEFAULT})", + ) - assignment = generate.add_argument_group("read assignment (stage 2; these settings key the cache)") + assignment = generate.add_argument_group("read assignment (stage 2)") assignment.add_argument( "--min-identity", type=probability, default=0.80, - help="minimum substitution-only identity to Ref or Alt tag; a gapped reading spends the same edit budget (default: 0.80)", + help=( + "minimum identity to the marker's Ref or Alt tag, counting substitutions only; a gapped " + "alignment uses the same edit limit (default: 0.80)" + ), ) assignment.add_argument( "--leading-anchor", @@ -5355,9 +5429,8 @@ def _add_generate_command(subparsers) -> None: choices=(0, 1), default=1, help=( - "allowed mismatches while assigning a left primer/anchor (default: 1; every " - "validated benchmark ran with 1 — exact-only matching silently drops Match-haplotype " - "reads whose distinguishing SNP falls inside the anchor prefix)" + "mismatches allowed when matching a read's start to a primer/anchor (default: 1). With 0, " + "reads from Match haplotypes whose distinguishing SNP falls inside the anchor are lost" ), ) assignment.add_argument( @@ -5365,26 +5438,27 @@ def _add_generate_command(subparsers) -> None: choices=("on", "off"), default="on", help=( - "offer reads that fail substitution-only identity at every candidate one gapped reading " - "under the same edit budget, so undesigned indel haplotypes are counted instead of dropped " - "(default: on; 'off' is the 0.9 path)" + "give reads that fail the substitution-only identity check for every marker one gapped " + "alignment under the same edit limit, so undesigned indels are counted instead of dropped " + "(default: on)" ), ) assignment.add_argument( "--max-indel-length", type=positive_int, default=GAP_BAND_DEFAULT, - help="widest gap (alignment band) a gapped reading may use (default: 10)", + help=f"widest gap (alignment band) a gapped alignment may use (default: {GAP_BAND_DEFAULT})", ) assignment.add_argument( "--partial-read-mode", - choices=("off", "assign"), - default="assign", + type=partial_read_mode, + choices=("on", "off"), + default="on", help=( - "handling of primer-assigned reads shorter than --tag-length: assign credits a truncated read " - "to the one learned full-length allele it is a prefix of, and only when it is longer than every " - "prefix shared between that locus's alleles, so it can support but never invent an allele or " - "shift an allele ratio (default: assign); off drops such reads" + "handling of reads shorter than --tag-length: on adds such a read to the one full-length " + "allele whose start it matches exactly, only when that allele already passes the evidence rules " + "and the read is longer than every prefix shared between that locus's alleles, so partial reads " + "cannot create an allele or change an allele ratio (default: on); off drops them" ), ) assignment.add_argument( @@ -5392,32 +5466,39 @@ def _add_generate_command(subparsers) -> None: type=positive_int, default=None, help=( - "minimum read length a partial read must have to be credited (--partial-read-mode assign) " - "(default: --tag-length minus 9, so at most 9 unobserved tail bases)" + "minimum length for a partial read to be counted with --partial-read-mode on " + "(default: --tag-length minus 9)" + ), + ) + assignment.add_argument( + "--cache-dir", + help=( + "folder for a reusable per-sample cache, so re-runs skip the read pass. The cache is rebuilt " + "when the panel or --tag-length, --min-identity, --leading-anchor, --primer-mismatches, " + "--gapped-rescue, --max-indel-length or --min-partial-length changes" ), ) - assignment.add_argument("--cache-dir", help="optional reusable per-sample count cache") _add_jobs_argument(assignment) - evidence = generate.add_argument_group("evidence (stage 4: when a novel sequence earns a row)") + evidence = generate.add_argument_group("evidence (stage 4: which new sequences are reported)") evidence.add_argument( "--min-global-count", type=positive_int, default=1, - help="minimum project-wide reads for a de novo discovered allele (default: 1)", + help="minimum project-wide reads for a newly discovered allele (default: 1)", ) evidence.add_argument( "--reproducible-min-reads", type=positive_int, default=10, - help="reproducibility evidence route: per-sample reads used by --reproducible-min-samples (default: 10)", + help="reproducibility route: reads per sample counted toward --reproducible-min-samples (default: 10)", ) evidence.add_argument( "--reproducible-min-samples", type=nonnegative_int, default=2, help=( - "reproducibility evidence route: a novel allele needs --reproducible-min-reads reads in at " + "reproducibility route: a new allele needs at least --reproducible-min-reads reads in at " "least this many samples; 0 disables the route (default: 2)" ), ) @@ -5426,9 +5507,9 @@ def _add_generate_command(subparsers) -> None: type=nonnegative_int, default=30, help=( - "depth evidence route: reads within a single sample that establish a novel allele on their " - "own, so private alleles in unreplicated designs remain discoverable; 0 disables the route " - "(default: 30)" + "depth route: reads within one sample that are enough on their own to report a new allele, so " + "alleles found in only one sample can still be reported when samples are not replicated; " + "0 disables the route (default: 30)" ), ) evidence.add_argument( @@ -5436,16 +5517,16 @@ def _add_generate_command(subparsers) -> None: type=int, choices=QUALITY_THRESHOLDS, default=20, - help="Phred score every tag base must reach for a read to count as clean on the quality route (default: 20)", + help="Phred score every tag base must reach for a read to count as high-quality on the quality route (default: 20)", ) evidence.add_argument( "--quality-min-clean-reads", type=nonnegative_int, default=0, help=( - "quality evidence route: project-wide clean reads that establish a novel allele on their own; " - "0 disables the route (default: 0). Whole-tag clean-read requirements veto sequences with one " - "systematically low-quality base position, so reproducibility is the primary route." + "quality route: project-wide high-quality reads that are enough on their own to report a new " + "allele; 0 disables the route (default: 0). Off by default because one systematically " + "low-quality base position would block real alleles; reproducibility is the main route" ), ) evidence.add_argument( @@ -5453,10 +5534,11 @@ def _add_generate_command(subparsers) -> None: type=proper_fraction, default=0.0, help=( - "minimum share of a sample's locus reads before that sample can support the reproducibility " - "or depth route; 0 disables the floor (default: 0). Consider around 0.01 when relying on " - "the depth route in unreplicated designs; note co-amplified paralogs deflate real alleles' " - "shares, so a floor trades recall for precision" + "minimum share of a sample's reads at the locus before that sample can count toward the " + "reproducibility or depth route; 0 disables the floor (default: 0). Use 0.10 at very high depth " + "(above ~500x per sample per marker, e.g. AgriSeq), or about 0.01 for unreplicated designs that " + "rely on the depth route. Otherwise leave it off: co-amplified paralogs lower real alleles' " + "shares, so a floor can drop them" ), ) @@ -5466,7 +5548,10 @@ def _add_generate_command(subparsers) -> None: type=int, choices=(0, 1, 2), default=1, - help="per-sample abundance-skew denoising distance (default: 1)", + help=( + "per-sample denoising: the most differing bases allowed between a rare sequence and the more " + "abundant parent it is merged into; 0 disables (default: 1)" + ), ) denoising.add_argument( "--global-denoise-distance", @@ -5474,12 +5559,16 @@ def _add_generate_command(subparsers) -> None: choices=(0, 1, 2), default=1, help=( - "pooled same-allele fold: an evidence-failing sequence within this edit distance of a protected " - "parent that carries the same target allele joins it in every sample; 0 disables (default: 1)" + "pooled same-allele merge: a sequence that failed the evidence rules and differs at no more " + "than this many bases from a protected parent carrying the same target allele is merged into " + "it in every sample; 0 disables (default: 1)" ), ) denoising.add_argument( - "--denoise-alpha", type=float, default=2.0, help="UNOISE-style abundance-skew exponent (default: 2.0)" + "--denoise-alpha", + type=float, + default=2.0, + help="UNOISE-style abundance-skew exponent: a parent must be at least 2^(alpha*d + 1) times more abundant (default: 2.0)", ) output = generate.add_argument_group("output (stages 6-7)") @@ -5490,21 +5579,24 @@ def _add_generate_command(subparsers) -> None: choices=("on", "off"), default="on", help=( - "annotate every row in two never-computed statistics columns and a sidecar: PICSnp becomes " + "annotate every row in two otherwise unused statistics columns and a sidecar file: PICSnp becomes " "Variants (the row's events against the Ref tag in 1-based tag coordinates, e.g. " "55:C>CAAAA;91:A>T;clip=4), AvgPIC becomes MarkerType ('.', 'snp', 'indel' or 'both' relative " "to the nearest canonical), and OUTPUT.variants.csv lists one event per line; 'off' keeps the " "exact vendor header and writes no sidecar (default: on)" ), ) - output.add_argument("--number-match-ids", action="store_true", help="append stable per-file serials to Match IDs") + output.add_argument( + "--number-match-ids", + action="store_true", + help="append a stable serial to each RefMatch/AltMatch/Other ID (e.g. |RefMatch_0001) so every AlleleID is unique", + ) output.add_argument( "--dart-preamble", action="store_true", help=( - "prepend the vendor-style seven-row starred preamble (service label, plate rows, and the " - "targetid row carrying the FASTQ-derived sample ids) so preamble-expecting consumers such as " - "HapApp accept the file" + "prepend the vendor-style seven-row starred preamble (service label, plate rows, and a " + "targetid row with the FASTQ-derived sample ids) that HapApp expects" ), ) output.add_argument( @@ -5515,16 +5607,46 @@ def _add_generate_command(subparsers) -> None: def _add_validate_command(subparsers) -> None: - validate = subparsers.add_parser("validate", help="compare generated and expected MADC tables") - validate.add_argument("--generated", required=True) - validate.add_argument("--expected", required=True) - validate.add_argument("--report", help="optional JSON report to write") + validate = subparsers.add_parser( + "validate", + help="compare a generated MADC with a reference MADC", + description=( + "Compare a generated MADC with a reference MADC, usually the vendor's delivery for the same " + "FASTQs: row precision and recall by allele category, classification agreement on shared rows, " + "and count agreement across shared cells. The report is printed as JSON. The reference MADC is " + "used only for this comparison." + ), + ) + validate.add_argument("--generated", required=True, help="MADC table written by generate") + validate.add_argument( + "--expected", + required=True, + help="reference MADC to compare against, with or without the vendor's starred preamble rows", + ) + validate.add_argument("--report", help="also write the JSON report to this file") validate.set_defaults(func=run_validate) def build_parser() -> argparse.ArgumentParser: parser = argparse.ArgumentParser( - description="Generate and validate MADC-compatible microhaplotype read-count tables from primer-started FASTQ files." + prog="madcats", + formatter_class=argparse.RawDescriptionHelpFormatter, + description=( + "MADCATS (Marker-Anchored Discovery and Counting of Alleles from Targeted\n" + "Sequencing) turns targeted-sequencing FASTQs into raw MADC tables: the\n" + "microhaplotypes at each marker and their read counts in every sample. It\n" + "counts the alleles each panel was designed for and discovers new ones." + ), + epilog=( + "typical workflow:\n" + " 1. build the panel design, or infer-design without a primer file\n" + " 2. prepare reads if needed merge-pairs, window-scan, normalize-reads\n" + " 3. count generate\n" + " 4. check (optional) validate\n" + "\n" + "Run 'madcats --help' for a command's options.\n" + "Documentation: https://madcats.readthedocs.io" + ), ) parser.add_argument("--version", action="version", version=f"%(prog)s {__version__}") subparsers = parser.add_subparsers(dest="command", required=True) diff --git a/pyproject.toml b/pyproject.toml index 20b2f62..0a4a751 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -48,7 +48,7 @@ Changelog = "https://madcats.readthedocs.io/en/latest/development/changelog.html madcats = "madcats:main" [tool.setuptools] -py-modules = ["madcats", "madc_pipeline"] +py-modules = ["madcats"] [tool.setuptools.dynamic] version = { attr = "madcats.__version__" } diff --git a/scripts/analyze_count_vector_matches.py b/scripts/analyze_count_vector_matches.py index 3501f29..61fbfb5 100644 --- a/scripts/analyze_count_vector_matches.py +++ b/scripts/analyze_count_vector_matches.py @@ -13,7 +13,7 @@ from typing import Optional sys.path.insert(0, str(Path(__file__).resolve().parents[1])) -from madc_pipeline import hamming_distance, natural_key, read_madc +from madcats import hamming_distance, natural_key, read_madc def parse_args() -> argparse.Namespace: diff --git a/scripts/analyze_dart_only.py b/scripts/analyze_dart_only.py index 5b2ca2b..5b79c55 100644 --- a/scripts/analyze_dart_only.py +++ b/scripts/analyze_dart_only.py @@ -13,7 +13,7 @@ from typing import Optional sys.path.insert(0, str(Path(__file__).resolve().parents[1])) -from madc_pipeline import QUALITY_THRESHOLDS, allele_type, compatible_mismatches, read_madc +from madcats import QUALITY_THRESHOLDS, allele_type, compatible_mismatches, read_madc def parse_args() -> argparse.Namespace: diff --git a/scripts/analyze_extra_quality.py b/scripts/analyze_extra_quality.py index 21aa9bb..ef5226a 100644 --- a/scripts/analyze_extra_quality.py +++ b/scripts/analyze_extra_quality.py @@ -12,7 +12,7 @@ from pathlib import Path sys.path.insert(0, str(Path(__file__).resolve().parents[1])) -from madc_pipeline import QUALITY_THRESHOLDS, allele_type, hamming_distance, read_madc +from madcats import QUALITY_THRESHOLDS, allele_type, hamming_distance, read_madc def parse_args() -> argparse.Namespace: diff --git a/scripts/madc_to_reference_tags.py b/scripts/madc_to_reference_tags.py index 95e6bbe..5d5ac7a 100644 --- a/scripts/madc_to_reference_tags.py +++ b/scripts/madc_to_reference_tags.py @@ -19,7 +19,7 @@ from pathlib import Path sys.path.insert(0, str(Path(__file__).resolve().parents[1])) -from madc_pipeline import Locus, allele_type, read_madc, write_reference_tags +from madcats import Locus, allele_type, read_madc, write_reference_tags # Why a marker cannot become a locus, in the words a user needs to act on it. diff --git a/scripts/validate_inferred_design.py b/scripts/validate_inferred_design.py index 4c5d666..5296aad 100644 --- a/scripts/validate_inferred_design.py +++ b/scripts/validate_inferred_design.py @@ -16,7 +16,7 @@ from pathlib import Path sys.path.insert(0, str(Path(__file__).resolve().parents[1])) -from madc_pipeline import load_reference_tags +from madcats import load_reference_tags def parse_args() -> argparse.Namespace: diff --git a/tests/test_docs.py b/tests/test_docs.py new file mode 100644 index 0000000..8fd5296 --- /dev/null +++ b/tests/test_docs.py @@ -0,0 +1,93 @@ +"""Keep the command reference pages in step with the command-line parser. + +Every option a subcommand accepts must appear in a table on its page in +docs/src/commands/, and every option those tables list must exist. Where a +table gives a choice set or a plain default, it must match the parser too. +""" +import argparse +import re +import unittest +from pathlib import Path + +import madcats + +COMMANDS_DIR = Path(__file__).resolve().parents[1] / "docs" / "src" / "commands" +OPTION = re.compile(r"`(--[a-z0-9-]+)") +CHOICES = re.compile(r"\{([^}]*)\}") + + +def subcommand_parsers(): + parser = madcats.build_parser() + for action in parser._actions: + if isinstance(action, argparse._SubParsersAction): + return action.choices + raise AssertionError("madcats has no subcommands") + + +def documented_rows(page: Path): + """Yield (options, first cell, default cell or None) for each table row naming options.""" + for line in page.read_text().splitlines(): + if not line.startswith("| `--"): + continue + cells = [cell.strip() for cell in re.split(r"(?[ACGT]{2}") + + def test_counts_follow_the_simulated_genotypes(self): + samples = [column for column in self.rows[0] if column.startswith("Line_")] + self.assertEqual(len(samples), 12) + names = {(row["MarkerName"], row["TagSequence"]): row["Haplotype"] for row in self.truth} + counts = defaultdict(lambda: defaultdict(int)) + for row in self.rows: + haplotype = names[(row["CloneID"], row["AlleleSequence"])] + for sample in samples: + counts[(sample, row["CloneID"])][haplotype] += int(row[sample]) + for call in read_csv(EXAMPLE / "truth_genotypes.csv"): + observed = counts[(call["Sample"], call["MarkerName"])] + total = sum(observed.values()) + carried = call["Genotype"].split("/") + for haplotype, reads in observed.items(): + share = reads / total + with self.subTest(sample=call["Sample"], marker=call["MarkerName"], haplotype=haplotype): + if haplotype not in carried: + self.assertLess(share, 0.05) + elif carried.count(haplotype) == 1: + self.assertTrue(0.2 < share < 0.8, share) + else: + self.assertGreater(share, 0.95) + + +if __name__ == "__main__": + unittest.main() diff --git a/tests/test_madc_pipeline.py b/tests/test_madcats.py similarity index 99% rename from tests/test_madc_pipeline.py rename to tests/test_madcats.py index 35dc768..bc51093 100644 --- a/tests/test_madc_pipeline.py +++ b/tests/test_madcats.py @@ -197,7 +197,7 @@ def test_partial_assign_never_invents_an_allele(self): updated, summary = mp.apply_partial_reads( [result], [locus], - mode="assign", + mode="on", rule=rule(clean_reads=5), ) self.assertEqual(summary["reads_unmapped"], 5) @@ -233,7 +233,7 @@ def _apply_partials(self, locus, counts, quality, partials): partial_quality_support={(0, partial): (n, n, n, n) for partial, n in partials.items()}, qc={}, ) - return mp.apply_partial_reads([result], [locus], mode="assign", rule=rule(clean_reads=5)) + return mp.apply_partial_reads([result], [locus], mode="on", rule=rule(clean_reads=5)) def test_partial_assign_locus_guard_drops_both_alleles_when_a_tail_variant_exists(self): # Ref, Alt and a RefMatch differing from Ref only at positions 107-109: @@ -277,7 +277,9 @@ def test_partial_shorter_than_locus_shared_prefix_is_never_credited(self): def test_partial_defaults_resolve_from_tag_length(self): parser = mp.build_parser() args = parser.parse_args(["generate", "--fastq-dir", "x", "--reference-tags", "y", "--tag-length", "109", "--output", "z"]) - self.assertEqual(args.partial_read_mode, "assign") + self.assertEqual(args.partial_read_mode, "on") + legacy = parser.parse_args(["generate", "--fastq-dir", "x", "--reference-tags", "y", "--output", "z", "--partial-read-mode", "assign"]) + self.assertEqual(legacy.partial_read_mode, "on") self.assertIsNone(args.min_partial_length) self.assertEqual(max(1, args.tag_length - mp.PARTIAL_TAIL_DEFAULT), 100)