design_primers.py · Scientific Agent Skills

1#!/usr/bin/env python3

2"""Design paired PCR primers with Primer3 and export auditable coordinates."""

3from __future__ import annotations

4

5import argparse

6import csv

7import json

8import sys

9from pathlib import Path

10

11sys.dont_write_bytecode = True

12from _common import (PAIR_FIELDS, check_output_paths, dna, finite_number, integer,

13 load_primer3, read_fasta, revcomp, sha256_file, write_json)

14

15# A deliberately explicit supported interface prevents misspelled tags being

16# silently ignored by Primer3. More advanced Primer3 tasks use its native API.

17GLOBAL_KEYS = set("""

18PRIMER_NUM_RETURN PRIMER_MIN_SIZE PRIMER_OPT_SIZE PRIMER_MAX_SIZE

19PRIMER_MIN_TM PRIMER_OPT_TM PRIMER_MAX_TM PRIMER_PAIR_MAX_DIFF_TM

20PRIMER_MIN_GC PRIMER_OPT_GC_PERCENT PRIMER_MAX_GC PRIMER_MAX_POLY_X

21PRIMER_GC_CLAMP PRIMER_MAX_END_GC PRIMER_MAX_END_STABILITY

22PRIMER_MAX_SELF_ANY_TH PRIMER_MAX_SELF_END_TH PRIMER_MAX_HAIRPIN_TH

23PRIMER_PAIR_MAX_COMPL_ANY_TH PRIMER_PAIR_MAX_COMPL_END_TH

24PRIMER_SALT_MONOVALENT PRIMER_SALT_DIVALENT PRIMER_DNTP_CONC

25PRIMER_DNA_CONC PRIMER_TM_FORMULA PRIMER_SALT_CORRECTIONS

26PRIMER_PRODUCT_SIZE_RANGE PRIMER_MIN_5_PRIME_OVERLAP_OF_JUNCTION

27PRIMER_MIN_3_PRIME_OVERLAP_OF_JUNCTION PRIMER_MAX_NS_ACCEPTED

28PRIMER_WT_TM_GT PRIMER_WT_TM_LT PRIMER_WT_SIZE_GT PRIMER_WT_SIZE_LT

29PRIMER_PAIR_WT_DIFF_TM PRIMER_PAIR_WT_PR_PENALTY

30PRIMER_MIN_QUALITY PRIMER_MIN_END_QUALITY PRIMER_QUALITY_RANGE_MIN

31PRIMER_QUALITY_RANGE_MAX PRIMER_WT_SEQ_QUAL PRIMER_WT_END_QUAL

32""".split())

33SEQUENCE_KEYS = set("""

34SEQUENCE_TARGET SEQUENCE_INCLUDED_REGION SEQUENCE_EXCLUDED_REGION

35SEQUENCE_PRIMER_PAIR_OK_REGION_LIST SEQUENCE_OVERLAP_JUNCTION_LIST

36SEQUENCE_PRIMER SEQUENCE_PRIMER_REVCOMP SEQUENCE_QUALITY

37""".split())

38

39

40def region(value, length, name):

41 if not isinstance(value, list) or len(value) != 2:

42 raise ValueError(f"{name} requires [start, length]")

43 start, size = value

44 integer(start, name + ' start')

45 integer(size, name + ' length', 1)

46 if start + size > length:

47 raise ValueError(f"{name} lies outside the selected FASTA record")

48

49

50def validate_config(config, length):

51 if not isinstance(config, dict) or set(config) - {"sequence_args", "global_args"}:

52 raise ValueError("Config accepts only sequence_args and global_args mappings")

53 sequence, global_args = config.get("sequence_args", {}), config.get("global_args", {})

54 if not isinstance(sequence, dict) or not isinstance(global_args, dict):

55 raise ValueError("sequence_args and global_args must be JSON objects")

56 for values, allowed, label in [(sequence, SEQUENCE_KEYS, "sequence_args"), (global_args, GLOBAL_KEYS, "global_args")]:

57 if set(values) - allowed:

58 raise ValueError(f"Unsupported {label}: {sorted(set(values) - allowed)}; see input-contract.md")

59 for key in ("SEQUENCE_TARGET", "SEQUENCE_EXCLUDED_REGION"):

60 if key in sequence:

61 if not isinstance(sequence[key], list):

62 raise ValueError(f"{key} must be a list of [start, length] intervals")

63 for value in sequence[key]:

64 region(value, length, key)

65 if "SEQUENCE_INCLUDED_REGION" in sequence:

66 region(sequence["SEQUENCE_INCLUDED_REGION"], length, "SEQUENCE_INCLUDED_REGION")

67 for key in ("SEQUENCE_PRIMER", "SEQUENCE_PRIMER_REVCOMP"):

