Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension


Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion .pre-commit-config.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,7 @@ repos:
hooks:
- id: pytest
name: pytest
entry: .venv/bin/python -m pytest -m "not integration" --no-header -q
entry: .venv/bin/python -m pytest -m "not integration" --no-header -q --no-cov
language: system
pass_filenames: false
always_run: true
Expand Down
14 changes: 14 additions & 0 deletions CHANGELOG.md
Original file line number Diff line number Diff line change
Expand Up @@ -5,13 +5,27 @@ 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.1.0] - 13-03-2026

### Changed

- **Vectorize off-target amplicon finding** (`blast/offtarget_finder.py`): Replace O(F × R) nested `iterrows()` + `query()` loop with sorted arrays and `np.searchsorted` for O(F log R) per-chromosome lookups. Eliminates heavy per-row Pandas overhead for large multiplex panels.
- **Vectorize BLAST annotation** (`blast/annotator.py`): Replace 4× `apply(axis=1)` Python lambdas with direct vectorized column operations for `from_3prime`, `length_pass_3prime`, `evalue_pass_3prime`, and `predicted_bound`.

### Added

- **HTML Visual QC Report (REPT-02)**: New self-contained HTML report (`panel_report.html`) generated alongside `panel_qc.json` after each pipeline run. Includes interactive Plotly.js charts for Tm distribution, amplicon size, GC content, and a cross-reactivity heatmap, plus tables for sequence flags, off-targets, failed junctions, and solution comparison. Reports inline Plotly.js (~1MB) for fully offline viewing. Added `plexus report` CLI command for standalone re-generation from existing pipeline output.
- **Jinja2 dependency** added for HTML template rendering.

## [1.0.2] - 04-03-2026

### Changed

- **Rich progress bars replace log output during pipeline runs**: The CLI now shows clean Rich progress bars on stderr (step-level + per-junction detail for primer design and SNP check) instead of a wall of log messages. All detailed logs still go to the file. Warnings and errors are printed above the progress bar. Multi-panel parallel mode shows a panel-level progress bar. Progress bars are only active when stderr is a TTY; non-interactive runs behave as before.
- **SNP and off-target filters now retain all tied least-affected pairs** (`snpcheck/checker.py`, `blast/specificity.py`): When all primer pairs for a junction overlap SNPs or have off-target products, the filters now keep every pair tied at the minimum count instead of arbitrarily picking one. This lets the downstream selector evaluate tied candidates on other properties (Tm, GC%, pair penalty, etc.).
- **Normalise cross-dimer penalty in multiplex cost function** (`selector/cost.py`): The cross-dimer penalty was a raw sum over all C(2n, 2) pairwise primer interactions, scaling quadratically with multiplex size. This caused it to dominate the cost function at higher plexities, effectively drowning out off-target and SNP penalties during selection. The penalty is now divided by the number of interactions, making it a per-interaction average. Weights are now directly comparable regardless of multiplex size.
- **Increase default cross-dimer weight** (`data/designer_default_config.json`): `wt_cross_dimer` 1.0 → 20.0 to better penalise adapter-driven cross-dimer formation in multiplex selection, particularly for panels with custom adapter tails.
- **Off-target CSV now includes alignment details** (`blast/offtarget_finder.py`, `designer/multiplexpanel.py`): `off_targets.csv` now reports per-primer alignment quality for each off-target product: percent identity, number of mismatches, alignment length, and E-value for both forward and reverse primer hits. Helps assess the likelihood of off-target amplification.

- **Separate warnings from errors in pipeline output** (`pipeline.py`, `cli.py`): Off-target and SNP fallback messages (where all pairs had issues but the least-affected were kept) were incorrectly reported as errors, causing the CLI to display "Some panels had errors" for panels that completed successfully. These are now reported as warnings. Errors are reserved for actual failures (e.g. design exceptions, BLAST unavailable). The CLI now shows a distinct warnings section below the success summary.
- **Update fallback message wording**: "least-affected pair kept" now reads "all least-affected pairs kept" to reflect the v1.0.2 change that retains all tied pairs.
Expand Down
65 changes: 60 additions & 5 deletions docs/ROADMAP.md
Original file line number Diff line number Diff line change
Expand Up @@ -589,6 +589,55 @@ the BLAST database may be incomplete.

---

### SPEC-01 · Assign cross-pair off-target amplicons to primer pairs

**Severity: Important · Files: `src/plexus/blast/specificity.py`, `src/plexus/selector/cost.py`**

