_specificity.py · Scientific Agent Skills

1"""Bounded, ungapped primer binding and inward-facing amplicon enumeration.

2

3Coordinates are zero-based, half-open on the reference's forward strand.

4For circular records, start is canonical and end may exceed the record length.

5Neither this mismatch model nor BLAST models PCR chemistry or proves specificity.

6"""

7

8from __future__ import annotations

9

10import csv

11import hashlib

12import shutil

13import subprocess

14import tempfile

15from collections import Counter, defaultdict

16from itertools import combinations

17from pathlib import Path

18

19from _common import revcomp

20

21

22IUPAC = {

23 "A": "A", "C": "C", "G": "G", "T": "T", "R": "AG", "Y": "CT",

24 "S": "GC", "W": "AT", "K": "GT", "M": "AC", "B": "CGT", "D": "AGT",

25 "H": "ACT", "V": "ACG", "N": "ACGT",

26}

27

28

29class SearchIncomplete(ValueError):

30 """The requested screen cannot support a completed-search result."""

31

32

33def source_oligos(pair):

34 """Retain the physical oligo identity behind each local F/R role."""

35 return pair.get("source_oligos") or {

36 role: {

37 "oligo_id": f"{pair['pair_id']}:{role}", "pair_id": pair["pair_id"],

38 "primer_role": role, "annealing_sequence": pair[key],

39 "tail_5prime": pair.get(key + "_tail", ""),

40 }

41 for role, key in (("F", "forward"), ("R", "reverse"))

42 }

43

44

45def expand_panel_pairs(pairs, max_combinations):

46 """Add all cross-pair oligo combinations without mutating original pairs."""

47 count = 2 * len(pairs) * (len(pairs) - 1)

48 if count > max_combinations:

49 raise SearchIncomplete(

50 f"Multiplex requires {count:,} cross-oligo combinations; cap is {max_combinations:,}. "

51 "No panel expansion or truncated screen was accepted."

52 )

53 prefix = "multiplex-cross-"

54 while any(pair["pair_id"].startswith(prefix) for pair in pairs):

55 prefix = "_" + prefix

56 generated = []

57 for first, second in combinations(sorted(pairs, key=lambda pair: pair["pair_id"]), 2):

58 for first_role in ("F", "R"):

59 for second_role in ("F", "R"):

60 first_oligo, second_oligo = source_oligos(first)[first_role], source_oligos(second)[second_role]

61 generated.append({

62 "pair_id": f"{prefix}{len(generated) + 1:06d}",

63 "forward": first_oligo["annealing_sequence"],

64 "reverse": second_oligo["annealing_sequence"],

65 "forward_tail": first_oligo["tail_5prime"],

66 "reverse_tail": second_oligo["tail_5prime"],

67 "cross_pair": True, "source_oligos": {"F": first_oligo, "R": second_oligo},

68 })

69 return list(pairs) + generated

70

71

72def read_expected(path, pairs, references, circular):

73 """Read exact intended product intervals; reject unverifiable identifiers."""

74 result = {pair["pair_id"]: [] for pair in pairs}

75 if path is None:

76 return result

77 with Path(path).open(newline="", encoding="utf-8-sig") as handle:

78 reader = csv.DictReader(handle, delimiter="\t")

79 required = {"pair_id", "record_id", "start", "end"}

80 header = reader.fieldnames or []

81 if set(header) != required or len(header) != len(required):

82 raise ValueError("Expected-target TSV requires exactly pair_id, record_id, start, end")

83 for line, row in enumerate(reader, 2):

84 if None in row or any(value is None for value in row.values()):

85 raise ValueError(f"Wrong expected-target column count, line {line}")

86 pair_id, record_id = row["pair_id"], row["record_id"]

87 if pair_id not in result or record_id not in references:

88 raise ValueError(f"Unknown pair_id or record_id in expected targets, line {line}")

89 try:

90 start, end = int(row["start"]), int(row["end"])

