File size: 10,698 Bytes
81ae663 | 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 | """Reaction execution, canonical outcomes, parent metadata, and route replay.
Templates declare graph transformations. Experimental yields and substrate
scope require separate chemical validation. The shipped examples use named
public compounds and are not an Enamine inventory.
"""
from pathlib import Path
import csv, json, re, hashlib
import numpy as np
from rdkit import Chem, DataStructs, RDLogger
from rdkit.Chem import AllChem, Descriptors, QED, rdFingerprintGenerator
from .graph import Graph, Node, Edge
RDLogger.DisableLog("rdApp.warning")
FP = rdFingerprintGenerator.GetMorganGenerator(radius=2, fpSize=128)
def canonical(smiles):
"""Return canonical isomeric SMILES for one valid connected molecule."""
mol = Chem.MolFromSmiles(str(smiles))
if mol is None:
raise ValueError(f"Invalid SMILES: {smiles}")
if len(Chem.GetMolFrags(mol)) != 1:
raise ValueError("A connected molecule is required")
Chem.SanitizeMol(mol)
return Chem.MolToSmiles(mol, isomericSmiles=True)
def features(smiles, depth, budget, terminal=False):
"""Return 131 features comprising 128 Morgan bits and state indicators."""
m = Chem.MolFromSmiles(smiles)
fp = np.zeros(128, dtype=np.float32)
DataStructs.ConvertToNumpyArray(FP.GetFingerprint(m), fp)
return fp.tolist() + [depth / max(budget, 1), float(terminal), 0.0]
def read_catalog(path, limit=None):
"""CSV/TSV input with smiles and id, or native SMILES and Enamine_ID."""
if limit is not None and (int(limit) != limit or limit < 1):
raise ValueError("Catalog row limit must be a positive integer")
p = Path(path)
import gzip
opener = gzip.open if p.suffix == ".gz" else open
rows = []
seen = {}
with opener(p, "rt", newline="") as f:
sample = f.read(8192)
f.seek(0)
if not sample.strip():
raise ValueError("Empty catalog file")
delimiter = "\t" if "\t" in sample.splitlines()[0] else ","
for i, r in enumerate(csv.DictReader(f, delimiter=delimiter)):
r = {k.lower(): v for k, v in r.items()}
smi = canonical(r.get("smiles", ""))
cid = r.get(
"id", r.get("enamine_id", r.get("catalog_id", r.get("code", "")))
)
if not cid:
raise ValueError(f"Missing catalog identifier at row {i+2}")
if cid in seen and seen[cid] != smi:
raise ValueError(f"Conflicting structures for {cid}")
seen[cid] = smi
meta = {k: v for k, v in r.items() if k != "smiles"}
match = re.fullmatch(r"([sm])_(\d+)_([\d_]+)", cid)
if match:
meta.update(
chemistry_class=match[1],
parent_reaction_id=match[2],
parent_reagent_ids=match[3].split("_"),
)
rows.append({"id": cid, "smiles": smi, "metadata": meta})
if limit is not None and len(rows) >= limit:
break
if not rows:
raise ValueError("Empty parent catalog")
# Canonical molecular identity avoids duplicated parent action descriptions.
grouped = {}
for r in rows:
if r["smiles"] in grouped:
grouped[r["smiles"]]["metadata"].setdefault("aliases", []).append(r["id"])
else:
grouped[r["smiles"]] = r
return list(grouped.values())
def reactions(path):
"""Read two-reactant SMARTS templates with unique IDs and finite costs."""
data = json.loads(Path(path).read_text())
items = []
if len({r["id"] for r in data}) != len(data):
raise ValueError("Reaction template IDs must be unique")
for r in data:
rxn = AllChem.ReactionFromSmarts(r["smarts"])
if rxn is None:
raise ValueError("Invalid reaction SMARTS")
if rxn.GetNumReactantTemplates() != 2:
raise ValueError("Templates require parent and one reagent")
cost = float(r.get("cost", 1))
if cost < 0 or not np.isfinite(cost):
raise ValueError("Invalid reaction cost")
items.append((r, rxn))
return items
def products(rxn, parent, reagent):
"""Enumerate unique sanitized products with molecular weight at most 750 Da."""
result = set()
for ps in rxn.RunReactants(
(Chem.MolFromSmiles(parent), Chem.MolFromSmiles(reagent))
):
if len(ps) != 1:
continue
try:
Chem.SanitizeMol(ps[0])
s = canonical(Chem.MolToSmiles(ps[0]))
if Descriptors.MolWt(ps[0]) <= 750:
result.add(s)
except (ValueError, RuntimeError):
continue
return sorted(result)
def build_graph(parents, reagents, templates, budget=2, max_nodes=100000):
"""Expand the supplied inventory through at most budget additional reactions."""
if int(budget) != budget or budget < 0:
raise ValueError("Reaction budget must be a nonnegative integer")
nodes = [Node("root", features=[0.0] * 130 + [1.0])]
edges = []
states = {}
outcomes = {}
def state(smi, k):
key = (smi, k)
if key not in states:
nid = "s" + str(len(states))
states[key] = nid
nodes.append(
Node(
nid,
features=features(smi, k, budget),
metadata={"smiles": smi, "steps": k},
)
)
if len(nodes) > max_nodes:
raise ValueError(
"Graph limit exceeded. Select fewer parents or increase max_nodes explicitly."
)
return states[key]
for p in parents:
sid = state(p["smiles"], 0)
edges.append(
Edge(
"e" + str(len(edges)),
"root",
sid,
0,
{
"kind": "parent",
"parent_id": p["id"],
"smiles": p["smiles"],
"metadata": p.get("metadata", {}),
},
)
)
for k in range(budget + 1):
current = [(s, nid) for (s, depth), nid in states.items() if depth == k]
for smi, sid in current:
if smi not in outcomes:
tid = "t" + str(len(outcomes))
outcomes[smi] = tid
nodes.append(Node(tid, smi, features(smi, 0, budget, True)))
if len(nodes) > max_nodes:
raise ValueError("Graph limit exceeded during terminal creation")
edges.append(
Edge("e" + str(len(edges)), sid, outcomes[smi], 0, {"kind": "stop"})
)
if k == budget:
continue
seen_actions = set()
for spec, rxn in templates:
for reagent in reagents:
for prod in products(rxn, smi, reagent["smiles"]):
if prod == smi:
continue
signature = (spec["id"], reagent["id"], prod)
if signature in seen_actions:
continue
seen_actions.add(signature)
tid = state(prod, k + 1)
edges.append(
Edge(
"e" + str(len(edges)),
sid,
tid,
float(spec.get("cost", 1)),
{
"kind": "reaction",
"template_id": spec["id"],
"smarts": spec["smarts"],
"reagent_id": reagent["id"],
"reagent_smiles": reagent["smiles"],
"input": smi,
"product": prod,
},
)
)
return Graph(
nodes,
edges,
metadata={
"kind": "reaction",
"budget": budget,
"parent_count": len(parents),
"reagent_count": len(reagents),
"template_ids": [x[0]["id"] for x in templates],
"scope": "Complete expansion of the supplied parents, reagents, templates and molecular-weight cap",
"max_molecular_weight": 750,
},
)
def replay(route, budget, parents=None):
"""Re-execute one parent-to-product route under an additional-reaction cap.
Actions must begin with one parent selection and end with one stop. Pass
parsed parent records to check the parent identifier against the inventory.
Malformed records and unsuccessful reaction replay return False.
"""
try:
actions = route["actions"]
if budget < 0 or int(budget) != budget or len(actions) < 2:
return False
if actions[0]["kind"] != "parent" or actions[-1]["kind"] != "stop":
return False
current = canonical(actions[0]["smiles"])
if parents is not None:
inventory = {p["id"]: p["smiles"] for p in parents}
if inventory.get(actions[0]["parent_id"]) != current:
return False
if len(actions) - 2 > budget:
return False
for action in actions[1:-1]:
if action["kind"] != "reaction" or current != action["input"]:
return False
rxn = AllChem.ReactionFromSmarts(action["smarts"])
if rxn is None or rxn.GetNumReactantTemplates() != 2:
return False
if action["product"] not in products(
rxn, current, action["reagent_smiles"]
):
return False
current = action["product"]
return current == canonical(route["outcome"])
except (KeyError, TypeError, ValueError, RuntimeError, IndexError):
return False
def descriptor_rewards(graph, weights=(2.0, 1.0, 1.0)):
"""Return log rewards from QED, logP proximity, and molecular-size utility."""
w = np.asarray(weights, dtype=float)
if w.shape != (3,) or np.any(w < 0) or not np.all(np.isfinite(w)):
raise ValueError("Three nonnegative weights required")
out = {}
for s in graph.terminals:
m = Chem.MolFromSmiles(s)
u = [
QED.qed(m),
np.exp(-(((Descriptors.MolLogP(m) - 2) / 2) ** 2)),
np.exp(-(((Descriptors.MolWt(m) - 300) / 200) ** 2)),
]
out[s] = float(w @ u)
return out
|