68 if key in sequence:

69 sequence[key] = dna(sequence[key], key)

70 if not 15 <= len(sequence[key]) <= 36:

71 raise ValueError(f"{key}: supported annealing core length is 15..36")

72 if "SEQUENCE_OVERLAP_JUNCTION_LIST" in sequence:

73 junctions = sequence["SEQUENCE_OVERLAP_JUNCTION_LIST"]

74 if not isinstance(junctions, list):

75 raise ValueError("SEQUENCE_OVERLAP_JUNCTION_LIST must be a list")

76 for value in junctions:

77 integer(value, "junction")

78 if value >= length - 1:

79 raise ValueError("Junction must name the base immediately left of an internal template boundary")

80 if "SEQUENCE_QUALITY" in sequence:

81 quality = sequence["SEQUENCE_QUALITY"]

82 if not isinstance(quality, list) or len(quality) != length:

83 raise ValueError("SEQUENCE_QUALITY needs one integer per template base")

84 for value in quality:

85 integer(value, "quality")

86 if "SEQUENCE_PRIMER_PAIR_OK_REGION_LIST" in sequence:

87 rows = sequence["SEQUENCE_PRIMER_PAIR_OK_REGION_LIST"]

88 if not isinstance(rows, list):

89 raise ValueError("SEQUENCE_PRIMER_PAIR_OK_REGION_LIST must be a list")

90 for row in rows:

91 if not isinstance(row, list) or len(row) != 4:

92 raise ValueError("OK region requires [left_start,left_length,right_start,right_length]")

93 for value in (row[:2], row[2:]):

94 if value != [-1, -1]:

95 region(value, length, "OK region")

96 for key, value in global_args.items():

97 if key == "PRIMER_PRODUCT_SIZE_RANGE":

98 if not isinstance(value, list) or not value:

99 raise ValueError("Product ranges must be a nonempty list of [minimum,maximum]")

100 for limits in value:

101 if not isinstance(limits, list) or len(limits) != 2:

102 raise ValueError("Product range requires [minimum,maximum]")

103 integer(limits[0], "minimum product", 1)

104 integer(limits[1], "maximum product", limits[0])

105 else:

106 finite_number(value, key)

107 if global_args.get("PRIMER_MAX_NS_ACCEPTED", 0) != 0:

108 raise ValueError("This workflow requires unambiguous primers; PRIMER_MAX_NS_ACCEPTED must be 0")

109 return sequence, global_args

110

111

112def bed_exclusions(path, records, record_id):

113 exclusions = []

114 with Path(path).open(encoding="utf-8") as handle:

115 for number, raw in enumerate(handle, 1):

116 if not raw.strip() or raw.startswith(("#", "track ", "browser ")):

117 continue

118 parts = raw.rstrip().split("\t")

119 if len(parts) < 3 or parts[0] not in records:

120 raise ValueError(f"BED line {number}: unknown record or fewer than 3 tab-separated columns")

121 start, end = int(parts[1]), int(parts[2])

122 region([start, end - start], len(records[parts[0]]), "BED interval")

123 if parts[0] == record_id:

124 exclusions.append([start, end - start])

125 return exclusions

126

127

128def design(args):

129 records = read_fasta(args.template)

130 record_id = args.record or (next(iter(records)) if len(records) == 1 else None)

131 if record_id not in records:

132 raise ValueError("Choose --record from the FASTA IDs when the template has multiple records")

133 original = records[record_id]

134 # Primer3's sequence alphabet is narrower than IUPAC. Unknown positions are

135 # made explicitly unavailable, never guessed into a concrete nucleotide.

136 template = "".join(base if base in "ACGT" else "N" for base in original)

137 config = json.loads(Path(args.config).read_text()) if args.config else {}

138 sequence_args, overrides = validate_config(config, len(template))

139 global_args = {

140 "PRIMER_TASK": "generic", "PRIMER_PICK_LEFT_PRIMER": 1,

141 "PRIMER_PICK_RIGHT_PRIMER": 1, "PRIMER_PICK_INTERNAL_OLIGO": 0,

142 "PRIMER_FIRST_BASE_INDEX": 0, "PRIMER_EXPLAIN_FLAG": 1,

143 "PRIMER_THERMODYNAMIC_OLIGO_ALIGNMENT": 1,

144 "PRIMER_NUM_RETURN": 5, "PRIMER_MIN_SIZE": 18, "PRIMER_OPT_SIZE": 20,

145 "PRIMER_MAX_SIZE": 25, "PRIMER_MIN_TM": 57.0, "PRIMER_OPT_TM": 60.0,

146 "PRIMER_MAX_TM": 63.0, "PRIMER_PAIR_MAX_DIFF_TM": 3.0,

147 "PRIMER_MIN_GC": 20.0, "PRIMER_MAX_GC": 80.0,

148 "PRIMER_MAX_POLY_X": 4, "PRIMER_MAX_NS_ACCEPTED": 0,

149 "PRIMER_SALT_MONOVALENT": 50.0, "PRIMER_SALT_DIVALENT": 1.5,

150 "PRIMER_DNTP_CONC": 0.6, "PRIMER_DNA_CONC": 50.0,

151 "PRIMER_TM_FORMULA": 1, "PRIMER_SALT_CORRECTIONS": 1,

152 "PRIMER_PRODUCT_SIZE_RANGE": [[70, 200]] if args.preset == "qpcr" else [[100, 1000]],

153 }