91 except (ValueError, TypeError) as exc:

92 raise ValueError(f"Non-integer expected coordinates, line {line}") from exc

93 length = len(references[record_id])

94 if not (0 <= start < length and start < end):

95 raise ValueError(f"Invalid expected interval, line {line}")

96 if end > (start + length if record_id in circular else length):

97 raise ValueError(f"Expected interval exceeds molecule, line {line}")

98 interval = {"record_id": record_id, "start": start, "end": end}

99 if interval in result[pair_id]:

100 raise ValueError(f"Duplicate expected interval, line {line}")

101 result[pair_id].append(interval)

102 return result

103

104

105def reference_window(sequence, start, length, circular=False):

106 if circular:

107 if length > len(sequence):

108 return None

109 start %= len(sequence)

110 end = start + length

111 return sequence[start:end] if end <= len(sequence) else sequence[start:] + sequence[:end - len(sequence)]

112 if start < 0 or start + length > len(sequence):

113 return None

114 return sequence[start:start + length]

115

116

117def score_window(primer, window, strand, settings):

118 """Conservative mismatch lower bounds; ambiguity never certifies a base."""

119 oriented = window if strand == "+" else revcomp(window)

120 mismatches, ambiguous = [], []

121 for index, (base, reference_base) in enumerate(zip(primer, oriented), 1):

122 if reference_base not in "ACGT":

123 ambiguous.append(index)

124 if base not in IUPAC[reference_base]:

125 mismatches.append(index)

126 boundary = len(primer) - settings["three_prime_bases"]

127 terminal = [index for index in mismatches if index > boundary]

128 if (len(mismatches) > settings["max_mismatches"] or

129 len(terminal) > settings["max_three_prime_mismatches"]):

130 return None

131 possible = sorted(set(mismatches) | set(ambiguous))

132 return {

133 "mismatches": len(mismatches),

134 "mismatch_positions_1based": mismatches,

135 "mismatch_positions_from_five_prime_1based": mismatches,

136 "mismatch_positions_from_three_prime_1based": sorted(len(primer) - i + 1 for i in mismatches),

137 "three_prime_mismatches": len(terminal),

138 "three_prime_mismatch_positions_1based": terminal,

139 "three_prime_mismatch_positions_from_three_prime_1based": sorted(len(primer) - i + 1 for i in terminal),

140 "mismatches_upper_bound": len(possible),

141 "three_prime_mismatches_upper_bound": sum(i > boundary for i in possible),

142 "ambiguous_positions_1based": ambiguous,

143 "ambiguous_positions_from_five_prime_1based": ambiguous,

144 "ambiguous_positions_from_three_prime_1based": sorted(len(primer) - i + 1 for i in ambiguous),

145 "unresolved_ambiguity": bool(ambiguous),

146 }

147

148

149def make_hit(primer, record_id, sequence, start, strand, circular, settings):

150 window = reference_window(sequence, start, len(primer), circular)

151 if window is None:

152 return None

153 score = score_window(primer, window, strand, settings)

154 if score is None:

155 return None

156 start = start % len(sequence) if circular else start

157 return {

158 "record_id": record_id, "start": start, "end": start + len(primer),

159 "strand": strand, "wraps_origin": start + len(primer) > len(sequence),

160 **score,

161 }

162

163

164def exhaustive_hits(primers, references, circular, settings):

165 """Search every full-length window on both strands within declared bounds."""

166 comparisons = sum(

167 (len(seq) if record in circular and len(primer) <= len(seq)

168 else max(0, len(seq) - len(primer) + 1)) * len(primer) * 2

169 for primer in primers for record, seq in references.items()

170 )

171 if comparisons > settings["max_comparisons"]:

172 raise SearchIncomplete(

173 f"Exhaustive search requires at most {comparisons:,} base comparisons; "

174 f"cap is {settings['max_comparisons']:,}. Increase deliberately or use "

175 "the explicitly heuristic BLAST engine. No truncated screen was accepted."

176 )

