check_thermodynamics.py · Scientific Agent Skills

1#!/usr/bin/env python3

2"""Assess annealing cores and full tailed oligos using explicit solution chemistry."""

3from __future__ import annotations

4

5import argparse

6import itertools

7import math

8import sys

9from pathlib import Path

10

11sys.dont_write_bytecode = True

12from _common import check_output_paths, finite_number, load_primer3, read_pairs, sha256_file, write_json

13

14

15def structure(result):

16 values = [result.tm, result.dg, result.dh, result.ds]

17 if any(not math.isfinite(v) for v in values):

18 raise ValueError("Thermodynamic engine returned nonfinite values")

19 return {"structure_found": bool(result.structure_found), "tm_c": result.tm,

20 "delta_g_kcal_per_mol": result.dg / 1000,

21 "delta_h_kcal_per_mol": result.dh / 1000,

22 "delta_s_cal_per_mol_k": result.ds,

23 "ascii_structure": result.ascii_structure}

24

25

26def assess(args):

27 pairs = read_pairs(args.pairs)

28 p3 = load_primer3()

29 chemistry = {key: finite_number(getattr(args, key), key, positive=(key == "dna_conc"))

30 for key in ("mv_conc", "dv_conc", "dntp_conc", "dna_conc")}

31 finite_number(args.temp_c, "temp_c", -20)

32 if args.temp_c > 100:

33 raise ValueError("temp_c must be between -20 and 100 Celsius")

34 if args.max_interactions < 1:

35 raise ValueError("max_interactions must be positive")

36 count = 2 * len(pairs)

37 interaction_count = count * (count - 1) // 2 if args.multiplex else len(pairs)

38 if interaction_count > args.max_interactions:

39 raise ValueError(f"{interaction_count} dimer comparisons exceed --max-interactions; partition deliberately or raise the limit")

40 thermo_args = {**chemistry, "temp_c": args.temp_c}

41 records, pair_summaries, warnings = [], [], []

42 if chemistry["dntp_conc"] >= chemistry["dv_conc"] and chemistry["dv_conc"] > 0:

43 warnings.append("Total dNTP is at least divalent concentration; verify the buffer/free-Mg assumptions used by the salt model.")

44 for pair in pairs:

45 entries = []

46 for role in ("forward", "reverse"):

47 core, tail = pair[role], pair[f"{role}_tail"]

48 full = tail + core

49 entry = {"oligo_id": f"{pair['pair_id']}:{role}", "pair_id": pair["pair_id"],

50 "role": role, "annealing_core": core, "tail": tail,

51 "order_sequence": full, "length": len(full), "core_length": len(core),

52 "core_gc_percent": 100 * sum(c in "GC" for c in core) / len(core)}

53 if not 2 <= len(core) <= 60:

54 raise ValueError(f"{entry['oligo_id']}: core length must be 2..60 for this nearest-neighbor workflow")

55 entry["annealing_tm_c"] = p3.calc_tm(core, **chemistry, tm_method=args.tm_method, salt_corrections_method=args.salt_corrections_method)

56 if not math.isfinite(entry["annealing_tm_c"]) or entry["annealing_tm_c"] <= -273.15:

57 raise ValueError(

58 "Thermodynamic engine returned an invalid core Tm (nonfinite or at/below absolute zero); "

59 "check salt, dNTP and DNA concentrations and their units"

60 )

61 entry["core_hairpin"] = structure(p3.calc_hairpin(core, **thermo_args, output_structure=True))

62 entry["core_homodimer"] = structure(p3.calc_homodimer(core, **thermo_args, output_structure=True))

63 entry["core_self_three_prime"] = structure(p3.calc_end_stability(core, core, **thermo_args))

64 if len(full) <= 60:

65 entry["full_oligo_status"] = "computed"

66 entry["full_hairpin"] = structure(p3.calc_hairpin(full, **thermo_args, output_structure=True))

67 entry["full_homodimer"] = structure(p3.calc_homodimer(full, **thermo_args, output_structure=True))

68 entry["full_self_three_prime"] = structure(p3.calc_end_stability(full, full, **thermo_args))

69 else:

70 entry["full_oligo_status"] = "unsupported_length"

71 warnings.append(f"{entry['oligo_id']}: full oligo exceeds 60 nt; full-length secondary structure is unresolved, not truncated.")