`AmpliconFinder.find_amplicons()` correctly discovers amplicons formed by primers from different
pairs (e.g., forward from pair A + reverse from pair B on the same chromosome). However, the
mapping step in `run_specificity_check()` (step 6, line 112+) only looks up amplicons keyed by
a single pair's own `(f_id, r_id)`. Cross-pair amplicons exist in `amplicon_map` but are never
retrieved — they silently fall through because no single `PrimerPair` owns both primer IDs.

**Impact:** For multiplex panels, cross-pair off-target products are the most likely source of
spurious amplification (many more primer combinations than same-pair). These products are
computed by BLAST and identified by `AmpliconFinder`, but never assigned to any pair. They do
not appear in `off_targets.csv`, do not affect `filter_offtarget_pairs()`, and do not contribute
to the `wt_off_target` cost term. The specificity check is therefore blind to the dominant class
of off-target risk in multiplex reactions.

**Required changes:**

1. After the per-pair mapping loop in `run_specificity_check()`, perform a second pass that
identifies cross-pair amplicons — entries in `amplicon_map` where `F_primer` belongs to one
pair and `R_primer` belongs to a different pair (or different junction).
2. Assign each cross-pair amplicon to **both** contributing pairs (the pair that owns the forward
primer and the pair that owns the reverse primer), so that both are penalised. Add a
`cross_pair: bool` flag to the product dict to distinguish these from same-pair off-targets.
3. Ensure cross-pair off-targets are surfaced in `off_targets.csv` with both contributing pair
IDs visible.
4. Verify that `filter_offtarget_pairs()` and `wt_off_target` in the cost function correctly
account for the increased off-target counts.

**Design consideration:** A cross-pair amplicon penalises two pairs simultaneously. If pair A's
forward primer participates in cross-products with 5 different reverse primers from other pairs,
replacing pair A eliminates all 5 — but the current per-pair count treats them independently.
Consider whether cross-pair off-targets should carry a reduced weight or whether the existing
`wt_off_target` is sufficient given that the selector can swap individual pairs.

**Relationship to other items:**
- ISPCR-04 (cross-target interaction matrix) is a reporting view of the same underlying data;
SPEC-01 makes the data actionable in the selector.
- ISPCR-02 (ΔG-weighted off-target cost) would naturally extend to cross-pair products once
they are tracked.

**Tests to add:** `tests/test_specificity.py` — construct a mock `bound_df` with primers from
two different pairs that form a cross-pair amplicon; assert that both pairs receive the
off-target product. Assert that the off-target CSV includes the cross-pair product with both
pair IDs.

---

### ISPCR-04 · Improve AmpliconFinder: cross-target interaction matrix
**File: `src/plexus/blast/offtarget_finder.py`**

Expand Down Expand Up @@ -660,13 +709,18 @@ consistently to junction coordinates, VCF queries, and BLAST sequence IDs.

---

### REPT-02 · Visual QC Report (HTML)
### ~~REPT-02 · Visual QC Report (HTML)~~ ✅ Implemented in v1.1.0

**Severity: Low**

Transform the `panel_qc.json` (REPT-01) into a visual, standalone HTML report. This should
include Plotly or Seaborn charts for distributions and a searchable heatmap for the
cross-reactivity matrix, facilitating rapid review by clinical lab staff.
Transform the `panel_qc.json` (REPT-01) into a visual, standalone HTML report with Plotly.js
charts for Tm distribution, amplicon size, GC content, and a cross-reactivity heatmap.

**Implementation:** Added `src/plexus/reporting/html_report.py` with `generate_html_report()`
and a Jinja2 template (`panel_report.html.j2`) that inlines Plotly.js basic for fully
offline-capable reports. Integrated into the pipeline after `panel_qc.json` write. Added
`plexus report` CLI command for standalone re-generation. 11 tests in
`tests/test_reporting_html.py`.

---