177 result = {primer: [] for primer in primers}

178 total = 0

179 for primer in primers:

180 for record, sequence in references.items():

181 is_circular = record in circular

182 windows = len(sequence) if is_circular else len(sequence) - len(primer) + 1

183 if len(primer) > len(sequence):

184 continue

185 for strand in ("+", "-"):

186 for start in range(windows):

187 hit = make_hit(primer, record, sequence, start, strand, is_circular, settings)

188 if hit is not None:

189 total += 1

190 if total > settings["max_hits"]:

191 raise SearchIncomplete("Binding-hit cap exceeded; screen discarded, not truncated")

192 result[primer].append(hit)

193 return result, {"base_comparisons_upper_bound": comparisons, "tools": {}}, []

194

195

196def run_tool(command, timeout):

197 try:

198 completed = subprocess.run(

199 command, check=True, capture_output=True, text=True, timeout=timeout,

200 )

201 except subprocess.TimeoutExpired as exc:

202 raise SearchIncomplete(f"{Path(command[0]).name} timed out after {timeout}s") from exc

203 except subprocess.CalledProcessError as exc:

204 raise SearchIncomplete(

205 f"{Path(command[0]).name} failed (exit {exc.returncode}): {exc.stderr[-2000:]}"

206 ) from exc

207 except OSError as exc:

208 raise SearchIncomplete(f"Cannot execute {Path(command[0]).name}: {exc}") from exc

209 return completed.stdout.strip()

210

211

212def blast_hits(primers, references, circular, settings):

213 """Use BLAST for candidate discovery, then recheck entire ungapped primers.

214

215 Discovery is heuristic even with unsaturated reporting limits. Both ends of

216 clipped HSPs are restored before mismatch counting on the original FASTA.

217 Synthetic database IDs prevent pipe/accession conventions from changing IDs.

218 """

219 if any(len(primer) < 7 for primer in primers):

220 raise SearchIncomplete("BLAST's configured 7-base seed requires annealing cores of at least 7 bases")

221 program_paths = {}

222 for name in ("blastn", "makeblastdb"):

223 executable = shutil.which(settings[name])

224 if not executable:

225 raise SearchIncomplete(f"{name} executable not found: {settings[name]}")

226 program_paths[name] = executable

227 timeout = settings["blast_timeout"]

228 tools = {name: {"executable": path, "version": run_tool([path, "-version"], timeout)}

229 for name, path in program_paths.items()}

230 primer_map = {f"q{i}": primer for i, primer in enumerate(primers)}

231 record_map = {f"r{i}": record for i, record in enumerate(references)}

232 hits = {primer: [] for primer in primers}

233 seen, target_counts, hsp_counts = set(), defaultdict(set), Counter()

234 total_hits, rows, comparisons = 0, 0, 0

235 max_length = max(map(len, primers))

236 with tempfile.TemporaryDirectory(prefix="primer-specificity-") as temporary:

237 temporary = Path(temporary)

238 database_fasta, query_fasta = temporary / "reference.fa", temporary / "queries.fa"

239 database_path, output = temporary / "reference_db", temporary / "hits.tsv"

240 with database_fasta.open("w", encoding="utf-8") as handle:

241 for synthetic, record in record_map.items():

242 sequence = references[record]

243 if record in circular:

244 sequence += sequence[:max_length - 1]

245 handle.write(f">{synthetic}\n{sequence}\n")

246 with query_fasta.open("w", encoding="utf-8") as handle:

247 for query, primer in primer_map.items():

248 handle.write(f">{query}\n{primer}\n")

249 make_command = [program_paths["makeblastdb"], "-in", str(database_fasta), "-dbtype", "nucl",

250 "-out", str(database_path), "-parse_seqids"]