72 records.append(entry)

73 entries.append(entry)

74 pair_summaries.append({"pair_id": pair["pair_id"],

75 "core_tm_difference_c": abs(entries[0]["annealing_tm_c"] - entries[1]["annealing_tm_c"])})

76 comparisons = itertools.combinations(records, 2) if args.multiplex else ((records[i], records[i + 1]) for i in range(0, len(records), 2))

77 interactions = []

78 for a, b in comparisons:

79 interaction = {"oligo_a": a["oligo_id"], "oligo_b": b["oligo_id"],

80 "cross_pair": a["pair_id"] != b["pair_id"]}

81 if max(a["length"], b["length"]) > 60:

82 interaction["status"] = "unsupported_length"

83 else:

84 seq_a, seq_b = a["order_sequence"], b["order_sequence"]

85 interaction.update(status="computed",

86 heterodimer=structure(p3.calc_heterodimer(seq_a, seq_b, **thermo_args, output_structure=True)),

87 a_three_prime_to_b=structure(p3.calc_end_stability(seq_a, seq_b, **thermo_args)),

88 b_three_prime_to_a=structure(p3.calc_end_stability(seq_b, seq_a, **thermo_args)))

89 interactions.append(interaction)

90 complete = all(r["full_oligo_status"] == "computed" for r in records)

91 report = {"schema_version": "1.0", "status": "computed" if complete else "incomplete",

92 "provenance": {"pairs_file": str(args.pairs), "pairs_sha256": sha256_file(args.pairs),

93 "primer3_py_version": p3.__version__,

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

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

96 "settings": {**chemistry, "temp_c": args.temp_c, "multiplex": args.multiplex,

97 "tm_method": args.tm_method, "salt_corrections_method": args.salt_corrections_method,

98 "secondary_structure_model": "libprimer3 thermodynamic alignment defaults",

99 "max_interactions": args.max_interactions},

100 "units": {"mv_conc": "mM monovalent cations", "dv_conc": "mM total divalent cations",

101 "dntp_conc": "mM total dNTPs", "dna_conc": "nM oligonucleotide concentration under Primer3 conventions",

102 "temp_c": "Celsius for delta-G; not annealing temperature selection"},

103 "oligos": records, "pairs": pair_summaries, "interactions": interactions, "warnings": warnings,

104 "limitations": ["Predictions do not establish PCR performance or a universal pass/fail cutoff.",

105 "Core Tm excludes untemplated tails; secondary structures use complete ordered oligos.",

106 "End-stability calculations are directional and are reported in both orders.",

107 "Modified bases, degenerate mixtures, unequal multiplex concentrations, and oligos >60 nt require a suitable external model.",

108 "Multiplex dimer checks do not screen cross-pair genomic amplification."]}

109 write_json(args.output, report)

110 return 0 if complete else 1

111

112

113def main(argv=None):

114 parser = argparse.ArgumentParser(description=__doc__)

115 parser.add_argument("--pairs", required=True, type=Path)

116 parser.add_argument("--output", required=True, type=Path)

117 parser.add_argument("--mv-conc", type=float, default=50.0, help="Monovalent cations, mM")

118 parser.add_argument("--dv-conc", type=float, default=1.5, help="Total divalent cations, mM")

119 parser.add_argument("--dntp-conc", type=float, default=0.6, help="Total dNTPs, mM (sum of all four)")

120 parser.add_argument("--dna-conc", type=float, default=50.0, help="Oligonucleotide concentration, nM")

121 parser.add_argument("--temp-c", type=float, default=37.0, help="Temperature at which delta-G is calculated")

122 parser.add_argument("--tm-method", choices=("breslauer", "santalucia"), default="santalucia")

123 parser.add_argument("--salt-corrections-method", choices=("schildkraut", "santalucia", "owczarzy"), default="santalucia")

124 parser.add_argument("--multiplex", action="store_true", help="Check all unordered oligo combinations, including between pairs")

125 parser.add_argument("--max-interactions", type=int, default=10000)

126 args = parser.parse_args(argv)

127 try:

128 check_output_paths([args.pairs], [args.output])

129 return assess(args)

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

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

132 return 2

133

134

135if __name__ == "__main__":

136 raise SystemExit(main())