Expand Down Expand Up @@ -724,11 +778,12 @@ project.
| ISPCR-01 | ntthal ΔG scoring for BLAST binding sites | v1.1 | Important | |
| ISPCR-02 | ΔG-weighted off-target cost in selector | v1.1 | Important | |
| ISPCR-03 | Template mispriming check | v1.1 | Low | |
| SPEC-01 | Assign cross-pair off-target amplicons to primer pairs | v1.1 | Important | |
| ISPCR-04 | Cross-target interaction matrix output | v1.1 | Low | |
| PERF-01 | Parallel thermodynamics by default | v1.1 | Low | |
| PERF-02 | Pre-filter candidates per junction before optimisation | v1.1 | Low | |
| EXT-01 | Additional genome presets | v1.1 | Low | |
| EXT-02 | Chromosome naming normalisation | v1.1 | Low | |
| REPT-02 | Visual QC Report (HTML) | v1.1 | Low | |
| ~~REPT-02~~ | ~~Visual QC Report (HTML)~~ | ~~v1.1~~ | ~~Low~~ | ✅ v1.1.0 |
| SPLIT-01 | Automated Panel Splitting | Future | Future | |
| ~~TEST-01~~ | ~~End-to-end integration test with real BLAST~~ | ~~v1.0~~ | ~~Important~~ | ✅ v0.4.0 |
1 change: 1 addition & 0 deletions pyproject.toml
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,7 @@ classifiers = [
"Programming Language :: Python :: 3.13",
]
dependencies = [
"Jinja2",
"loguru",
"pandas",
"primer3-py>=2.3.0",
Expand Down
45 changes: 25 additions & 20 deletions src/plexus/blast/annotator.py
Original file line number Diff line number Diff line change
Expand Up @@ -38,33 +38,38 @@ def build_annotation_dict(
at or near the 3' end will have ``qend < qlen``. A tolerance
of 2-3 catches these cases without accepting 5'-only hits.
"""
self.annotation_params = {
"length_threshold": length_threshold,
"evalue_threshold": evalue_threshold,
"max_mismatches": max_mismatches,
"three_prime_tolerance": three_prime_tolerance,
}
# Keep annotation names for summarise_by_primer
self.annotations = {
"from_3prime": lambda row: (row["qlen"] - row["qend"])
<= three_prime_tolerance,
"length_pass_3prime": lambda row: (
row["length"] >= length_threshold
and row["mismatch"] <= max_mismatches
and row["from_3prime"]
),
"evalue_pass_3prime": lambda row: (
row["evalue"] < evalue_threshold and row["from_3prime"]
),
"predicted_bound": lambda row: (
row["length_pass_3prime"] or row["evalue_pass_3prime"]
),
"from_3prime": None,
"length_pass_3prime": None,
"evalue_pass_3prime": None,
"predicted_bound": None,
}

def add_annotations(self):
"""
Add annotations to `blast_df`
Add annotations to `blast_df` using vectorized column operations.

"""
for annot_name, annot_func in self.annotations.items():
self.blast_df.insert(
self.blast_df.shape[1],
annot_name,
self.blast_df.apply(func=annot_func, axis=1),
)
p = self.annotation_params
df = self.blast_df

df["from_3prime"] = (df["qlen"] - df["qend"]) <= p["three_prime_tolerance"]
df["length_pass_3prime"] = (
(df["length"] >= p["length_threshold"])
& (df["mismatch"] <= p["max_mismatches"])
& df["from_3prime"]
)
df["evalue_pass_3prime"] = (df["evalue"] < p["evalue_threshold"]) & df[
"from_3prime"
]
df["predicted_bound"] = df["length_pass_3prime"] | df["evalue_pass_3prime"]

def get_predicted_bound(self):
"""
Expand Down
175 changes: 124 additions & 51 deletions src/plexus/blast/offtarget_finder.py
Original file line number Diff line number Diff line change
@@ -1,5 +1,4 @@
from collections import namedtuple

import numpy as np
import pandas as pd


Expand Down Expand Up @@ -30,6 +29,14 @@ class AmpliconFinder:
- F_start: 5' genomic position of the forward primer
- R_start: 5' genomic position of the reverse primer (on minus strand)
- product_bp: Predicted size of the amplicon in base pairs (inclusive)
- F_pident: Percent identity of the forward primer alignment
- R_pident: Percent identity of the reverse primer alignment
- F_mismatch: Number of mismatches in the forward primer alignment
- R_mismatch: Number of mismatches in the reverse primer alignment
- F_align_len: Alignment length of the forward primer (bp)
- R_align_len: Alignment length of the reverse primer (bp)
- F_evalue: E-value of the forward primer alignment
- R_evalue: E-value of the reverse primer alignment
"""

def __init__(self, bound_df, target_map=None):
Expand All @@ -48,63 +55,129 @@ def __init__(self, bound_df, target_map=None):
# Generated
self.amplicon_df = None

def _get_col_or_none(self, df, col):
"""Return column values as array, or None if column is missing."""
if col in df.columns:
return df[col].values
return None

def find_amplicons(self, max_size_bp=6000):
"""
Identify all potential amplicons within max_size_bp.

Uses sorted arrays and np.searchsorted for O(F log R) per chromosome
instead of O(F × R) nested iteration.

Args:
max_size_bp (int): Maximum predicted product size to consider.
"""

# Storage -- can make flexible
Amplicon = namedtuple(
"Amplicon",
[
"chrom",
"F_target",
"R_target",
"F_primer",
"R_primer",
"F_start",
"R_start",
"product_bp",
],
)

# Iterate and find amplicons
amplicons = []
amplicon_columns = [
"chrom",
"F_target",
"R_target",
"F_primer",
"R_primer",
"F_start",
"R_start",
"product_bp",
"F_pident",
"R_pident",
"F_mismatch",
"R_mismatch",
"F_align_len",
"R_align_len",
"F_evalue",
"R_evalue",
]
result_chunks = []
target_map = self.target_map

for chrom, chrom_df in self.bound_df.groupby("sseqid"):
for _, F_df in chrom_df.query("sstrand == 'plus'").iterrows():
# Get start position of the forward primer
# 5' position
F_start = int(F_df["sstart"])

# Get proximal reverse primers, if they exist
# - Very important to get directionality of query correct
R_df = chrom_df.query(
"sstrand == 'minus' and 0 < (sstart - @F_start) < @max_size_bp"
)

# Check for any pairs
if R_df.shape[0] > 0:
F_amplicons = [
Amplicon(
chrom=chrom,
F_target=self.target_map.get(
F_df["qseqid"], F_df["qseqid"].split("_")[0]
),
R_target=self.target_map.get(
row["qseqid"], row["qseqid"].split("_")[0]
),
F_primer=F_df["qseqid"],
R_primer=row["qseqid"],
F_start=F_start,
R_start=row["sstart"],
product_bp=row["sstart"] - F_start + 1,
)
for _, row in R_df.iterrows()
]
amplicons.extend(F_amplicons)
fwd = chrom_df[chrom_df["sstrand"] == "plus"]
rev = chrom_df[chrom_df["sstrand"] == "minus"]

if fwd.empty or rev.empty:
continue

# Sort reverse hits by sstart for searchsorted
rev = rev.sort_values("sstart")
r_starts = rev["sstart"].values
r_qseqids = rev["qseqid"].values
r_pident = self._get_col_or_none(rev, "pident")
r_mismatch = self._get_col_or_none(rev, "mismatch")
r_length = self._get_col_or_none(rev, "length")
r_evalue = self._get_col_or_none(rev, "evalue")

f_starts = fwd["sstart"].values
f_qseqids = fwd["qseqid"].values
f_pident = self._get_col_or_none(fwd, "pident")
f_mismatch = self._get_col_or_none(fwd, "mismatch")
f_length = self._get_col_or_none(fwd, "length")
f_evalue = self._get_col_or_none(fwd, "evalue")

for i in range(len(f_starts)):
f_start = f_starts[i]
# Find reverse hits where 0 < (r_start - f_start) < max_size_bp
lo = np.searchsorted(
r_starts, f_start, side="right"
) # r_start > f_start
hi = np.searchsorted(
r_starts, f_start + max_size_bp, side="left"
) # r_start < f_start + max_size_bp

if lo >= hi:
continue

n_matches = hi - lo
f_qseqid = f_qseqids[i]
f_target = target_map.get(f_qseqid, f_qseqid.split("_")[0])

matched_r_starts = r_starts[lo:hi]
matched_r_qseqids = r_qseqids[lo:hi]

r_targets = [
target_map.get(rq, rq.split("_")[0]) for rq in matched_r_qseqids
]

chunk = {
"chrom": np.full(n_matches, chrom),
"F_target": np.full(n_matches, f_target),
"R_target": r_targets,
"F_primer": np.full(n_matches, f_qseqid),
"R_primer": matched_r_qseqids,
"F_start": np.full(n_matches, f_start, dtype=int),
"R_start": matched_r_starts,
"product_bp": matched_r_starts - f_start + 1,
"F_pident": np.full(
n_matches, f_pident[i] if f_pident is not None else None
),
"R_pident": r_pident[lo:hi]
if r_pident is not None
else np.full(n_matches, None),
"F_mismatch": np.full(
n_matches, f_mismatch[i] if f_mismatch is not None else None
),
"R_mismatch": r_mismatch[lo:hi]
if r_mismatch is not None
else np.full(n_matches, None),
"F_align_len": np.full(
n_matches, f_length[i] if f_length is not None else None
),
"R_align_len": r_length[lo:hi]
if r_length is not None
else np.full(n_matches, None),
"F_evalue": np.full(
n_matches, f_evalue[i] if f_evalue is not None else None
),
"R_evalue": r_evalue[lo:hi]
if r_evalue is not None
else np.full(n_matches, None),
}
result_chunks.append(pd.DataFrame(chunk, columns=amplicon_columns))

# Store
self.amplicon_df = pd.DataFrame(amplicons)
if result_chunks:
self.amplicon_df = pd.concat(result_chunks, ignore_index=True)
else:
self.amplicon_df = pd.DataFrame(columns=amplicon_columns)
Loading
Loading