251 blast_command = [

252 program_paths["blastn"], "-task", "blastn-short", "-query", str(query_fasta),

253 "-db", str(database_path), "-out", str(output), "-strand", "both", "-ungapped",

254 "-dust", "no", "-soft_masking", "false", "-word_size", "7", "-evalue", "1000",

255 "-max_target_seqs", str(settings["blast_max_target_seqs"]),

256 "-max_hsps", str(settings["blast_max_hsps"]),

257 "-outfmt", "6 qseqid sseqid qstart qend sstart send gaps",

258 ]

259 run_tool(make_command, timeout)

260 run_tool(blast_command, timeout)

261 with output.open(encoding="utf-8") as handle:

262 for line in handle:

263 rows += 1

264 if rows > settings["max_blast_rows"]:

265 raise SearchIncomplete("BLAST output-row cap exceeded; screen discarded")

266 fields = line.rstrip("\n").split("\t")

267 if len(fields) != 7:

268 raise SearchIncomplete("Unexpected BLAST output format")

269 query, subject = fields[:2]

270 # BLAST may decorate local identifiers when -parse_seqids is used.

271 subject = subject.removeprefix("lcl|")

272 if query not in primer_map or subject not in record_map:

273 raise SearchIncomplete("Unrecognized BLAST query/reference identifier")

274 qstart, qend, sstart, send, gaps = map(int, fields[2:])

275 primer, record = primer_map[query], record_map[subject]

276 if gaps or not 1 <= qstart <= qend <= len(primer):

277 raise SearchIncomplete("Unsupported gapped/reversed-query BLAST alignment")

278 target_counts[query].add(subject)

279 hsp_counts[query, subject] += 1

280 strand = "+" if send >= sstart else "-"

281 start = sstart - qstart if strand == "+" else sstart + qstart - 1 - len(primer)

282 if record in circular:

283 start %= len(references[record])

284 key = (primer, record, start, strand)

285 if key in seen:

286 continue

287 seen.add(key)

288 comparisons += len(primer)

289 if comparisons > settings["max_comparisons"]:

290 raise SearchIncomplete("BLAST candidate realignment comparison cap exceeded")

291 hit = make_hit(primer, record, references[record], start, strand,

292 record in circular, settings)

293 if hit is not None:

294 total_hits += 1

295 if total_hits > settings["max_hits"]:

296 raise SearchIncomplete("Binding-hit cap exceeded; screen discarded")

297 hits[primer].append(hit)

298 issues = []

299 if any(len(records) >= settings["blast_max_target_seqs"] and

300 len(records) < len(references) for records in target_counts.values()):

301 issues.append("BLAST max_target_seqs reporting limit reached; targets may be omitted")

302 if any(count >= settings["blast_max_hsps"] for count in hsp_counts.values()):

303 issues.append("BLAST max_hsps reporting limit reached; binding sites may be omitted")

304 provenance = {

305 "tools": tools, "commands": [make_command, blast_command],

306 "blast_rows": rows, "candidate_intervals": len(seen),

307 "realignment_base_comparisons": comparisons,

308 }

309 return hits, provenance, issues

310

311

312def enumerate_products(pair, hits_by_sequence, references, circular, settings, expected):

313 cross_pair = pair.get("cross_pair", False)

314 sources = source_oligos(pair)

315 hits = []

316 for role, key in (("F", "forward"), ("R", "reverse")):

317 for hit in sorted(hits_by_sequence[pair[key]],

318 key=lambda h: (h["record_id"], h["start"], h["strand"])):

319 hits.append({

320 **hit, "primer_role": role, "hit_id": f"{pair['pair_id']}:h{len(hits) + 1}",

321 "source_oligo_id": sources[role]["oligo_id"],

322 "source_pair_id": sources[role]["pair_id"],

323 "source_primer_role": sources[role]["primer_role"],

324 })

325 by_record = defaultdict(lambda: {"+": [], "-": []})

326 for hit in hits:

327 by_record[hit["record_id"]][hit["strand"]].append(hit)

