Skills › Scientific Agent Skills › primer-design › screen_specificity.py
screen_specificity.py
180 lines · skills/primer-design/scripts/screen_specificity.py at 5f2bc03
1#!/usr/bin/env python3
2"""Screen primers against an explicit local FASTA with reproducible search bounds."""
4from __future__ import annotations
6import argparse
7import platform
8import sys
9from datetime import datetime, timezone
10from pathlib import Path
12sys.dont_write_bytecode = True
14from _common import check_output_paths, read_fasta, read_pairs, sha256_file, write_json
15from _specificity import blast_hits, enumerate_products, exhaustive_hits, expand_panel_pairs, read_expected
18def build_parser():
19 parser = argparse.ArgumentParser(description=__doc__)
20 parser.add_argument("--pairs", required=True, type=Path, help="TSV: pair_id, forward, reverse; optional forward_tail, reverse_tail")
21 parser.add_argument("--reference", required=True, type=Path, help="FASTA with unique record IDs")
22 parser.add_argument("--output", required=True, type=Path, help="JSON report; incomplete runs exit 2")
23 parser.add_argument("--expected", type=Path, help="TSV: pair_id, record_id, start, end (0-based half-open)")
24 parser.add_argument("--engine", choices=("exhaustive", "blast"), default="exhaustive")
25 parser.add_argument("--multiplex", action="store_true", help="Also screen every cross-pair oligo combination; all cross-pair products are potential off-targets")
26 parser.add_argument("--max-panel-combinations", type=int, default=10_000, help="Cross-pair oligo combination cap, checked before panel expansion")
27 parser.add_argument("--circular", action="append", default=[], metavar="RECORD_ID", help="Repeat for circular records; intervals can end beyond record length")
28 parser.add_argument("--min-product", type=int, default=40)
29 parser.add_argument("--max-product", type=int, default=2000)
30 parser.add_argument("--max-mismatches", type=int, default=2)
31 parser.add_argument("--three-prime-bases", type=int, default=5)
32 parser.add_argument("--max-three-prime-mismatches", type=int, default=0)
33 parser.add_argument("--max-comparisons", type=int, default=100_000_000, help="Base comparison cap; exceeding a cap fails closed")
34 parser.add_argument("--max-hits", type=int, default=100_000)
35 parser.add_argument("--max-products", type=int, default=100_000, help="Per-pair product cap")
36 parser.add_argument("--max-product-combinations", type=int, default=1_000_000, help="Per-pair inward-facing candidate comparison cap")
37 parser.add_argument("--max-blast-rows", type=int, default=1_000_000)
38 parser.add_argument("--blastn", default="blastn", help="Local BLAST+ blastn executable")
39 parser.add_argument("--makeblastdb", default="makeblastdb", help="Local BLAST+ makeblastdb executable")
40 parser.add_argument("--blast-max-target-seqs", type=int, default=100_000)
41 parser.add_argument("--blast-max-hsps", type=int, default=100_000)
42 parser.add_argument("--blast-timeout", type=int, default=300, help="Seconds per BLAST+ invocation")
43 return parser
46def run(args):
47 settings = {key: value for key, value in vars(args).items()
48 if key not in {"pairs", "reference", "expected", "output"}}
49 report = {
50 "schema_version": "1.0", "analysis": "local_primer_specificity", "status": "incomplete",
51 "complete_within_model": False, "exhaustive": False,
52 "created_utc": datetime.now(timezone.utc).isoformat(),
53 "coordinate_convention": "0-based half-open forward-reference; circular end may exceed record length",
54 "strand_convention": "+ matches reference sequence and extends right; - matches its reverse complement and extends left",
55 "settings": settings, "provenance": {"python": platform.python_version()},
56 "limitations": [
57 "Potential amplification under a bounded ungapped mismatch model; not experimental specificity.",
58 "Only supplied reference records, mismatch thresholds, 3-prime window and product-size range are screened.",
59 "Unrepresented variants, contaminants, insertions/deletions, and chemistry are not screened.",
60 "Reference ambiguity is conservatively compatible when a primer base is possible; mismatches are lower bounds.",
61 "Only annealing cores are searched. Tails affect final length but tail-created later-cycle binding is not searched.",
62 "F/F, R/R and both F/R orientations are enumerated; overlapping binding sites are excluded.",
63 ],
64 "issues": [], "pairs": [],
65 }
66 report["limitations"].append(
67 "Multiplex mode treats all cross-pair products as unintended; predeclared intentional cross-pair "
68 "designs should be assessed as explicit standalone pairs. Same-oligo products remain in original-pair results."
69 if args.multiplex else "Cross-pair multiplex products are not screened; enable --multiplex for a panel."
70 )
71 try:
72 input_paths = [args.pairs, args.reference] + ([args.expected] if args.expected else [])
73 check_output_paths(input_paths, [args.output])
74 for key in ("min_product", "max_product", "three_prime_bases", "max_comparisons", "max_hits",
75 "max_products", "max_product_combinations", "max_blast_rows", "blast_max_target_seqs",
76 "blast_max_hsps", "blast_timeout", "max_panel_combinations"):
77 if settings[key] <= 0:
78 raise ValueError(f"{key} must be positive")
79 if args.min_product > args.max_product:
80 raise ValueError("min_product must not exceed max_product")
81 if args.max_mismatches < 0 or args.max_three_prime_mismatches < 0:
82 raise ValueError("Mismatch thresholds must be nonnegative")
83 if args.max_three_prime_mismatches > args.three_prime_bases:
84 raise ValueError("3-prime mismatch threshold exceeds its window")
85 pairs, references = read_pairs(args.pairs), read_fasta(args.reference)
86 circular = set(args.circular)
87 if circular - set(references):
88 raise ValueError("Unknown circular record IDs: " + ", ".join(sorted(circular - set(references))))
89 primers = sorted({pair[key] for pair in pairs for key in ("forward", "reverse")})
90 if not primers or not references:
91 raise ValueError("Nonempty primer pairs and reference records are required")
92 if any(args.three_prime_bases > len(primer) for primer in primers):
93 raise ValueError("3-prime window exceeds a primer's annealing-core length")
94 if any(set(primer) - set("ACGT") for primer in primers):
95 raise ValueError("Only unambiguous ACGT annealing cores are supported")
96 report["provenance"].update({
97 "implementation_sha256": {
98 name: sha256_file(Path(__file__).parent / name)
99 for name in ("screen_specificity.py", "_specificity.py", "_common.py")
100 },
101 "pairs": {"path": str(args.pairs), "sha256": sha256_file(args.pairs)},
102 "reference": {"path": str(args.reference), "sha256": sha256_file(args.reference),
103 "record_count": len(references), "total_bases": sum(map(len, references.values())),
104 "record_lengths": {key: len(value) for key, value in references.items()},
105 "ambiguous_base_count": sum(base not in "ACGT" for sequence in references.values() for base in sequence)},
106 "expected": {"path": str(args.expected), "sha256": sha256_file(args.expected)} if args.expected else None,
107 })
108 expected = read_expected(args.expected, pairs, references, circular)
109 original_pair_count = len(pairs)
110 if args.multiplex:
111 pairs = expand_panel_pairs(pairs, args.max_panel_combinations)
112 report["provenance"]["panel"] = {
113 "enabled": args.multiplex, "original_pair_count": original_pair_count,
114 "source_oligo_count": original_pair_count * 2,
115 "cross_pair_count": len(pairs) - original_pair_count,
116 "unique_annealing_cores_searched": len(primers),
117 "generated_pair_ids": [pair["pair_id"] for pair in pairs if pair.get("cross_pair")],
118 }
119 search = exhaustive_hits if args.engine == "exhaustive" else blast_hits
120 if args.engine == "blast":
121 report["limitations"].append(
122 "BLAST discovery is heuristic, uses exact 7-base seeds and E-value 1000; full-length ungapped "
123 "realignment validates only discovered sites. No-hit output never proves absence of binding sites."
124 )
125 hits, engine_provenance, issues = search(primers, references, circular, settings)
126 report["provenance"].update(engine_provenance)
127 report["issues"].extend(issues)
128 for pair in pairs:
129 result = enumerate_products(pair, hits, references, circular, settings, expected.get(pair["pair_id"], []))
130 if result["unresolved_ambiguity"]:
131 report["issues"].append(
132 f"Pair {pair['pair_id']}: ambiguous reference bases at potential binding sites; "
133 "binding and 3-prime mismatch counts remain unresolved"
134 )
135 if issues:
136 result.setdefault("status_before_incomplete", result["status"])
137 result["status"] = "incomplete"
138 report["pairs"].append(result)
139 report["complete_within_model"] = args.engine == "exhaustive" and not report["issues"]
140 report["exhaustive"] = args.engine == "exhaustive" and not report["issues"]
141 report["discovery_completed"] = not issues
142 statuses = {pair["status"] for pair in report["pairs"]}
143 report["status"] = next(status for status in (
144 "incomplete", "intended_target_not_found", "potential_off_target", "no_expected_target",
145 "no_off_target_found_within_search_scope",
146 ) if status in statuses)
147 return report, 2 if report["issues"] else 0
148 except (ValueError, OSError, KeyError) as exc:
149 report["issues"].append(str(exc))
150 report["status"] = "incomplete"
151 report["complete_within_model"] = False
152 report["exhaustive"] = False
153 report["discovery_completed"] = False
154 for pair in report["pairs"]:
155 pair.setdefault("status_before_incomplete", pair["status"])
156 pair["status"] = "incomplete"
157 return report, 2
160def main(argv=None):
161 args = build_parser().parse_args(argv)
162 # Never replace an input with an error report, either.
163 inputs = [args.pairs, args.reference] + ([args.expected] if args.expected else [])
164 try:
165 check_output_paths(inputs, [args.output])
166 except (ValueError, OSError) as exc:
167 print(f"error: {exc}", file=sys.stderr)
168 return 2
169 report, exit_code = run(args)
170 try:
171 write_json(args.output, report)
172 except OSError as exc:
173 print(f"Cannot write report: {exc}", file=sys.stderr)
174 return 2
175 print(f"{report['status']}: {args.output}")
176 return exit_code
179if __name__ == "__main__":
180 raise SystemExit(main())