Fixed large file cleanup from BLAST and updated config files

CHANGELOG.md

@@ -5,6 +5,18 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
+## [1.0.0b4] - 03-03-2026
+
+### Changed
+
+- **Automatic BLAST archive cleanup (`specificity.py`)**: The BLAST archive file (outfmt 11, can be 20 GB+) is now automatically deleted after the tabular output is produced. Wrapped in `try/finally` so the archive is removed even if an error occurs during the run or reformat step.
+- **Extracted init wizard and Docker wrapper from `cli.py`**: Split two self-contained blocks into their own modules to improve maintainability. `cli_init_wizard.py` contains `_is_interactive()`, `_prompt_path()`, and `_run_init_wizard()`. `cli_docker.py` contains the `docker` command, which registers itself on the shared Typer app via side-effect import. `cli.py` drops from 989 to 666 lines; no public API changes.
+
+### Added
+
+- **Off-target hard filter (`specificity.py`)**: New `filter_offtarget_pairs()` removes primer pairs with BLAST off-target amplicons before multiplex optimization. Per junction, pairs with zero off-targets are kept; when all pairs have off-targets, the pair with fewest is retained as a fallback. Called automatically from `run_pipeline()` after `run_specificity_check()`. Previously the only mechanism was a soft cost penalty (`wt_off_target * count`), which was often insufficient to outweigh thermodynamic scores, allowing off-target-contaminated pairs into the final panel. 5 new tests in `test_blast_specificity.py`.
+- **Configurable BLAST specificity thresholds (`config.py`, `specificity.py`)**: New `BlastParameters` config section with 5 fields that were previously hardcoded: `length_threshold` (default 15), `evalue_threshold` (default 10.0), `max_mismatches` (default 2), `max_amplicon_size` (default 2000), and `ontarget_tolerance` (default 5). These control which BLAST hits are classified as "predicted bound" and how on-target vs off-target amplicons are distinguished. Added to both `designer_default_config.json` and `designer_lenient_config.json` (identical values — physics-based, not stringency-based). `run_specificity_check()` and `_is_on_target()` now accept these as parameters, threaded from `pipeline.py` via `config.blast_parameters`. 9 new tests across `test_config.py` and `test_blast_specificity.py`.
+
 ## [1.0.0b3] - 03-03-2026
 
 ### Fixed

config/designer_default_config.json

@@ -70,6 +70,13 @@
         "snp_strict": false,
         "snp_af_weight": 1.0
     },
+    "blast_parameters": {
+        "length_threshold": 15,
+        "evalue_threshold": 10.0,
+        "max_mismatches": 2,
+        "max_amplicon_size": 2000,
+        "ontarget_tolerance": 5
+    },
     "multiplex_picker_parameters": {
         "initial_solutions": 100,
         "top_solutions_to_keep": 4,

config/designer_lenient_config.json

@@ -70,6 +70,13 @@
         "snp_strict": false,
         "snp_af_weight": 0.5
     },