154 global_args.update(overrides)

155 for key in ("PRIMER_NUM_RETURN", "PRIMER_MIN_SIZE", "PRIMER_OPT_SIZE", "PRIMER_MAX_SIZE", "PRIMER_MAX_POLY_X", "PRIMER_GC_CLAMP", "PRIMER_MAX_END_GC", "PRIMER_TM_FORMULA", "PRIMER_SALT_CORRECTIONS", "PRIMER_MIN_5_PRIME_OVERLAP_OF_JUNCTION", "PRIMER_MIN_3_PRIME_OVERLAP_OF_JUNCTION", "PRIMER_MAX_NS_ACCEPTED", "PRIMER_MIN_QUALITY", "PRIMER_MIN_END_QUALITY", "PRIMER_QUALITY_RANGE_MIN", "PRIMER_QUALITY_RANGE_MAX"):

156 if key in global_args:

157 integer(global_args[key], key)

158 if global_args["PRIMER_TM_FORMULA"] not in (0, 1):

159 raise ValueError("PRIMER_TM_FORMULA must be 0 (Breslauer) or 1 (SantaLucia)")

160 if global_args["PRIMER_SALT_CORRECTIONS"] not in (0, 1, 2):

161 raise ValueError("PRIMER_SALT_CORRECTIONS must be 0 (Schildkraut), 1 (SantaLucia), or 2 (Owczarzy)")

162 if not 1 <= global_args["PRIMER_NUM_RETURN"] <= 100:

163 raise ValueError("PRIMER_NUM_RETURN must be 1..100")

164 if not 15 <= global_args["PRIMER_MIN_SIZE"] <= global_args["PRIMER_OPT_SIZE"] <= global_args["PRIMER_MAX_SIZE"] <= 36:

165 raise ValueError("Require 15 <= MIN_SIZE <= OPT_SIZE <= MAX_SIZE <= 36")

166 if not global_args["PRIMER_MIN_TM"] <= global_args["PRIMER_OPT_TM"] <= global_args["PRIMER_MAX_TM"]:

167 raise ValueError("Require MIN_TM <= OPT_TM <= MAX_TM")

168 if not 0 <= global_args["PRIMER_MIN_GC"] <= global_args["PRIMER_MAX_GC"] <= 100:

169 raise ValueError("Require 0 <= MIN_GC <= MAX_GC <= 100")

170 finite_number(global_args["PRIMER_DNA_CONC"], "PRIMER_DNA_CONC", positive=True)

171 if args.mask_bed:

172 sequence_args.setdefault("SEQUENCE_EXCLUDED_REGION", []).extend(bed_exclusions(args.mask_bed, records, record_id))

173 sequence_args.update(SEQUENCE_ID=record_id, SEQUENCE_TEMPLATE=template)

174 tails = {name: dna(getattr(args, name), name, empty=True) for name in ("forward_tail", "reverse_tail")}

175 p3 = load_primer3()

176 result = p3.bindings.design_primers(sequence_args, global_args)

177 if result.get("PRIMER_ERROR"):

178 raise ValueError(result["PRIMER_ERROR"])

179 pairs = []

180 for i in range(result.get("PRIMER_PAIR_NUM_RETURNED", 0)):

181 left, left_length = result[f"PRIMER_LEFT_{i}"]

182 right, right_length = result[f"PRIMER_RIGHT_{i}"]

183 forward, reverse = result[f"PRIMER_LEFT_{i}_SEQUENCE"], result[f"PRIMER_RIGHT_{i}_SEQUENCE"]

184 if template[left:left + left_length] != forward or revcomp(template[right - right_length + 1:right + 1]) != reverse:

185 raise ValueError("Primer3 output failed template/orientation verification")

186 product = template[left:right + 1]

187 if len(product) != result[f"PRIMER_PAIR_{i}_PRODUCT_SIZE"]:

188 raise ValueError("Primer3 output failed product-length verification")