328 combinations = sum(len(group["+"]) * len(group["-"]) for group in by_record.values())

329 if combinations > settings["max_product_combinations"]:

330 raise SearchIncomplete(f"Pair {pair['pair_id']}: product-combination cap exceeded")

331 products = []

332 intended_coordinates = {(item["record_id"], item["start"], item["end"]) for item in expected}

333 found_expected = set()

334 for record, strands in by_record.items():

335 sequence = references[record]

336 for plus in strands["+"]:

337 for minus in strands["-"]:

338 if cross_pair and plus["primer_role"] == minus["primer_role"]:

339 continue # Single-source-oligo products are already enumerated under its original pair.

340 minus_start, end = minus["start"], minus["end"]

341 if record in circular and minus_start < plus["end"]:

342 minus_start += len(sequence)

343 end += len(sequence)

344 if minus_start < plus["end"]:

345 continue # Overlapping primer sites are not accepted as PCR products.

346 length = end - plus["start"]

347 if record in circular and length > len(sequence):

348 continue

349 if not settings["min_product"] <= length <= settings["max_product"]:

350 continue

351 kind = plus["primer_role"] + "/" + minus["primer_role"]

352 coordinates = (record, plus["start"], end)

353 intended = not cross_pair and coordinates in intended_coordinates and kind in ("F/R", "R/F")

354 if intended:

355 found_expected.add(coordinates)

356 classification = "intended" if intended else ("potential_off_target" if expected or cross_pair else "unclassified")

357 tails = {"F": pair.get("forward_tail", ""), "R": pair.get("reverse_tail", "")}

358 product_sequence = reference_window(sequence, plus["start"], length, record in circular)

359 products.append({

360 "product_id": f"{pair['pair_id']}:p{len(products) + 1}",

361 "record_id": record, "start": plus["start"], "end": end,

362 "length": length, "wraps_origin": end > len(sequence),

363 "primer_roles": kind, "left_hit_id": plus["hit_id"], "right_hit_id": minus["hit_id"],

364 "left_source_oligo_id": plus["source_oligo_id"], "right_source_oligo_id": minus["source_oligo_id"],

365 "source_primer_roles": plus["source_primer_role"] + "/" + minus["source_primer_role"],

366 "cross_pair": cross_pair,

367 "classification": classification,

368 "length_with_tails": length + len(tails[plus["primer_role"]]) + len(tails[minus["primer_role"]]),

369 "reference_product_sha256": hashlib.sha256(product_sequence.encode("ascii")).hexdigest(),

370 "unresolved_ambiguity": plus["unresolved_ambiguity"] or minus["unresolved_ambiguity"],

371 })

372 if len(products) > settings["max_products"]:

373 raise SearchIncomplete(f"Pair {pair['pair_id']}: product cap exceeded; screen discarded")

374 missing = [item for item in expected if (item["record_id"], item["start"], item["end"]) not in found_expected]

375 if cross_pair:

376 status = "potential_off_target" if products else "no_off_target_found_within_search_scope"

377 elif not expected:

378 status = "no_expected_target"

379 elif missing:

380 status = "intended_target_not_found"

381 elif any(item["classification"] == "potential_off_target" for item in products):

382 status = "potential_off_target"

383 else:

384 status = "no_off_target_found_within_search_scope"

385 unresolved = any(h["unresolved_ambiguity"] for h in hits)

386 return {

387 "pair_id": pair["pair_id"], "status": status,

388 "cross_pair": cross_pair, "source_oligos": sources,

389 "hit_count": len(hits), "product_count": len(products), "hits": hits,

390 "products": products, "expected": expected, "missing_expected": missing,

391 "potential_off_target_count": sum(p["classification"] == "potential_off_target" for p in products),

392 "unresolved_ambiguity": unresolved,

393 "product_combinations_checked": combinations,

394 **({"status_before_incomplete": status, "status": "incomplete"} if unresolved else {}),

395 }