+    "blast_parameters": {
+        "length_threshold": 15,
+        "evalue_threshold": 10.0,
+        "max_mismatches": 2,
+        "max_amplicon_size": 2000,
+        "ontarget_tolerance": 5
+    },
     "multiplex_picker_parameters": {
         "initial_solutions": 100,
         "top_solutions_to_keep": 4,

docs/CONFIG_GUIDE.md

@@ -0,0 +1,122 @@
+# Configuration Guide for Plexus
+
+Plexus uses a hierarchical JSON configuration to control every aspect of the multiplex PCR primer design process. This guide explains the parameters, their impact on the design, and provides sensible ranges for typical applications.
+
+## Configuration Structure
+
+The configuration is divided into six main sections:
+
+1. `singleplex_design_parameters`: Individual primer properties.
+2. `primer_pair_parameters`: Properties of the forward/reverse pair (e.g., amplicon size).
+3. `pcr_conditions`: Thermodynamic environment (salt, concentrations).
+4. `snp_check_parameters`: How variants (SNPs) affect primer scoring.
+5. `blast_parameters`: Specificity and off-target detection thresholds.
+6. `multiplex_picker_parameters`: Weights for the final panel optimization.
+
+---
+
+## 1. Singleplex Design Parameters
+
+These parameters define what makes a "good" individual primer.
+
+| Parameter | Default | Sensible Range | Description |
+| :--- | :--- | :--- | :--- |
+| `PRIMER_OPT_TM` | 60.0 | 58.0 - 62.0 | The ideal melting temperature (°C). |
+| `PRIMER_MIN_TM` | 57.0 | 55.0 - 59.0 | Minimum allowed Tm. |
+| `PRIMER_MAX_TM` | 63.0 | 61.0 - 68.0 | Maximum allowed Tm. |
+| `PRIMER_OPT_SIZE` | 22 | 18 - 26 | The ideal primer length (bp). |
+| `primer_min_length`| 18 | 15 - 20 | Minimum primer length. |
+| `primer_max_length`| 28 | 25 - 35 | Maximum primer length. |
+| `primer_min_gc` | 30 | 20 - 40 | Minimum GC content (%). |
+| `primer_max_gc` | 70 | 60 - 80 | Maximum GC content (%). |
+| `primer_gc_clamp` | 1 | 0 or 1 | If 1, requires 1-3 G/C bases in the last 5 bases of the 3' end. |
+| `primer_max_poly_x`| 5 | 3 - 6 | Max allowed homopolymer length (e.g., AAAAA). |
+| `PRIMER_MAX_HAIRPIN_TH` | 24.0 | 20.0 - 30.0 | Max allowed Tm (°C) for internal hairpin structures. |
+| `PRIMER_MAX_SELF_ANY_TH` | 45.0 | 35.0 - 50.0 | Max allowed Tm (°C) for any self-dimer. |
+| `PRIMER_MAX_SELF_END_TH` | 35.0 | 25.0 - 40.0 | Max allowed Tm (°C) for 3'-anchored self-dimers. |
+
+---
+
+## 2. Primer Pair Parameters
+
+These parameters control the relationship between the forward and reverse primers.
+
+| Parameter | Default | Sensible Range | Description |
+| :--- | :--- | :--- | :--- |
+| `PRIMER_PAIR_MAX_DIFF_TM` | 3.0 | 1.0 - 5.0 | Max difference in Tm between F and R primers. |
+| `PRIMER_PRODUCT_OPT_SIZE` | 60 | 50 - 150 | Ideal total amplicon length (including primers). |
+| `PRIMER_PRODUCT_MAX_SIZE` | 120 | 80 - 250 | Maximum total amplicon length. |
+| `PRIMER_PRODUCT_MIN_INSERT_SIZE` | 20 | 10 - 50 | Minimum bases between the two primers. |
+
+---
+
+## 3. PCR Conditions
+
+These values are used for thermodynamic calculations (SantaLucia 1998). Ensure these match your actual lab protocol for accurate Tm and Dimer predictions.
+
+| Parameter | Default | Description |
+| :--- | :--- | :--- |
+| `annealing_temperature` | 60.0 | The temperature used for the annealing step (°C). |
+| `mv_concentration` | 50.0 | Monovalent cation concentration (mM), usually KCl. |
+| `dv_concentration` | 1.5 | Divalent cation concentration (mM), usually MgCl2. |
+| `dntp_concentration` | 0.6 | Total dNTP concentration (mM). |
+| `primer_concentration` | 50.0 | Concentration of each individual primer (nM). |
+
+---
+
+## 4. SNP Check Parameters
+
+Plexus queries VCF files to avoid designing primers over common variants.
+
+| Parameter | Default | Sensible Range | Description |
+| :--- | :--- | :--- | :--- |
+| `af_threshold` | 0.01 | 0.001 - 0.05 | Minimum Allele Frequency (AF) to consider a SNP problematic. |
+| `snp_penalty_weight` | 10.0 | 5.0 - 20.0 | Base penalty added to a primer for each overlapping SNP. |
+| `snp_3prime_window` | 5 | 3 - 8 | Bases from the 3' end where a SNP is considered high-impact. |
+| `snp_3prime_multiplier` | 3.0 | 2.0 - 5.0 | Multiplier for the penalty if the SNP is in the 3' window. |
+| `snp_strict` | false | true/false | If true, any primer with a SNP > `af_threshold` is discarded immediately. |
+| `snp_af_weight` | 0.0 | 0.0 - 1.0 | Exponent for AF scaling. `0.5` (sqrt) is recommended to penalize common SNPs more than rare ones. |
+
+---
+
+## 5. BLAST Parameters
+
+Used to identify off-target products that might cause non-specific amplification.
+
+| Parameter | Default | Description |
+| :--- | :--- | :--- |
+| `length_threshold` | 15 | Min 3'-anchored alignment length (bp) to predict binding. |
+| `max_mismatches` | 2 | Max mismatches allowed in a 3'-anchored alignment. |
+| `max_amplicon_size` | 2000 | Max distance between two hits to be considered a potential off-target amplicon. |
+| `ontarget_tolerance` | 5 | BP tolerance when verifying if a BLAST hit matches the intended target. |
+
+---
+
+## 6. Multiplex Picker Parameters
+
+These weights control the final "cost" of a panel. The optimizer tries to minimize this cost.
+
+| Parameter | Default | Description |
+| :--- | :--- | :--- |
+| `wt_pair_penalty` | 1.0 | Weight for the basic primer quality penalty (Tm deviation, length, etc). |
+| `wt_off_target` | 5.0 | Penalty for each predicted off-target product. Usually high to avoid mispriming. |
+| `wt_cross_dimer` | 1.0 | Weight for penalties arising from primers in DIFFERENT pairs forming dimers. |
+| `wt_pair_dimer` | 1.0 | Weight for the dimer score of the F/R primers within the SAME pair. |
+| `wt_snp_penalty` | 3.0 | Weight for the accumulated SNP penalty from the `snpcheck` step. |
+| `initial_solutions`| 100 | Number of iterations for stochastic selectors (Greedy, Simulated Annealing). |
+
+---
+
+## Best Practices
+
+### For Small, High-Quality Panels (< 24 targets)
+
+* Use `preset: "default"`.
+* Set `snp_strict: true` to ensure no variants interfere with your clinical targets.
+* Increase `initial_solutions` to `1000` for a more thorough search.
+
+### For Large Discovery Panels (> 100 targets)
+
+* Use `preset: "lenient"` to allow for a wider range of Tms.
+* Keep `snp_strict: false` but use `snp_af_weight: 0.5` to prioritize the best available sites.
+* Use the `Greedy` or `SimulatedAnnealing` selectors; avoid `BruteForce` or `DFS`.

docs/COORDINATE_SYSTEMS.md

@@ -0,0 +1,72 @@
+# Coordinate Systems in Plexus
+
+For a tool where off-by-one errors are safety-critical, understanding the coordinate systems used for inputs, internal processing, and outputs is essential.
+
+## Summary
+
+| Context | Coordinate System | Base | Inclusive |
+|---------|-------------------|------|-----------|
+| **Input CSV** | Genomic | 1 | Yes |
+| **Output CSV** | Genomic | 1 | Yes |
+| **Internal Offsets** | Design-Region Relative | 0 | Yes |
+| **Reference FASTA** | Standard (pysam handles) | 1 (text) | Yes |
+| **VCF/gnomAD** | Genomic | 1 | Yes |
+
+---
+
+## Genomic Coordinates (1-based)
+
+Plexus follows the standard convention used by genome browsers (IGV, UCSC), VCF files, and GFF/GTF annotations.
+
+### Input: `junctions.csv`
+The columns `Five_Prime_Coordinate` and `Three_Prime_Coordinate` must be provided as **1-based inclusive** genomic coordinates.
+- For a single-base target (e.g., an SNV), both coordinates should be the same.
+- For a range (e.g., a deletion or a whole exon), provide the first and last base of the target.
+
+### Output: `selected_multiplex.csv` / `candidate_pairs.csv`
+All genomic coordinate columns in the output files are **1-based inclusive**:
+- `Forward_Genomic_Start` / `Forward_Genomic_End`
+- `Reverse_Genomic_Start` / `Reverse_Genomic_End`
+- `Junction_Start` / `Junction_End`
+
+This ensures that you can directly paste these coordinates into IGV or a genome browser to verify the primer positions.
+
+---
+
+## Design-Region Offsets (0-based)
+
+Internally, Plexus extracts a padded genomic sequence around your target, known as the **Design Region**.
+
+### `design_start`
+The attribute `junction.design_start` stores the **1-based genomic coordinate** of the first base in the extracted sequence.
+
+### Primer Offsets
+The `start` attribute of a `Primer` object is a **0-based index** into the `design_region` string.
+- A primer starting at the very beginning of the design region has a `start` of `0`.
+- The genomic start of a primer is calculated as: `design_start + primer.start`.
+- The genomic end of a primer is calculated as: `design_start + primer.start + primer.length - 1`.
+
+---
+
+## External Tools and Libraries
+
+Plexus handles the conversion between its internal 1-based system and the 0-based systems used by some underlying libraries:
+
+### pysam (VCF/FASTA access)
+`pysam` uses 0-based half-open coordinates (like the BAM/BED formats).
+- **Extraction**: When Plexus fetches a sequence from the FASTA, it converts the 1-based `[start, end]` range to 0-based `[start-1, end)`.
+- **SNP Checking**: When querying a VCF, Plexus converts the 1-based primer range `[P_start, P_end]` to the 0-based query `fetch(chrom, P_start - 1, P_end)`.
+
+### Primer3
+`primer3-py` uses 0-based coordinates for sequence offsets. Plexus passes the `design_region` string and receives 0-based offsets, which are stored directly in the `Primer` objects.
+
+### BLAST
+The BLAST specificity check uses 1-based coordinates for its output, which Plexus parses and maps back to the original genomic coordinates for on-target/off-target classification.
+
+---
+
+## Developer Note: Coordinate Math
+
+When performing arithmetic on coordinates, always verify the base:
+- **Genomic to Design-Region**: `junction.start - junction.design_start` (both 1-based) yields a 0-based offset.
+- **Design-Region to Genomic**: `junction.design_start + primer.offset` (1-based + 0-based) yields a 1-based genomic coordinate.

docs/USER_GUIDE.md

@@ -11,6 +11,7 @@
 7. [Output Interpretation](#output-interpretation)
 8. [Troubleshooting](#troubleshooting)
 9. [Technical Reference](#technical-reference)
+    - [Coordinate Systems](#coordinate-systems-reference)
 10. [Compliance and Clinical Use](#compliance-and-clinical-use)
     - [Two Deployment Paths](#two-compliance-deployment-paths)
     - [The Compliance Manifest](#the-compliance-manifest)
@@ -147,6 +148,8 @@ plexus run -i junctions.csv -g hg38 -o results/
 
 A **Junction** represents a single genomic target for primer design:
 
+> **Important**: Plexus uses **1-based inclusive genomic coordinates** for all target inputs and final outputs. See the [Coordinate Systems Guide](COORDINATE_SYSTEMS.md) for full details.
+
 ```python
 Junction(
     name="EGFR_T790M",
@@ -766,6 +769,18 @@ plexus status
 
 ## Technical Reference
 
+### Coordinate Systems Reference
+
+For a tool where off-by-one errors are safety-critical, understanding the coordinate systems used for inputs, internal processing, and outputs is essential.
+
+| Context | Coordinate System | Base | Inclusive |
+|---------|-------------------|------|-----------|
+| **Input CSV** | Genomic | 1 | Yes |
+| **Output CSV** | Genomic | 1 | Yes |
+| **Internal Offsets** | Design-Region Relative | 0 | Yes |
+
+For more detailed coordinate system specifications and internal handling, please see the [Coordinate Systems Guide](COORDINATE_SYSTEMS.md).
+
 ### CLI Command Reference
 
 **Main Commands:**

src/plexus/blast/specificity.py

@@ -9,7 +9,16 @@
 
 
 def run_specificity_check(
-    panel: MultiplexPanel, work_dir: str, genome_fasta: str, num_threads: int = 1
+    panel: MultiplexPanel,
+    work_dir: str,
+    genome_fasta: str,
+    num_threads: int = 1,
+    *,
+    length_threshold: int = 15,
+    evalue_threshold: float = 10.0,
+    max_mismatches: int = 2,
+    max_amplicon_size: int = 2000,
+    ontarget_tolerance: int = 5,
 ):
     """
     Run BLAST on all candidate primers in the panel to check for specificity
@@ -21,6 +30,12 @@ def run_specificity_check(
         panel: The MultiplexPanel object containing junctions and primer designs.
         work_dir: Directory to store temporary BLAST files.
         genome_fasta: Path to the reference genome FASTA file.
+        num_threads: Number of BLAST threads.
+        length_threshold: Minimum 3'-anchored alignment length (bp) to predict binding.
+        evalue_threshold: E-value cutoff for predicted binding.
+        max_mismatches: Maximum mismatches in a 3'-anchored alignment.
+        max_amplicon_size: Maximum amplicon size (bp) to consider.
+        ontarget_tolerance: Coordinate tolerance (bp) for on-target classification.
     """
     logger.info("Starting specificity check (BLAST)...")
     os.makedirs(work_dir, exist_ok=True)
@@ -43,8 +58,13 @@ def run_specificity_check(
 
     runner = BlastRunner(input_fasta, genome_fasta)
     runner.create_database()
-    runner.run(output_archive=blast_archive, num_threads=num_threads)
-    runner.reformat_output_as_table(blast_table)
+    try:
+        runner.run(output_archive=blast_archive, num_threads=num_threads)
+        runner.reformat_output_as_table(blast_table)
+    finally:
+        if os.path.exists(blast_archive):
+            os.remove(blast_archive)
+            logger.debug(f"Removed BLAST archive: {blast_archive}")
 
     blast_df = runner.get_dataframe()
 
@@ -56,14 +76,16 @@ def run_specificity_check(
     target_map = getattr(panel, "primer_target_map", {})
     annotator = BlastResultsAnnotator(blast_df, target_map=target_map)
     annotator.build_annotation_dict(
-        length_threshold=15, evalue_threshold=10, max_mismatches=2
+        length_threshold=length_threshold,
+        evalue_threshold=evalue_threshold,
+        max_mismatches=max_mismatches,
     )
     annotator.add_annotations()
 
     # 5. Find Off-Target Amplicons
     bound_df = annotator.get_predicted_bound()
     finder = AmpliconFinder(bound_df, target_map=target_map)
-    finder.find_amplicons(max_size_bp=2000)
+    finder.find_amplicons(max_size_bp=max_amplicon_size)
 
     all_amplicons_df = finder.amplicon_df
 
@@ -106,7 +128,7 @@ def run_specificity_check(
             off_targets = []
             on_targets = []
             for prod in potential_products:
-                if _is_on_target(prod, junction, pair):
+                if _is_on_target(prod, junction, pair, tolerance=ontarget_tolerance):
                     on_targets.append(prod)
                 else:
                     off_targets.append(prod)
@@ -191,7 +213,7 @@ def filter_offtarget_pairs(panel: MultiplexPanel) -> tuple[int, list[str]]:
     return total_removed, fallback_junctions
 
 
-def _is_on_target(prod: dict, junction, pair) -> bool:
+def _is_on_target(prod: dict, junction, pair, tolerance: int = 5) -> bool:
     """Check if a BLAST amplicon overlaps the intended target region.
 
     Compares the BLAST hit genomic coordinates against the expected
@@ -224,9 +246,8 @@ def _is_on_target(prod: dict, junction, pair) -> bool:
     expected_rev_start = design_start + pair.reverse.start + pair.reverse.length - 1
 
     # Allow small tolerance for BLAST coordinate alignment differences.
-    # 5 bp chosen to absorb minor alignment shifts while still distinguishing
+    # Default 5 bp absorbs minor alignment shifts while still distinguishing
     # on-target hits from nearby off-target loci.
-    tolerance = 5  # bp
     fwd_match = abs(prod["F_start"] - expected_fwd_start) <= tolerance
     rev_match = abs(prod["R_start"] - expected_rev_start) <= tolerance