189 pairs.append({

190 "pair_id": f"pair_{i + 1}", "forward": forward, "reverse": reverse, **tails,

191 "forward_order_sequence": tails["forward_tail"] + forward,

192 "reverse_order_sequence": tails["reverse_tail"] + reverse,

193 "record_id": record_id, "forward_interval": [left, left + left_length],

194 "reverse_interval": [right - right_length + 1, right + 1],

195 "product_interval": [left, right + 1], "product_size": len(product),

196 "tailed_product_size": len(product) + sum(map(len, tails.values())),

197 "product_sequence": product,

198 "tailed_product_sequence": tails["forward_tail"] + product + revcomp(tails["reverse_tail"]),

199 "forward_tm_c": result[f"PRIMER_LEFT_{i}_TM"],

200 "reverse_tm_c": result[f"PRIMER_RIGHT_{i}_TM"],

201 "penalty": result[f"PRIMER_PAIR_{i}_PENALTY"],

202 "specificity_status": "not_screened",

203 })

204 report = {

205 "schema_version": "1.0", "status": "candidates_generated" if pairs else "no_candidates",

206 "coordinate_system": "0-based half-open on supplied FASTA record; reverse oligos are 5-prime to 3-prime",

207 "provenance": {"template_file": str(args.template), "template_sha256": sha256_file(args.template),

208 "record_id": record_id, "template_length": len(template),

209 "primer3_py_version": p3.__version__,

210 "libprimer3_version": p3.thermoanalysis.get_libprimer3_version(),

211 "python_version": sys.version.split()[0],

212 "config_sha256": sha256_file(args.config) if args.config else None,

213 "mask_bed_sha256": sha256_file(args.mask_bed) if args.mask_bed else None},

214 "settings": {"preset": args.preset, "sequence_args": sequence_args, "global_args": global_args, "tails": tails},

215 "ambiguous_template_bases_masked": sum(c not in "ACGT" for c in original),

216 "pairs": pairs,

217 "engine_explanations": {key: value for key, value in result.items() if key.endswith("EXPLAIN") or key in ("PRIMER_WARNING", "PRIMER_ERROR")},

218 "limitations": ["Primer3 penalty is a design ranking, not experimental validation or genomic specificity.",

219 "Coordinates refer to the supplied template, not a genome unless the FASTA itself is that reference.",

220 "Defaults are starting conditions; record the actual polymerase/buffer and validate experimentally.",

221 "Tails were appended after design; assess full oligos separately for secondary structures."],

222 }

223 write_json(args.output, report)

224 for path in (args.pairs_out, args.expected_out):

225 if path:

226 Path(path).parent.mkdir(parents=True, exist_ok=True)

227 if args.pairs_out:

228 with Path(args.pairs_out).open("w", newline="", encoding="utf-8") as handle:

229 writer = csv.DictWriter(handle, fieldnames=PAIR_FIELDS, delimiter="\t", extrasaction="ignore")

230 writer.writeheader()

231 writer.writerows(pairs)

232 if args.expected_out:

233 with Path(args.expected_out).open("w", newline="", encoding="utf-8") as handle:

234 writer = csv.writer(handle, delimiter="\t")

235 writer.writerow(["pair_id", "record_id", "start", "end"])

236 writer.writerows([p["pair_id"], record_id, *p["product_interval"]] for p in pairs)

237 return 0 if pairs else 1

238

239

240def main(argv=None):

241 parser = argparse.ArgumentParser(description=__doc__)

242 parser.add_argument("--template", required=True, type=Path, help="Local FASTA; no remote retrieval")

243 parser.add_argument("--record", help="Exact FASTA ID; required for multi-record input")

244 parser.add_argument("--preset", choices=("pcr", "qpcr"), default="pcr")

245 parser.add_argument("--config", type=Path, help="JSON sequence_args/global_args overrides")

246 parser.add_argument("--mask-bed", type=Path, help="Exclude primer binding in BED intervals in template coordinates")

247 parser.add_argument("--forward-tail", default="", help="5-prime tail appended to forward core")

248 parser.add_argument("--reverse-tail", default="", help="5-prime tail appended to reverse core")

249 parser.add_argument("--output", required=True, type=Path, help="JSON design report")

250 parser.add_argument("--pairs-out", type=Path, help="TSV for downstream thermodynamics/screening")

251 parser.add_argument("--expected-out", type=Path, help="TSV intended product coordinates on this template")

252 args = parser.parse_args(argv)

253 try:

254 check_output_paths([args.template, args.config, args.mask_bed], [args.output, args.pairs_out, args.expected_out])

255 return design(args)

256 except (ValueError, OSError, RuntimeError, TypeError, OverflowError) as exc:

257 print(f"error: {exc}", file=sys.stderr)

258 return 2

259

260

261if __name__ == "__main__":

262 raise SystemExit(main())