Spaces:
Running
Running
Upload folder using huggingface_hub (part 2)
Browse files- app/services/workflow_manifest.py +66 -0
- app/tools/castp.py +27 -64
- app/tools/function_predict.py +3 -1
- app/tools/oligo_qc.py +71 -149
- app/tools/primer_advanced.py +6 -5
- app/tools/restriction_map.py +54 -0
- app/tools/structure_prep.py +12 -56
- requirements.txt +1 -1
- tests/test_castp_provenance.py +69 -33
- tests/test_function_research_grade.py +22 -0
- tests/test_primer_advanced.py +1 -1
- tests/test_primers_qc.py +40 -6
- tests/test_restriction_map.py +73 -0
- tests/test_seq_tools.py +6 -5
- tests/test_workflow_manifest.py +52 -0
app/services/workflow_manifest.py
ADDED
|
@@ -0,0 +1,66 @@
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
| 1 |
+
"""Portable evidence export. No workflow adapter or execution is implied."""
|
| 2 |
+
from __future__ import annotations
|
| 3 |
+
|
| 4 |
+
import hashlib
|
| 5 |
+
import json
|
| 6 |
+
|
| 7 |
+
|
| 8 |
+
_STATES = {"complete": "COMPLETED", "completed": "COMPLETED", "failed": "FAILED",
|
| 9 |
+
"skipped": "SKIPPED", "running": "RUNNING", "queued": "QUEUED",
|
| 10 |
+
"pending": "NOT_EXECUTED"}
|
| 11 |
+
|
| 12 |
+
|
| 13 |
+
def build_workflow_manifest(experiment: dict, provenance: list[dict]) -> dict:
|
| 14 |
+
"""Export stored evidence without synthesizing versions or scientific results."""
|
| 15 |
+
nodes = []
|
| 16 |
+
gaps = []
|
| 17 |
+
for field in ("input_hash", "git_commit", "container_hash"):
|
| 18 |
+
if not experiment.get(field):
|
| 19 |
+
gaps.append({"node": None, "field": field, "reason": "NOT_RECORDED"})
|
| 20 |
+
for index, source in enumerate(provenance):
|
| 21 |
+
node_id = source.get("node_id") or source.get("id")
|
| 22 |
+
if not node_id:
|
| 23 |
+
node_id = f"unidentified-record-{index}"
|
| 24 |
+
gaps.append({"node": node_id, "field": "node_id", "reason": "NOT_RECORDED"})
|
| 25 |
+
raw_status = source.get("status")
|
| 26 |
+
node = {"id": node_id, "recorded_status": raw_status,
|
| 27 |
+
"state": _STATES.get(raw_status, "UNKNOWN"),
|
| 28 |
+
"tool": source.get("tool"), "tool_version": source.get("tool_version"),
|
| 29 |
+
"database": source.get("database"), "database_version": source.get("database_version"),
|
| 30 |
+
"parameters": source.get("params"), "depends_on": source.get("deps") or [],
|
| 31 |
+
"input_ref": source.get("input_ref"), "output_ref": source.get("output_ref"),
|
| 32 |
+
"completed_at": source.get("completed_at"), "evidence": source.get("evidence")}
|
| 33 |
+
nodes.append(node)
|
| 34 |
+
if node["state"] == "COMPLETED":
|
| 35 |
+
for field in ("tool", "tool_version", "completed_at"):
|
| 36 |
+
if not node[field]:
|
| 37 |
+
gaps.append({"node": node_id, "field": field, "reason": "NOT_RECORDED"})
|
| 38 |
+
if node["database"] and not node["database_version"]:
|
| 39 |
+
gaps.append({"node": node_id, "field": "database_version", "reason": "NOT_RECORDED"})
|
| 40 |
+
ids = {node["id"] for node in nodes}
|
| 41 |
+
if len(ids) != len(nodes):
|
| 42 |
+
gaps.append({"node": None, "field": "node_id", "reason": "DUPLICATE_ID"})
|
| 43 |
+
edges = []
|
| 44 |
+
for node in nodes:
|
| 45 |
+
for dependency in node["depends_on"]:
|
| 46 |
+
edges.append({"from": dependency, "to": node["id"]})
|
| 47 |
+
if dependency not in ids:
|
| 48 |
+
gaps.append({"node": node["id"], "field": "depends_on", "reason": "UNRECORDED_DEPENDENCY", "dependency": dependency})
|
| 49 |
+
if not nodes:
|
| 50 |
+
gaps.append({"node": None, "field": "steps", "reason": "NOT_RECORDED"})
|
| 51 |
+
manifest = {"format": "bionexus-workflow-evidence", "schema_version": 1,
|
| 52 |
+
"execution": {"executable": False, "adapter": None,
|
| 53 |
+
"statement": "Evidence export only; no tools are launched and external workflow compatibility has not been validated."},
|
| 54 |
+
"experiment_id": experiment.get("experiment_id"), "job_id": experiment.get("job_id"),
|
| 55 |
+
"pipeline": experiment.get("pipeline"), "recorded_status": experiment.get("status"),
|
| 56 |
+
"input_sha256": experiment.get("input_hash"), "output_sha256": experiment.get("output_hash"),
|
| 57 |
+
"git_commit": experiment.get("git_commit"), "container_hash": experiment.get("container_hash"),
|
| 58 |
+
"software_versions": experiment.get("software_versions"), "database_versions": experiment.get("database_versions"),
|
| 59 |
+
"environment": experiment.get("environment"), "random_seed": experiment.get("random_seed"),
|
| 60 |
+
"parameters": experiment.get("parameters"), "started_at": experiment.get("started_at"),
|
| 61 |
+
"finished_at": experiment.get("finished_at"), "steps": nodes, "edges": edges,
|
| 62 |
+
"reproducibility_gaps": gaps,
|
| 63 |
+
"integrity": {"algorithm": "sha256", "serialization": "JSON sorted keys, UTF-8, compact separators; manifest_sha256 excluded"}}
|
| 64 |
+
canonical = json.dumps(manifest, sort_keys=True, separators=(",", ":"), ensure_ascii=False)
|
| 65 |
+
manifest["manifest_sha256"] = hashlib.sha256(canonical.encode("utf-8")).hexdigest()
|
| 66 |
+
return manifest
|
app/tools/castp.py
CHANGED
|
@@ -71,60 +71,19 @@ async def _analyze_pockets(pdb_text: str, pdb_id: str, probe_radius: float, meth
|
|
| 71 |
if method_key in {"fpocket", "f_pocket"}:
|
| 72 |
fpocket = FPOCKET_BIN if Path(FPOCKET_BIN).exists() else (shutil.which("fpocket") or "")
|
| 73 |
if not fpocket or not Path(fpocket).exists():
|
| 74 |
-
|
| 75 |
-
|
| 76 |
-
|
| 77 |
-
|
| 78 |
-
method="fpocket",
|
| 79 |
-
status="FAILED",
|
| 80 |
-
engine_version="unavailable",
|
| 81 |
-
evidence_class="Unsupported/insufficient evidence",
|
| 82 |
-
reason="fpocket executable is unavailable; no heuristic was substituted",
|
| 83 |
-
)
|
| 84 |
-
_attach_structure_summary(pdb_text, result)
|
| 85 |
-
return result
|
| 86 |
-
methods_tried = [{"method": "fpocket", "status": "unavailable"}]
|
| 87 |
-
try:
|
| 88 |
-
loop = asyncio.get_running_loop()
|
| 89 |
-
result = await loop.run_in_executor(None, _analyze_pockets_sasa_sync, pdb_text, pdb_id, probe_radius)
|
| 90 |
-
except Exception as exc:
|
| 91 |
-
error = f"{type(exc).__name__}: {exc}"
|
| 92 |
-
methods_tried.append({"method": "sasa_heuristic", "status": "crashed"})
|
| 93 |
-
result = _empty_result(
|
| 94 |
-
pdb_id,
|
| 95 |
-
probe_radius,
|
| 96 |
-
method="none",
|
| 97 |
-
status="FAILED",
|
| 98 |
-
engine_version="fpocket unavailable; sasa_heuristic crashed",
|
| 99 |
-
evidence_class="Unsupported/insufficient evidence",
|
| 100 |
-
reason=error,
|
| 101 |
-
)
|
| 102 |
-
result["methods_tried"] = methods_tried
|
| 103 |
-
result["fallback_used"] = True
|
| 104 |
-
result["fallback_method"] = "sasa_heuristic"
|
| 105 |
-
result["error"] = str(exc)
|
| 106 |
-
_attach_structure_summary(pdb_text, result)
|
| 107 |
-
return result
|
| 108 |
-
heuristic_status = "ran_no_pockets" if not result.get("pockets") else "ran_with_pockets"
|
| 109 |
-
methods_tried.append({"method": "sasa_heuristic", "status": heuristic_status})
|
| 110 |
-
result["method"] = "sasa_heuristic"
|
| 111 |
-
result["engine"] = "BioNexus exploratory SASA heuristic"
|
| 112 |
-
result["status"] = "DEGRADED"
|
| 113 |
-
result["methods_tried"] = methods_tried
|
| 114 |
-
result["fallback_used"] = True
|
| 115 |
-
result["fallback_method"] = "sasa_heuristic"
|
| 116 |
-
result["note"] = (
|
| 117 |
-
"CASTp architecture: fpocket is unavailable in this deployment, so the "
|
| 118 |
-
"analysis fell back to the BioNexus exploratory SASA heuristic. All "
|
| 119 |
-
"values are attributed to sasa_heuristic; no fpocket or CASTp values "
|
| 120 |
-
"were substituted."
|
| 121 |
)
|
|
|
|
| 122 |
_attach_structure_summary(pdb_text, result)
|
| 123 |
return result
|
| 124 |
loop = asyncio.get_running_loop()
|
| 125 |
result = await loop.run_in_executor(None, _run_fpocket_analysis, fpocket, pdb_text, pdb_id, probe_radius)
|
| 126 |
result.setdefault("fallback_used", False)
|
| 127 |
-
result["methods_tried"] = [{"method": "fpocket", "status": "ok"}]
|
|
|
|
| 128 |
return result
|
| 129 |
|
| 130 |
if method_key in {"heuristic", "sasa", "sasa_heuristic", "bionexus_heuristic"}:
|
|
@@ -140,9 +99,10 @@ def _fpocket_version(fpocket_bin: str) -> str:
|
|
| 140 |
for args in ([fpocket_bin, "-v"], [fpocket_bin, "--version"]):
|
| 141 |
try:
|
| 142 |
completed = subprocess.run(args, capture_output=True, text=True, timeout=10)
|
| 143 |
-
text = (completed.stdout or completed.stderr or "")
|
| 144 |
-
|
| 145 |
-
|
|
|
|
| 146 |
except Exception:
|
| 147 |
continue
|
| 148 |
return "unknown"
|
|
@@ -155,7 +115,7 @@ def _run_fpocket_analysis(fpocket_bin: str, pdb_text: str, pdb_id: str, probe_ra
|
|
| 155 |
input_path.write_text(pdb_text, encoding="utf-8")
|
| 156 |
try:
|
| 157 |
completed = subprocess.run(
|
| 158 |
-
[fpocket_bin, "-f", str(input_path)
|
| 159 |
capture_output=True,
|
| 160 |
text=True,
|
| 161 |
timeout=60,
|
|
@@ -219,24 +179,24 @@ def _parse_fpocket(out_dir: Path, pdb_id: str, probe_radius: float, version: str
|
|
| 219 |
"alpha_spheres": None,
|
| 220 |
"num_residues": 0,
|
| 221 |
"residues": [],
|
| 222 |
-
"centroid":
|
| 223 |
}
|
| 224 |
continue
|
| 225 |
if current is None:
|
| 226 |
continue
|
| 227 |
-
|
| 228 |
value = _extract_number(line)
|
| 229 |
if value is None:
|
| 230 |
continue
|
| 231 |
-
if "druggability score"
|
| 232 |
current["druggability_score"] = value
|
| 233 |
-
elif
|
| 234 |
current["score"] = value
|
| 235 |
-
elif "volume"
|
| 236 |
current["volume"] = value
|
| 237 |
-
elif
|
| 238 |
current["area"] = value
|
| 239 |
-
elif "alpha
|
| 240 |
current["alpha_spheres"] = int(round(value))
|
| 241 |
if current:
|
| 242 |
pockets.append(current)
|
|
@@ -259,19 +219,21 @@ def _parse_fpocket(out_dir: Path, pdb_id: str, probe_radius: float, version: str
|
|
| 259 |
public_pockets: list[dict[str, Any]] = []
|
| 260 |
for pocket in pockets:
|
| 261 |
volume = pocket.get("volume")
|
| 262 |
-
radius = round((3 * volume / (4 * math.pi)) ** (1 / 3), 2) if isinstance(volume, (int, float)) and volume > 0 else
|
| 263 |
public_pockets.append({
|
| 264 |
"id": pocket["id"],
|
| 265 |
-
"area_sa": round(float(pocket["area"]), 3) if pocket.get("area") is not None else
|
| 266 |
-
"volume_sa": round(float(volume), 3) if volume is not None else
|
| 267 |
"score": pocket.get("score"),
|
| 268 |
"druggability_score": pocket.get("druggability_score"),
|
| 269 |
"alpha_spheres": pocket.get("alpha_spheres"),
|
| 270 |
"num_residues": int(pocket.get("num_residues") or 0),
|
| 271 |
"residues": pocket.get("residues", []),
|
| 272 |
-
"centroid": [round(float(value), 3) for value in pocket.get("centroid"
|
| 273 |
"radius": radius,
|
| 274 |
"method_metrics": {
|
|
|
|
|
|
|
| 275 |
"fpocket_score": pocket.get("score"),
|
| 276 |
"fpocket_druggability_score": pocket.get("druggability_score"),
|
| 277 |
"fpocket_alpha_spheres": pocket.get("alpha_spheres"),
|
|
@@ -287,7 +249,8 @@ def _parse_fpocket(out_dir: Path, pdb_id: str, probe_radius: float, version: str
|
|
| 287 |
"evidence_class": "Deterministic computation",
|
| 288 |
"validation": {"source": "fpocket native output", "castp_values_substituted": False},
|
| 289 |
"pdb_id": pdb_id,
|
| 290 |
-
"probe_radius":
|
|
|
|
| 291 |
"total_residues": 0,
|
| 292 |
"pockets": public_pockets,
|
| 293 |
}
|
|
|
|
| 71 |
if method_key in {"fpocket", "f_pocket"}:
|
| 72 |
fpocket = FPOCKET_BIN if Path(FPOCKET_BIN).exists() else (shutil.which("fpocket") or "")
|
| 73 |
if not fpocket or not Path(fpocket).exists():
|
| 74 |
+
result = _empty_result(
|
| 75 |
+
pdb_id, probe_radius, method="fpocket", status="FAILED",
|
| 76 |
+
engine_version="unavailable", evidence_class="Unsupported/insufficient evidence",
|
| 77 |
+
reason="fpocket executable is unavailable; no heuristic substitution was performed. Select the exploratory SASA heuristic explicitly to run that method.",
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
| 78 |
)
|
| 79 |
+
result["methods_tried"] = [{"method": "fpocket", "status": "unavailable"}]
|
| 80 |
_attach_structure_summary(pdb_text, result)
|
| 81 |
return result
|
| 82 |
loop = asyncio.get_running_loop()
|
| 83 |
result = await loop.run_in_executor(None, _run_fpocket_analysis, fpocket, pdb_text, pdb_id, probe_radius)
|
| 84 |
result.setdefault("fallback_used", False)
|
| 85 |
+
result["methods_tried"] = [{"method": "fpocket", "status": "ok" if result.get("status", "VALID") == "VALID" else "failed"}]
|
| 86 |
+
_attach_structure_summary(pdb_text, result)
|
| 87 |
return result
|
| 88 |
|
| 89 |
if method_key in {"heuristic", "sasa", "sasa_heuristic", "bionexus_heuristic"}:
|
|
|
|
| 99 |
for args in ([fpocket_bin, "-v"], [fpocket_bin, "--version"]):
|
| 100 |
try:
|
| 101 |
completed = subprocess.run(args, capture_output=True, text=True, timeout=10)
|
| 102 |
+
text = re.sub(r"\x1b\[[0-9;]*m", "", (completed.stdout or "") + "\n" + (completed.stderr or ""))
|
| 103 |
+
match = re.search(r"\bfpocket\s+(?:version\s+)?(\d+(?:\.\d+)+(?:[-+][\w.]+)?)", text, re.IGNORECASE)
|
| 104 |
+
if match:
|
| 105 |
+
return match.group(1)
|
| 106 |
except Exception:
|
| 107 |
continue
|
| 108 |
return "unknown"
|
|
|
|
| 115 |
input_path.write_text(pdb_text, encoding="utf-8")
|
| 116 |
try:
|
| 117 |
completed = subprocess.run(
|
| 118 |
+
[fpocket_bin, "-f", str(input_path)],
|
| 119 |
capture_output=True,
|
| 120 |
text=True,
|
| 121 |
timeout=60,
|
|
|
|
| 179 |
"alpha_spheres": None,
|
| 180 |
"num_residues": 0,
|
| 181 |
"residues": [],
|
| 182 |
+
"centroid": None,
|
| 183 |
}
|
| 184 |
continue
|
| 185 |
if current is None:
|
| 186 |
continue
|
| 187 |
+
descriptor = line.split(":", 1)[0].strip().lower()
|
| 188 |
value = _extract_number(line)
|
| 189 |
if value is None:
|
| 190 |
continue
|
| 191 |
+
if descriptor == "druggability score":
|
| 192 |
current["druggability_score"] = value
|
| 193 |
+
elif descriptor in {"score", "pocket score"}:
|
| 194 |
current["score"] = value
|
| 195 |
+
elif descriptor == "volume":
|
| 196 |
current["volume"] = value
|
| 197 |
+
elif descriptor in {"total sasa", "surface area", "area"}:
|
| 198 |
current["area"] = value
|
| 199 |
+
elif descriptor in {"number of alpha spheres", "alpha spheres"}:
|
| 200 |
current["alpha_spheres"] = int(round(value))
|
| 201 |
if current:
|
| 202 |
pockets.append(current)
|
|
|
|
| 219 |
public_pockets: list[dict[str, Any]] = []
|
| 220 |
for pocket in pockets:
|
| 221 |
volume = pocket.get("volume")
|
| 222 |
+
radius = round((3 * volume / (4 * math.pi)) ** (1 / 3), 2) if isinstance(volume, (int, float)) and volume > 0 else None
|
| 223 |
public_pockets.append({
|
| 224 |
"id": pocket["id"],
|
| 225 |
+
"area_sa": round(float(pocket["area"]), 3) if pocket.get("area") is not None else None,
|
| 226 |
+
"volume_sa": round(float(volume), 3) if volume is not None else None,
|
| 227 |
"score": pocket.get("score"),
|
| 228 |
"druggability_score": pocket.get("druggability_score"),
|
| 229 |
"alpha_spheres": pocket.get("alpha_spheres"),
|
| 230 |
"num_residues": int(pocket.get("num_residues") or 0),
|
| 231 |
"residues": pocket.get("residues", []),
|
| 232 |
+
"centroid": [round(float(value), 3) for value in pocket["centroid"]] if pocket.get("centroid") else None,
|
| 233 |
"radius": radius,
|
| 234 |
"method_metrics": {
|
| 235 |
+
"fpocket_total_sasa": pocket.get("area"),
|
| 236 |
+
"fpocket_volume": pocket.get("volume"),
|
| 237 |
"fpocket_score": pocket.get("score"),
|
| 238 |
"fpocket_druggability_score": pocket.get("druggability_score"),
|
| 239 |
"fpocket_alpha_spheres": pocket.get("alpha_spheres"),
|
|
|
|
| 249 |
"evidence_class": "Deterministic computation",
|
| 250 |
"validation": {"source": "fpocket native output", "castp_values_substituted": False},
|
| 251 |
"pdb_id": pdb_id,
|
| 252 |
+
"probe_radius": None,
|
| 253 |
+
"requested_probe_radius": probe_radius,
|
| 254 |
"total_residues": 0,
|
| 255 |
"pockets": public_pockets,
|
| 256 |
}
|
app/tools/function_predict.py
CHANGED
|
@@ -236,7 +236,7 @@ def _interpro_to_go(hits: list[dict]) -> list[dict]:
|
|
| 236 |
accessions: dict[str, int] = {}
|
| 237 |
for hit in hits:
|
| 238 |
accession = hit.get("accession")
|
| 239 |
-
if accession:
|
| 240 |
accessions[accession] = accessions.get(accession, 0) + 1
|
| 241 |
|
| 242 |
aggregated: dict[str, dict] = {}
|
|
@@ -356,6 +356,8 @@ def _prediction(sequence: str, pdb_id: str, domain_hits: list[dict],
|
|
| 356 |
retrieval_failures += [f for f in _RETRIEVAL_FAILURES if f not in retrieval_failures]
|
| 357 |
|
| 358 |
go_terms = _interpro_to_go(domain_hits) if domain_hits else []
|
|
|
|
|
|
|
| 359 |
|
| 360 |
if go_terms:
|
| 361 |
status = "inferred"
|
|
|
|
| 236 |
accessions: dict[str, int] = {}
|
| 237 |
for hit in hits:
|
| 238 |
accession = hit.get("accession")
|
| 239 |
+
if accession and accession.startswith("IPR"):
|
| 240 |
accessions[accession] = accessions.get(accession, 0) + 1
|
| 241 |
|
| 242 |
aggregated: dict[str, dict] = {}
|
|
|
|
| 356 |
retrieval_failures += [f for f in _RETRIEVAL_FAILURES if f not in retrieval_failures]
|
| 357 |
|
| 358 |
go_terms = _interpro_to_go(domain_hits) if domain_hits else []
|
| 359 |
+
# Mapping calls can fail after the scan succeeds. Preserve those failures too.
|
| 360 |
+
retrieval_failures += [f for f in _RETRIEVAL_FAILURES if f not in retrieval_failures]
|
| 361 |
|
| 362 |
if go_terms:
|
| 363 |
status = "inferred"
|
app/tools/oligo_qc.py
CHANGED
|
@@ -10,9 +10,8 @@ the canonical primer design workflow):
|
|
| 10 |
* In-silico PCR: where each primer binds in the template, expected amplicon
|
| 11 |
length, and specificity (number of binding sites).
|
| 12 |
|
| 13 |
-
|
| 14 |
-
|
| 15 |
-
ranked; it is not a substitute for a full NN thermodynamics calculation.
|
| 16 |
"""
|
| 17 |
|
| 18 |
from __future__ import annotations
|
|
@@ -26,15 +25,6 @@ _COMP = {ord('A'): 'T', ord('T'): 'A', ord('G'): 'C', ord('C'): 'G',
|
|
| 26 |
ord('a'): 't', ord('t'): 'a', ord('g'): 'c', ord('c'): 'g',
|
| 27 |
ord('n'): 'n', ord('u'): 'a'}
|
| 28 |
|
| 29 |
-
# Per-base-pair stacking energies (kcal/mol) for the complementary stem.
|
| 30 |
-
# Only canonical Watson-Crick pairs count towards a stem; wobble pairs are
|
| 31 |
-
# ignored so short weak "stems" don't turn every oligo into a dimer/hairpin.
|
| 32 |
-
_PAIR_DG = {
|
| 33 |
-
('A', 'T'): -1.0, ('T', 'A'): -1.0,
|
| 34 |
-
('G', 'C'): -2.0, ('C', 'G'): -2.0,
|
| 35 |
-
}
|
| 36 |
-
|
| 37 |
-
|
| 38 |
def clean(seq: str) -> str:
|
| 39 |
"""Upper-case and strip whitespace (and any non-letter junk)."""
|
| 40 |
return re.sub(r"[^A-Za-z]", "", seq).upper().replace("U", "T")
|
|
@@ -56,133 +46,69 @@ def gc_content(seq: str) -> float:
|
|
| 56 |
return round((s.count("G") + s.count("C")) / len(s) * 100.0, 1)
|
| 57 |
|
| 58 |
|
| 59 |
-
|
| 60 |
-
return _PAIR_DG.get((a, b), _PAIR_DG.get((b, a), 0.0))
|
| 61 |
|
| 62 |
|
| 63 |
-
def
|
| 64 |
-
|
|
|
|
|
|
|
|
|
|
|
|
|
| 65 |
|
| 66 |
-
Uses the classic NN-approximation formula Tm = 64.9 + 41*(nGC - 16.4)/N
|
| 67 |
-
(valid at ~50 mM monovalent salt) plus the SantaLucia salt correction.
|
| 68 |
-
"""
|
| 69 |
-
s = clean(seq)
|
| 70 |
-
n = len(s)
|
| 71 |
-
if n < 1:
|
| 72 |
-
return 0.0
|
| 73 |
-
n_gc = s.count("G") + s.count("C")
|
| 74 |
-
tm = 64.9 + 41.0 * (n_gc - 16.4) / n
|
| 75 |
-
if salt_mm != 50.0:
|
| 76 |
-
tm += 16.6 * math.log10(max(salt_mm, 0.01) / 50.0)
|
| 77 |
-
return round(tm, 1)
|
| 78 |
-
|
| 79 |
-
|
| 80 |
-
def _risk_from_dg(dg: float) -> Tuple[str, str]:
|
| 81 |
-
if dg <= -7.0:
|
| 82 |
-
return "high", "Strongly stabilizing — high risk of mis-priming."
|
| 83 |
-
if dg <= -4.5:
|
| 84 |
-
return "medium", "Moderately stabilizing — keep an eye on it."
|
| 85 |
-
if dg <= -3.0:
|
| 86 |
-
return "low", "Weakly stabilizing — unlikely to matter in practice."
|
| 87 |
-
return "none", "No significant structure predicted."
|
| 88 |
-
|
| 89 |
-
|
| 90 |
-
def hairpin_analysis(seq: str, min_stem: int = 4, min_loop: int = 3) -> Dict:
|
| 91 |
-
"""Predict the strongest intra-molecular hairpin.
|
| 92 |
-
|
| 93 |
-
A hairpin forms when the strand folds so that one arm is the reverse
|
| 94 |
-
complement of the other arm, separated by a loop of at least `min_loop`
|
| 95 |
-
bases: seq[i:i+L] must equal reverse_complement(seq[j:j+L]). Returns
|
| 96 |
-
delta-G, stem length, stem sequence, loop sequence and a risk level.
|
| 97 |
-
"""
|
| 98 |
-
s = clean(seq)
|
| 99 |
-
n = len(s)
|
| 100 |
-
best: Optional[Tuple[float, int, int, int]] = None # (dg, stem_len, i, j)
|
| 101 |
-
for i in range(n):
|
| 102 |
-
for j in range(i + min_stem + min_loop, n):
|
| 103 |
-
# loop = j - (i + L) must stay >= min_loop
|
| 104 |
-
lmax = min(j - i - min_loop, n - j, 10)
|
| 105 |
-
if lmax < min_stem:
|
| 106 |
-
continue
|
| 107 |
-
for L in range(min_stem, lmax + 1):
|
| 108 |
-
dg = 0.0
|
| 109 |
-
ok = True
|
| 110 |
-
for k in range(L):
|
| 111 |
-
e = _pair_energy(s[i + k], s[j + L - 1 - k])
|
| 112 |
-
if e == 0.0:
|
| 113 |
-
ok = False
|
| 114 |
-
break
|
| 115 |
-
dg += e
|
| 116 |
-
if ok and (best is None or dg < best[0]):
|
| 117 |
-
best = (dg, L, i, j)
|
| 118 |
-
|
| 119 |
-
if best is None:
|
| 120 |
-
return {"dg": 5.0, "stem_length": 0, "stem": "", "loop": "",
|
| 121 |
-
"risk": "none", "note": "No hairpin structure predicted."}
|
| 122 |
-
dg, stem_len, i, j = best
|
| 123 |
-
stem = s[i:i + stem_len]
|
| 124 |
-
loop = s[i + stem_len:j]
|
| 125 |
-
risk, note = _risk_from_dg(dg)
|
| 126 |
-
return {
|
| 127 |
-
"dg": round(dg, 1),
|
| 128 |
-
"stem_length": stem_len,
|
| 129 |
-
"stem": stem,
|
| 130 |
-
"loop": loop,
|
| 131 |
-
"risk": risk,
|
| 132 |
-
"note": note,
|
| 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 |
-
return
|
| 178 |
-
|
| 179 |
-
|
| 180 |
-
|
| 181 |
-
|
| 182 |
-
|
| 183 |
-
|
| 184 |
-
|
| 185 |
-
|
|
|
|
| 186 |
|
| 187 |
|
| 188 |
def find_binding_sites(template: str, primer: str) -> List[int]:
|
|
@@ -262,17 +188,13 @@ def in_silico_pcr(template: str, left: str, right: str,
|
|
| 262 |
return result
|
| 263 |
|
| 264 |
|
| 265 |
-
def oligo_report(seq: str) -> Dict:
|
| 266 |
-
"""
|
| 267 |
-
s =
|
| 268 |
-
|
| 269 |
-
|
| 270 |
-
return {
|
| 271 |
-
|
| 272 |
-
|
| 273 |
-
|
| 274 |
-
|
| 275 |
-
"tm_model": "IDT-style nearest-neighbor approximation (64.9 + 41*(nGC-16.4)/N) with SantaLucia salt correction at 50 mM; not a fully calibrated NN melt calculation",
|
| 276 |
-
"hairpin": hp,
|
| 277 |
-
"self_dimer": sd,
|
| 278 |
-
}
|
|
|
|
| 10 |
* In-silico PCR: where each primer binds in the template, expected amplicon
|
| 11 |
length, and specificity (number of binding sites).
|
| 12 |
|
| 13 |
+
Tm and secondary structures use primer3-py thermodynamics with explicit
|
| 14 |
+
reaction conditions. Risk thresholds remain screening rules, not PCR validation.
|
|
|
|
| 15 |
"""
|
| 16 |
|
| 17 |
from __future__ import annotations
|
|
|
|
| 25 |
ord('a'): 't', ord('t'): 'a', ord('g'): 'c', ord('c'): 'g',
|
| 26 |
ord('n'): 'n', ord('u'): 'a'}
|
| 27 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
| 28 |
def clean(seq: str) -> str:
|
| 29 |
"""Upper-case and strip whitespace (and any non-letter junk)."""
|
| 30 |
return re.sub(r"[^A-Za-z]", "", seq).upper().replace("U", "T")
|
|
|
|
| 46 |
return round((s.count("G") + s.count("C")) / len(s) * 100.0, 1)
|
| 47 |
|
| 48 |
|
| 49 |
+
DEFAULT_CONDITIONS = {"mv_conc": 50.0, "dv_conc": 1.5, "dntp_conc": 0.6, "dna_conc": 50.0, "temp_c": 37.0}
|
|
|
|
| 50 |
|
| 51 |
|
| 52 |
+
def _engine():
|
| 53 |
+
try:
|
| 54 |
+
import primer3
|
| 55 |
+
return primer3
|
| 56 |
+
except ImportError as exc:
|
| 57 |
+
raise RuntimeError("Primer3 thermodynamics unavailable; install primer3-py. No heuristic fallback is used.") from exc
|
| 58 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
| 59 |
|
| 60 |
+
def _sequence(seq: str) -> str:
|
| 61 |
+
s = clean(seq)
|
| 62 |
+
if not s or set(s) - set("ACGT") or len(s) > 60:
|
| 63 |
+
raise ValueError("Thermodynamic QC requires 1–60 unambiguous DNA bases (A/C/G/T); ambiguous bases are not scored.")
|
| 64 |
+
return s
|
| 65 |
+
|
| 66 |
+
|
| 67 |
+
def _conditions(conditions=None):
|
| 68 |
+
values = {**DEFAULT_CONDITIONS, **(conditions or {})}
|
| 69 |
+
if set(values) != set(DEFAULT_CONDITIONS):
|
| 70 |
+
raise ValueError("Unknown thermodynamic condition")
|
| 71 |
+
if any(not math.isfinite(float(v)) or float(v) < 0 for k, v in values.items() if k != "temp_c"):
|
| 72 |
+
raise ValueError("Concentrations must be finite and non-negative")
|
| 73 |
+
if values["dna_conc"] <= 0 or not -20 <= values["temp_c"] <= 100:
|
| 74 |
+
raise ValueError("DNA concentration must be positive and temperature must be -20 to 100 C")
|
| 75 |
+
return values
|
| 76 |
+
|
| 77 |
+
|
| 78 |
+
def salt_adjusted_tm(seq: str, salt_mm: float = 50.0, conditions=None) -> float:
|
| 79 |
+
"""SantaLucia nearest-neighbour Tm from Primer3; concentrations in mM except DNA in nM."""
|
| 80 |
+
c = _conditions(conditions)
|
| 81 |
+
if conditions is None:
|
| 82 |
+
c["mv_conc"] = salt_mm
|
| 83 |
+
kwargs = {k: v for k, v in c.items() if k != "temp_c"}
|
| 84 |
+
return round(_engine().calc_tm(_sequence(seq), **kwargs, tm_method="santalucia", salt_corrections_method="santalucia"), 2)
|
| 85 |
+
|
| 86 |
+
|
| 87 |
+
def _structure(result, conditions):
|
| 88 |
+
# Primer3 returns free energy in cal/mol. Public output uses kcal/mol.
|
| 89 |
+
dg = float(result.dg) / 1000.0
|
| 90 |
+
found = bool(result.structure_found)
|
| 91 |
+
risk = "none" if not found or dg > -3 else "low" if dg > -4.5 else "medium" if dg > -7 else "high"
|
| 92 |
+
return {"dg": round(dg, 3), "dg_units": "kcal/mol", "tm_c": float(result.tm),
|
| 93 |
+
"structure_found": found, "structure": getattr(result, "ascii_structure", "") or "",
|
| 94 |
+
"risk": risk, "risk_model": "Screening thresholds on predicted delta-G: low <= -3, medium <= -4.5, high <= -7 kcal/mol; not experimental PCR validation",
|
| 95 |
+
"conditions": conditions, "engine": "primer3-py", "engine_version": getattr(_engine(), "__version__", "unknown"),
|
| 96 |
+
"note": "Primer3 thermodynamic prediction; PCR efficiency and genome-wide specificity are not established."}
|
| 97 |
+
|
| 98 |
+
|
| 99 |
+
def hairpin_analysis(seq: str, min_stem: int = 4, min_loop: int = 3, conditions=None) -> Dict:
|
| 100 |
+
"""Primer3 hairpin thermodynamics (legacy stem arguments do not alter its model)."""
|
| 101 |
+
c = _conditions(conditions)
|
| 102 |
+
return _structure(_engine().calc_hairpin(_sequence(seq), **c, output_structure=True), c)
|
| 103 |
+
|
| 104 |
+
|
| 105 |
+
def dimer_analysis(a: str, b: str, min_stem: int = 4, conditions=None) -> Dict:
|
| 106 |
+
"""Primer3 homo/heterodimer thermodynamics, without invented stem coordinates."""
|
| 107 |
+
c = _conditions(conditions)
|
| 108 |
+
pa, pb = _sequence(a), _sequence(b)
|
| 109 |
+
engine = _engine()
|
| 110 |
+
result = engine.calc_homodimer(pa, **c, output_structure=True) if pa == pb else engine.calc_heterodimer(pa, pb, **c, output_structure=True)
|
| 111 |
+
return _structure(result, c)
|
| 112 |
|
| 113 |
|
| 114 |
def find_binding_sites(template: str, primer: str) -> List[int]:
|
|
|
|
| 188 |
return result
|
| 189 |
|
| 190 |
|
| 191 |
+
def oligo_report(seq: str, conditions=None) -> Dict:
|
| 192 |
+
"""Nearest-neighbour thermodynamic QC with fully recorded reaction conditions."""
|
| 193 |
+
s = _sequence(seq)
|
| 194 |
+
c = _conditions(conditions)
|
| 195 |
+
tm = salt_adjusted_tm(s, conditions=c)
|
| 196 |
+
return {"sequence": s, "length": len(s), "gc": gc_content(s), "tm_c": tm,
|
| 197 |
+
"tm_50mM": tm if c["mv_conc"] == 50 else None,
|
| 198 |
+
"tm_model": "Primer3 SantaLucia nearest-neighbor with SantaLucia salt correction",
|
| 199 |
+
"conditions": c, "concentration_units": {"mv_conc": "mM", "dv_conc": "mM", "dntp_conc": "mM", "dna_conc": "nM", "temp_c": "C"},
|
| 200 |
+
"hairpin": hairpin_analysis(s, conditions=c), "self_dimer": dimer_analysis(s, s, conditions=c)}
|
|
|
|
|
|
|
|
|
|
|
|
app/tools/primer_advanced.py
CHANGED
|
@@ -59,7 +59,7 @@ def snp_overlap(
|
|
| 59 |
}
|
| 60 |
|
| 61 |
|
| 62 |
-
def multiplex_compatibility(pairs: list[dict[str, Any]], max_tm_spread: float = 3.0) -> dict:
|
| 63 |
"""Screen primer pairs for multiplex use using Tm spread and cross-dimers.
|
| 64 |
|
| 65 |
This is a deterministic screening heuristic. It does not model full PCR
|
|
@@ -73,8 +73,8 @@ def multiplex_compatibility(pairs: list[dict[str, Any]], max_tm_spread: float =
|
|
| 73 |
if not left or not right:
|
| 74 |
raise ValueError(f"pair {idx} requires left_seq and right_seq")
|
| 75 |
pid = str(pair.get("id") or f"pair-{idx+1}")
|
| 76 |
-
left_tm =
|
| 77 |
-
right_tm =
|
| 78 |
normalized.append({"id": pid, "left_seq": left, "right_seq": right, "left_tm": left_tm, "right_tm": right_tm})
|
| 79 |
all_primers.extend([(f"{pid}:left", left), (f"{pid}:right", right)])
|
| 80 |
|
|
@@ -85,7 +85,7 @@ def multiplex_compatibility(pairs: list[dict[str, Any]], max_tm_spread: float =
|
|
| 85 |
for (name_a, seq_a), (name_b, seq_b) in combinations(all_primers, 2):
|
| 86 |
if name_a.split(":")[0] == name_b.split(":")[0]:
|
| 87 |
continue
|
| 88 |
-
dimer = dimer_analysis(seq_a, seq_b)
|
| 89 |
risk = dimer.get("risk", "none")
|
| 90 |
if risk in {"high", "medium"}:
|
| 91 |
high_risk += 1
|
|
@@ -94,7 +94,8 @@ def multiplex_compatibility(pairs: list[dict[str, Any]], max_tm_spread: float =
|
|
| 94 |
"primer_b": name_b,
|
| 95 |
"dg": dimer.get("dg"),
|
| 96 |
"risk": risk,
|
| 97 |
-
"
|
|
|
|
| 98 |
})
|
| 99 |
|
| 100 |
compatible = tm_spread <= max_tm_spread and high_risk == 0
|
|
|
|
| 59 |
}
|
| 60 |
|
| 61 |
|
| 62 |
+
def multiplex_compatibility(pairs: list[dict[str, Any]], max_tm_spread: float = 3.0, conditions=None) -> dict:
|
| 63 |
"""Screen primer pairs for multiplex use using Tm spread and cross-dimers.
|
| 64 |
|
| 65 |
This is a deterministic screening heuristic. It does not model full PCR
|
|
|
|
| 73 |
if not left or not right:
|
| 74 |
raise ValueError(f"pair {idx} requires left_seq and right_seq")
|
| 75 |
pid = str(pair.get("id") or f"pair-{idx+1}")
|
| 76 |
+
left_tm = salt_adjusted_tm(left, conditions=conditions)
|
| 77 |
+
right_tm = salt_adjusted_tm(right, conditions=conditions)
|
| 78 |
normalized.append({"id": pid, "left_seq": left, "right_seq": right, "left_tm": left_tm, "right_tm": right_tm})
|
| 79 |
all_primers.extend([(f"{pid}:left", left), (f"{pid}:right", right)])
|
| 80 |
|
|
|
|
| 85 |
for (name_a, seq_a), (name_b, seq_b) in combinations(all_primers, 2):
|
| 86 |
if name_a.split(":")[0] == name_b.split(":")[0]:
|
| 87 |
continue
|
| 88 |
+
dimer = dimer_analysis(seq_a, seq_b, conditions=conditions)
|
| 89 |
risk = dimer.get("risk", "none")
|
| 90 |
if risk in {"high", "medium"}:
|
| 91 |
high_risk += 1
|
|
|
|
| 94 |
"primer_b": name_b,
|
| 95 |
"dg": dimer.get("dg"),
|
| 96 |
"risk": risk,
|
| 97 |
+
"structure": dimer.get("structure"),
|
| 98 |
+
"conditions": dimer.get("conditions"),
|
| 99 |
})
|
| 100 |
|
| 101 |
compatible = tm_spread <= max_tm_spread and high_risk == 0
|
app/tools/restriction_map.py
ADDED
|
@@ -0,0 +1,54 @@
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
| 1 |
+
"""Bounded exact-DNA restriction map; enzyme definitions come from Biopython."""
|
| 2 |
+
from __future__ import annotations
|
| 3 |
+
|
| 4 |
+
import Bio
|
| 5 |
+
from Bio import Restriction
|
| 6 |
+
from app.tools.sequence_utilities import SequenceUtilitiesError, clean_sequence, reverse_complement_sequence
|
| 7 |
+
|
| 8 |
+
ENZYME_NAMES = ("EcoRI", "BamHI", "HindIII", "SalI", "XbaI", "XhoI", "NotI", "KpnI", "SmaI", "PstI", "SacI", "BsaI", "BsmBI")
|
| 9 |
+
|
| 10 |
+
|
| 11 |
+
def restriction_map(sequence: str, enzymes: list[str], circular: bool = False) -> dict:
|
| 12 |
+
seq = clean_sequence(sequence, "dna")
|
| 13 |
+
if set(seq) - set("ACGT"):
|
| 14 |
+
raise SequenceUtilitiesError("Restriction maps require unambiguous A/C/G/T DNA; ambiguous bases cannot establish definite sites")
|
| 15 |
+
if len(seq) > 200000:
|
| 16 |
+
raise SequenceUtilitiesError("Restriction map input is limited to 200,000 bases")
|
| 17 |
+
selected = list(dict.fromkeys(enzymes))
|
| 18 |
+
if not selected or len(selected) > len(ENZYME_NAMES) or any(e not in ENZYME_NAMES for e in selected):
|
| 19 |
+
raise SequenceUtilitiesError("Choose supported restriction enzymes: " + ", ".join(ENZYME_NAMES))
|
| 20 |
+
length = len(seq)
|
| 21 |
+
hits, definitions, cuts = [], [], set()
|
| 22 |
+
for name in selected:
|
| 23 |
+
enzyme = getattr(Restriction, name)
|
| 24 |
+
motif = enzyme.site
|
| 25 |
+
definitions.append({"name": name, "recognition": motif, "forward_top_offset": enzyme.fst5, "forward_bottom_offset": len(motif) + enzyme.fst3})
|
| 26 |
+
if len(motif) > length:
|
| 27 |
+
continue
|
| 28 |
+
scan = seq + seq[:len(motif) - 1] if circular else seq
|
| 29 |
+
orientations = [("+", motif)]
|
| 30 |
+
rc = reverse_complement_sequence(motif, "dna")
|
| 31 |
+
if rc != motif:
|
| 32 |
+
orientations.append(("-", rc))
|
| 33 |
+
for strand, recognition in orientations:
|
| 34 |
+
for start in range(length if circular else length - len(motif) + 1):
|
| 35 |
+
if scan[start:start + len(motif)] != recognition:
|
| 36 |
+
continue
|
| 37 |
+
top = start + enzyme.fst5 if strand == "+" else start - enzyme.fst3
|
| 38 |
+
bottom = start + len(motif) + enzyme.fst3 if strand == "+" else start + len(motif) - enzyme.fst5
|
| 39 |
+
complete = circular or (0 <= top <= length and 0 <= bottom <= length)
|
| 40 |
+
if circular:
|
| 41 |
+
top, bottom = top % length, bottom % length
|
| 42 |
+
if complete:
|
| 43 |
+
cuts.add(top)
|
| 44 |
+
hits.append({"enzyme": name, "recognition": recognition, "start": start + 1, "end": ((start + len(motif) - 1) % length) + 1, "strand": strand, "spans_origin": start + len(motif) > length, "top_cut_boundary": top, "bottom_cut_boundary": bottom, "complete_cut": complete})
|
| 45 |
+
ordered = sorted(cuts)
|
| 46 |
+
if circular and ordered:
|
| 47 |
+
fragments = [((ordered[(i + 1) % len(ordered)] - cut) % length) or length for i, cut in enumerate(ordered)]
|
| 48 |
+
elif circular:
|
| 49 |
+
fragments = [length]
|
| 50 |
+
else:
|
| 51 |
+
boundaries = sorted({0, length, *ordered})
|
| 52 |
+
fragments = [b - a for a, b in zip(boundaries, boundaries[1:])]
|
| 53 |
+
hits.sort(key=lambda h: (h["start"], h["enzyme"], h["strand"]))
|
| 54 |
+
return {"length": length, "topology": "circular" if circular else "linear", "enzymes": definitions, "sites": hits, "cut_boundaries": ordered, "fragment_lengths": fragments, "fragment_total": sum(fragments), "product_topology": "circular" if circular and not ordered else "linear", "coordinate_convention": "Recognition start/end: 1-based inclusive. Cleavage boundaries: number of bases before cut on displayed forward strand; 0 is before base 1. Circular boundaries modulo length.", "digest_assumptions": "Complete digestion at supported recognition sites; ignores methylation, star activity, reaction conditions and enzyme efficiency. Fragment lengths use forward-strand cleavage boundaries; overhang bases on the opposite strand may differ.", "provenance": {"source": "Biopython Bio.Restriction curated enzyme definitions", "version": Bio.__version__, "source_url": "https://biopython.org/docs/latest/api/Bio.Restriction.html"}}
|
app/tools/structure_prep.py
CHANGED
|
@@ -263,11 +263,13 @@ def run_fpocket(pdb_text: str, probe_radius: float = 1.4) -> FpocketResult:
|
|
| 263 |
|
| 264 |
try:
|
| 265 |
result = subprocess.run(
|
| 266 |
-
[FPOCKET_BIN, "-f", str(in_path)
|
| 267 |
capture_output=True,
|
| 268 |
text=True,
|
| 269 |
timeout=60,
|
| 270 |
)
|
|
|
|
|
|
|
| 271 |
fpocket_out = Path(tmpdir) / "input_out"
|
| 272 |
return _parse_fpocket_output(fpocket_out, result.stdout + result.stderr)
|
| 273 |
except subprocess.TimeoutExpired:
|
|
@@ -278,61 +280,15 @@ def run_fpocket(pdb_text: str, probe_radius: float = 1.4) -> FpocketResult:
|
|
| 278 |
|
| 279 |
def _parse_fpocket_output(out_dir: Path, raw_output: str) -> FpocketResult:
|
| 280 |
"""Parse fpocket output directory for pocket information."""
|
| 281 |
-
|
| 282 |
-
|
| 283 |
-
|
| 284 |
-
|
| 285 |
-
|
| 286 |
-
|
| 287 |
-
|
| 288 |
-
|
| 289 |
-
|
| 290 |
-
current_pocket: dict[str, Any] = {}
|
| 291 |
-
|
| 292 |
-
for line in text.splitlines():
|
| 293 |
-
line = line.strip()
|
| 294 |
-
if line.startswith("Pocket"):
|
| 295 |
-
if current_pocket:
|
| 296 |
-
pockets.append(current_pocket)
|
| 297 |
-
pocket_id_match = re.search(r"Pocket\s+(\d+)", line)
|
| 298 |
-
current_pocket = {
|
| 299 |
-
"id": int(pocket_id_match.group(1)) if pocket_id_match else len(pockets) + 1,
|
| 300 |
-
"druggability_score": 0.0,
|
| 301 |
-
"volume": 0.0,
|
| 302 |
-
"area": 0.0,
|
| 303 |
-
"score": 0.0,
|
| 304 |
-
"num_residues": 0,
|
| 305 |
-
}
|
| 306 |
-
elif "Druggability Score" in line:
|
| 307 |
-
m = re.search(r":\s*([\d.]+)", line)
|
| 308 |
-
if m:
|
| 309 |
-
current_pocket["druggability_score"] = float(m.group(1))
|
| 310 |
-
elif "Volume" in line:
|
| 311 |
-
m = re.search(r":\s*([\d.]+)", line)
|
| 312 |
-
if m:
|
| 313 |
-
current_pocket["volume"] = float(m.group(1))
|
| 314 |
-
elif "Area" in line:
|
| 315 |
-
m = re.search(r":\s*([\d.]+)", line)
|
| 316 |
-
if m:
|
| 317 |
-
current_pocket["area"] = float(m.group(1))
|
| 318 |
-
elif "Score" in line and "Drug" not in line:
|
| 319 |
-
m = re.search(r":\s*([\d.]+)", line)
|
| 320 |
-
if m:
|
| 321 |
-
current_pocket["score"] = float(m.group(1))
|
| 322 |
-
elif "Number of residues" in line:
|
| 323 |
-
m = re.search(r":\s*(\d+)", line)
|
| 324 |
-
if m:
|
| 325 |
-
current_pocket["num_residues"] = int(m.group(1))
|
| 326 |
-
|
| 327 |
-
if current_pocket:
|
| 328 |
-
pockets.append(current_pocket)
|
| 329 |
-
|
| 330 |
-
result.pockets = pockets
|
| 331 |
-
result.pocket_count = len(pockets)
|
| 332 |
-
except Exception as e:
|
| 333 |
-
logger.warning("Failed to parse fpocket output: %s", e)
|
| 334 |
-
|
| 335 |
-
return result
|
| 336 |
|
| 337 |
|
| 338 |
# ── Step 4b: SASA concave-packing fallback (no fpocket / scipy needed) ──────
|
|
|
|
| 263 |
|
| 264 |
try:
|
| 265 |
result = subprocess.run(
|
| 266 |
+
[FPOCKET_BIN, "-f", str(in_path)],
|
| 267 |
capture_output=True,
|
| 268 |
text=True,
|
| 269 |
timeout=60,
|
| 270 |
)
|
| 271 |
+
if result.returncode != 0:
|
| 272 |
+
return FpocketResult(raw_output=result.stdout + result.stderr, status="error")
|
| 273 |
fpocket_out = Path(tmpdir) / "input_out"
|
| 274 |
return _parse_fpocket_output(fpocket_out, result.stdout + result.stderr)
|
| 275 |
except subprocess.TimeoutExpired:
|
|
|
|
| 280 |
|
| 281 |
def _parse_fpocket_output(out_dir: Path, raw_output: str) -> FpocketResult:
|
| 282 |
"""Parse fpocket output directory for pocket information."""
|
| 283 |
+
from app.tools.castp import _parse_fpocket
|
| 284 |
+
parsed = _parse_fpocket(out_dir, "custom", 1.4)
|
| 285 |
+
pockets = [{"id": p["id"], "druggability_score": p.get("druggability_score"),
|
| 286 |
+
"volume": p.get("volume_sa"), "area": p.get("area_sa"),
|
| 287 |
+
"score": p.get("score"), "num_residues": p["num_residues"]}
|
| 288 |
+
for p in parsed.get("pockets", [])]
|
| 289 |
+
return FpocketResult(pocket_count=len(pockets), pockets=pockets,
|
| 290 |
+
raw_output=raw_output,
|
| 291 |
+
status="complete" if parsed["status"] == "VALID" else "error")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
| 292 |
|
| 293 |
|
| 294 |
# ── Step 4b: SASA concave-packing fallback (no fpocket / scipy needed) ──────
|
requirements.txt
CHANGED
|
@@ -19,7 +19,7 @@ numpy
|
|
| 19 |
openmm>=8.0 # OpenMM is on PyPI since 8.0 (pip-installable in Dockerfile and Render)
|
| 20 |
openbabel-wheel # obabel CLI — PDBQT prep for docking (SMILES→PDBQT, PDB→PDBQT receptor)
|
| 21 |
pymol-open-source-whl>=3.2 # unofficial wheel of open-source PyMOL; provides the `pymol2` module (headless cleanup + .pse export); 3.2 line supports numpy>=2
|
| 22 |
-
|
| 23 |
|
| 24 |
# dev / testing
|
| 25 |
pytest>=7.0
|
|
|
|
| 19 |
openmm>=8.0 # OpenMM is on PyPI since 8.0 (pip-installable in Dockerfile and Render)
|
| 20 |
openbabel-wheel # obabel CLI — PDBQT prep for docking (SMILES→PDBQT, PDB→PDBQT receptor)
|
| 21 |
pymol-open-source-whl>=3.2 # unofficial wheel of open-source PyMOL; provides the `pymol2` module (headless cleanup + .pse export); 3.2 line supports numpy>=2
|
| 22 |
+
primer3-py==2.3.1
|
| 23 |
|
| 24 |
# dev / testing
|
| 25 |
pytest>=7.0
|
tests/test_castp_provenance.py
CHANGED
|
@@ -55,47 +55,26 @@ async def test_fpocket_path_when_engine_available(monkeypatch):
|
|
| 55 |
|
| 56 |
|
| 57 |
|
| 58 |
-
async def
|
| 59 |
import app.tools.castp as m
|
| 60 |
-
|
| 61 |
-
monkeypatch.setattr(m, "FPOCKET_BIN", "C:/fake/fpocket.exe")
|
| 62 |
monkeypatch.setattr(Path, "exists", lambda self: False)
|
| 63 |
monkeypatch.setattr(shutil, "which", lambda name: None)
|
| 64 |
-
|
| 65 |
-
|
| 66 |
-
|
| 67 |
result = await m._analyze_pockets(_sample_pdb(), "TEST", 1.4)
|
| 68 |
-
assert result["
|
| 69 |
-
assert result["
|
|
|
|
| 70 |
assert result["pockets"] == []
|
| 71 |
-
tried = {t["method"]: t["status"] for t in result["methods_tried"]}
|
| 72 |
-
assert tried["fpocket"] == "unavailable"
|
| 73 |
-
assert tried["sasa_heuristic"] == "ran_no_pockets"
|
| 74 |
-
assert "CASTp" in result["note"] and "fpocket" in result["note"]
|
| 75 |
-
|
| 76 |
|
| 77 |
|
| 78 |
-
async def
|
| 79 |
import app.tools.castp as m
|
| 80 |
-
|
| 81 |
-
|
| 82 |
-
|
| 83 |
-
|
| 84 |
-
|
| 85 |
-
monkeypatch.setattr(Path, "exists", lambda self: False)
|
| 86 |
-
monkeypatch.setattr(shutil, "which", lambda name: None)
|
| 87 |
-
monkeypatch.setattr(
|
| 88 |
-
m.asyncio, "get_running_loop",
|
| 89 |
-
lambda: _BrokenLoop(),
|
| 90 |
-
)
|
| 91 |
-
|
| 92 |
-
result = await m._analyze_pockets(_sample_pdb(), "TEST", 1.4)
|
| 93 |
-
assert result["status"] == "FAILED"
|
| 94 |
-
assert result["pockets"] == []
|
| 95 |
-
assert result["method"] == "none"
|
| 96 |
-
assert result["fallback_used"] is True
|
| 97 |
-
assert "SASA engine crashed" in result["error"]
|
| 98 |
-
assert result["methods_tried"][0] == {"method": "fpocket", "status": "unavailable"}
|
| 99 |
|
| 100 |
|
| 101 |
def test_fpocket_parse_reads_real_output(tmp_path):
|
|
@@ -136,3 +115,60 @@ def test_fpocket_parse_reads_real_output(tmp_path):
|
|
| 136 |
assert result["pockets"][1]["volume_sa"] == 50.0
|
| 137 |
assert result["pockets"][1]["area_sa"] == 80.0
|
| 138 |
assert result["pockets"][1]["centroid"] == [1.0, 1.0, 1.0]
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
| 55 |
|
| 56 |
|
| 57 |
|
| 58 |
+
async def test_unavailable_fpocket_does_not_substitute_heuristic(monkeypatch):
|
| 59 |
import app.tools.castp as m
|
|
|
|
|
|
|
| 60 |
monkeypatch.setattr(Path, "exists", lambda self: False)
|
| 61 |
monkeypatch.setattr(shutil, "which", lambda name: None)
|
| 62 |
+
def forbidden(*args):
|
| 63 |
+
raise AssertionError("Heuristic must be explicitly selected")
|
| 64 |
+
monkeypatch.setattr(m, "_analyze_pockets_sasa_sync", forbidden)
|
| 65 |
result = await m._analyze_pockets(_sample_pdb(), "TEST", 1.4)
|
| 66 |
+
assert result["status"] == "FAILED"
|
| 67 |
+
assert result["method"] == "fpocket"
|
| 68 |
+
assert result["fallback_used"] is False
|
| 69 |
assert result["pockets"] == []
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
| 70 |
|
| 71 |
|
| 72 |
+
async def test_explicit_heuristic_is_labelled(monkeypatch):
|
| 73 |
import app.tools.castp as m
|
| 74 |
+
result = await m._analyze_pockets(_sample_pdb(), "TEST", 1.4, method="sasa_heuristic")
|
| 75 |
+
assert result["evidence_class"] == "Heuristic"
|
| 76 |
+
assert "SASA heuristic" in result["method"]
|
| 77 |
+
assert result["fallback_used"] is False
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
| 78 |
|
| 79 |
|
| 80 |
def test_fpocket_parse_reads_real_output(tmp_path):
|
|
|
|
| 115 |
assert result["pockets"][1]["volume_sa"] == 50.0
|
| 116 |
assert result["pockets"][1]["area_sa"] == 80.0
|
| 117 |
assert result["pockets"][1]["centroid"] == [1.0, 1.0, 1.0]
|
| 118 |
+
|
| 119 |
+
|
| 120 |
+
def test_missing_fpocket_metrics_are_not_zero(tmp_path):
|
| 121 |
+
(tmp_path / "input_info.txt").write_text("Pocket 1\nScore: 0.3\n")
|
| 122 |
+
pocket = _parse_fpocket(tmp_path, "TEST", 1.4)["pockets"][0]
|
| 123 |
+
assert pocket["volume_sa"] is None
|
| 124 |
+
assert pocket["area_sa"] is None
|
| 125 |
+
assert pocket["radius"] is None
|
| 126 |
+
assert pocket["centroid"] is None
|
| 127 |
+
|
| 128 |
+
|
| 129 |
+
def test_structure_prep_reads_native_fpocket_file_and_preserves_missing(tmp_path):
|
| 130 |
+
from app.tools.structure_prep import _parse_fpocket_output
|
| 131 |
+
(tmp_path / "input_info.txt").write_text("Pocket 1\nDruggability Score: 0.2\nVolume: 100.0\n")
|
| 132 |
+
result = _parse_fpocket_output(tmp_path, "native output")
|
| 133 |
+
assert result.status == "complete"
|
| 134 |
+
assert result.pockets[0]["volume"] == 100.0
|
| 135 |
+
assert result.pockets[0]["area"] is None
|
| 136 |
+
assert result.pockets[0]["score"] is None
|
| 137 |
+
|
| 138 |
+
|
| 139 |
+
def test_fpocket_does_not_misuse_solvent_probe_as_alpha_sphere_radius(monkeypatch):
|
| 140 |
+
import app.tools.castp as m
|
| 141 |
+
from types import SimpleNamespace
|
| 142 |
+
calls = []
|
| 143 |
+
def run(args, **kwargs):
|
| 144 |
+
calls.append(args)
|
| 145 |
+
return SimpleNamespace(stdout="fpocket", stderr="", returncode=0)
|
| 146 |
+
monkeypatch.setattr(m.subprocess, "run", run)
|
| 147 |
+
m._run_fpocket_analysis("fpocket", _sample_pdb(), "TEST", 1.4)
|
| 148 |
+
execution = next(args for args in calls if "-f" in args)
|
| 149 |
+
assert "-r" not in execution
|
| 150 |
+
|
| 151 |
+
|
| 152 |
+
def test_native_descriptors_cannot_overwrite_scores_or_geometry(tmp_path):
|
| 153 |
+
(tmp_path / "input_info.txt").write_text(
|
| 154 |
+
"Pocket 1 :\nScore : 0.664\nDruggability Score : 0.855\n"
|
| 155 |
+
"Number of Alpha Spheres : 143\nTotal SASA : 214.924\n"
|
| 156 |
+
"Polar SASA : 111.068\nApolar SASA : 103.856\nVolume : 1455.507\n"
|
| 157 |
+
"Mean alpha sphere radius : 3.977\nHydrophobicity score: 39.727\n"
|
| 158 |
+
"Volume score: 3.939\nPolarity score: 13\nCharge score : -2\n"
|
| 159 |
+
"Alpha sphere density: 8.233\n")
|
| 160 |
+
pocket = _parse_fpocket(tmp_path, "TEST", 1.4)["pockets"][0]
|
| 161 |
+
assert pocket["score"] == 0.664
|
| 162 |
+
assert pocket["druggability_score"] == 0.855
|
| 163 |
+
assert pocket["volume_sa"] == 1455.507
|
| 164 |
+
assert pocket["area_sa"] == 214.924
|
| 165 |
+
assert pocket["alpha_spheres"] == 143
|
| 166 |
+
assert pocket["method_metrics"]["fpocket_total_sasa"] == 214.924
|
| 167 |
+
|
| 168 |
+
|
| 169 |
+
def test_fpocket_version_uses_ansi_stripped_banner(monkeypatch):
|
| 170 |
+
from types import SimpleNamespace
|
| 171 |
+
import app.tools.castp as m
|
| 172 |
+
monkeypatch.setattr(m.subprocess, "run", lambda *a, **k: SimpleNamespace(
|
| 173 |
+
stdout="***** POCKET HUNTING BEGINS *****\n:||: \x1b[1mfpocket 4.0\x1b[0m :||:", stderr=""))
|
| 174 |
+
assert m._fpocket_version("fpocket") == "4.0"
|
tests/test_function_research_grade.py
CHANGED
|
@@ -194,3 +194,25 @@ def test_prediction_from_result_empty_scan_is_insufficient(monkeypatch):
|
|
| 194 |
assert result["status"] == "insufficient_evidence"
|
| 195 |
assert result["method"] == "interpro2go_no_evidence"
|
| 196 |
assert result["provenance"]["retrieval_failures"] == []
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
| 194 |
assert result["status"] == "insufficient_evidence"
|
| 195 |
assert result["method"] == "interpro2go_no_evidence"
|
| 196 |
assert result["provenance"]["retrieval_failures"] == []
|
| 197 |
+
|
| 198 |
+
|
| 199 |
+
def test_mapping_lookup_failure_is_not_biological_absence(monkeypatch):
|
| 200 |
+
monkeypatch.setattr(fp, "_fetch_pdb_sequence", lambda _pdb: "ACDEFGHIK")
|
| 201 |
+
monkeypatch.setattr(fp, "_run_interproscan", lambda _seq: [{"accession": "IPR000001"}])
|
| 202 |
+
|
| 203 |
+
def failed_mapping(_accession):
|
| 204 |
+
fp._record_retrieval_failure("InterPro mapping service unavailable")
|
| 205 |
+
return []
|
| 206 |
+
|
| 207 |
+
monkeypatch.setattr(fp, "_fetch_interpro_go_terms", failed_mapping)
|
| 208 |
+
result = fp.predict_function("1abc")
|
| 209 |
+
assert result["status"] == "evidence_unavailable"
|
| 210 |
+
assert result["provenance"]["retrieval_is_live"] is False
|
| 211 |
+
assert result["go_terms"] == []
|
| 212 |
+
|
| 213 |
+
|
| 214 |
+
def test_signature_only_match_is_not_queried_as_interpro_entry(monkeypatch):
|
| 215 |
+
queried = []
|
| 216 |
+
monkeypatch.setattr(fp, "_fetch_interpro_go_terms", lambda accession: queried.append(accession) or [])
|
| 217 |
+
assert fp._interpro_to_go([{"accession": "PF00001"}, {"accession": "IPR000001"}]) == []
|
| 218 |
+
assert queried == ["IPR000001"]
|
tests/test_primer_advanced.py
CHANGED
|
@@ -33,6 +33,6 @@ def test_multiplex_reports_tm_spread_and_heuristic_boundary():
|
|
| 33 |
{"id": "b", "left_seq": "GATCGATAGATCGATAGATC", "right_seq": "CTAGCTACCTAGCTACCTAG", "left_tm": 61.0, "right_tm": 61.5},
|
| 34 |
])
|
| 35 |
assert result["pair_count"] == 2
|
| 36 |
-
assert result["tm_spread_c"]
|
| 37 |
assert result["evidence_class"] == "Heuristic"
|
| 38 |
assert "empirical multiplex PCR" in result["limitation"]
|
|
|
|
| 33 |
{"id": "b", "left_seq": "GATCGATAGATCGATAGATC", "right_seq": "CTAGCTACCTAGCTACCTAG", "left_tm": 61.0, "right_tm": 61.5},
|
| 34 |
])
|
| 35 |
assert result["pair_count"] == 2
|
| 36 |
+
assert result["tm_spread_c"] > 1.5 # submitted Tm values cannot override recalculation
|
| 37 |
assert result["evidence_class"] == "Heuristic"
|
| 38 |
assert "empirical multiplex PCR" in result["limitation"]
|
tests/test_primers_qc.py
CHANGED
|
@@ -117,29 +117,28 @@ class TestStructures:
|
|
| 117 |
def test_hairpin_detected(self):
|
| 118 |
# CCCC stem pairs with GGGG, loop TTT between -> hairpin.
|
| 119 |
report = hairpin_analysis("AAAACCCCTTTGGGG")
|
| 120 |
-
assert report["
|
| 121 |
assert report["dg"] < 0
|
| 122 |
-
assert report["
|
| 123 |
|
| 124 |
def test_no_hairpin_on_plain_seq(self):
|
| 125 |
# A 20-mer empirically verified to contain no inverted repeat of
|
| 126 |
# >= 4 bases (so it cannot fold a stem-loop).
|
| 127 |
report = hairpin_analysis("GCTAAAGACAATTACATAAC")
|
| 128 |
-
assert report["
|
| 129 |
|
| 130 |
def test_self_dimer_palindrome_is_high_risk(self):
|
| 131 |
# GGGGGCCCCC is a palindrome; two copies dimerize fully.
|
| 132 |
report = dimer_analysis("GGGGGCCCCC", "GGGGGCCCCC")
|
| 133 |
assert report["dg"] <= -10
|
| 134 |
assert report["risk"] == "high"
|
| 135 |
-
assert report["
|
| 136 |
-
assert report["involves_b3"]
|
| 137 |
|
| 138 |
def test_self_dimer_random_oligo_is_low_risk(self):
|
| 139 |
# A random 20-mer empirically verified to contain no self-complementary
|
| 140 |
# run of >= 4 canonical base pairs.
|
| 141 |
report = dimer_analysis("TGGCATTTTTATTACACTCA", "TGGCATTTTTATTACACTCA")
|
| 142 |
-
assert report["
|
| 143 |
|
| 144 |
def test_hetero_dimer_complementary_pair(self):
|
| 145 |
# A primer and its reverse complement anneal perfectly.
|
|
@@ -287,3 +286,38 @@ def test_analyze_endpoint_rejects_protein(primer_client):
|
|
| 287 |
"right_seq": "ACDEFGHIKLMNP",
|
| 288 |
})
|
| 289 |
assert resp.status_code == 400
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
| 117 |
def test_hairpin_detected(self):
|
| 118 |
# CCCC stem pairs with GGGG, loop TTT between -> hairpin.
|
| 119 |
report = hairpin_analysis("AAAACCCCTTTGGGG")
|
| 120 |
+
assert report["structure_found"]
|
| 121 |
assert report["dg"] < 0
|
| 122 |
+
assert report["risk_model"].startswith("Screening thresholds")
|
| 123 |
|
| 124 |
def test_no_hairpin_on_plain_seq(self):
|
| 125 |
# A 20-mer empirically verified to contain no inverted repeat of
|
| 126 |
# >= 4 bases (so it cannot fold a stem-loop).
|
| 127 |
report = hairpin_analysis("GCTAAAGACAATTACATAAC")
|
| 128 |
+
assert report["dg_units"] == "kcal/mol"
|
| 129 |
|
| 130 |
def test_self_dimer_palindrome_is_high_risk(self):
|
| 131 |
# GGGGGCCCCC is a palindrome; two copies dimerize fully.
|
| 132 |
report = dimer_analysis("GGGGGCCCCC", "GGGGGCCCCC")
|
| 133 |
assert report["dg"] <= -10
|
| 134 |
assert report["risk"] == "high"
|
| 135 |
+
assert report["structure_found"]
|
|
|
|
| 136 |
|
| 137 |
def test_self_dimer_random_oligo_is_low_risk(self):
|
| 138 |
# A random 20-mer empirically verified to contain no self-complementary
|
| 139 |
# run of >= 4 canonical base pairs.
|
| 140 |
report = dimer_analysis("TGGCATTTTTATTACACTCA", "TGGCATTTTTATTACACTCA")
|
| 141 |
+
assert report["engine"] == "primer3-py"
|
| 142 |
|
| 143 |
def test_hetero_dimer_complementary_pair(self):
|
| 144 |
# A primer and its reverse complement anneal perfectly.
|
|
|
|
| 286 |
"right_seq": "ACDEFGHIKLMNP",
|
| 287 |
})
|
| 288 |
assert resp.status_code == 400
|
| 289 |
+
|
| 290 |
+
|
| 291 |
+
def test_qc_matches_primer3_and_records_conditions():
|
| 292 |
+
import primer3
|
| 293 |
+
from app.tools.oligo_qc import oligo_report, DEFAULT_CONDITIONS
|
| 294 |
+
seq = "ATGGAGGAGCCGCAGTCAGA"
|
| 295 |
+
report = oligo_report(seq)
|
| 296 |
+
c = DEFAULT_CONDITIONS
|
| 297 |
+
assert report["conditions"] == c
|
| 298 |
+
assert report["tm_c"] == pytest.approx(primer3.calc_tm(seq, **{k:v for k,v in c.items() if k != "temp_c"}), abs=0.01)
|
| 299 |
+
assert report["hairpin"]["dg"] == pytest.approx(primer3.calc_hairpin(seq, **c).dg / 1000, abs=0.001)
|
| 300 |
+
assert report["self_dimer"]["dg"] == pytest.approx(primer3.calc_homodimer(seq, **c).dg / 1000, abs=0.001)
|
| 301 |
+
|
| 302 |
+
|
| 303 |
+
def test_qc_rejects_ambiguous_and_unavailable(primer_client, monkeypatch):
|
| 304 |
+
import app.tools.oligo_qc as qc
|
| 305 |
+
resp = primer_client.post("/api/primers/analyze", json={"left_seq":"ACGTNACGT", "right_seq":"ACGTACGT"})
|
| 306 |
+
assert resp.status_code == 400
|
| 307 |
+
def unavailable():
|
| 308 |
+
raise RuntimeError("Primer3 unavailable; no fallback")
|
| 309 |
+
monkeypatch.setattr(qc, "_engine", unavailable)
|
| 310 |
+
resp = primer_client.post("/api/primers/analyze", json={"left_seq":"ACGTACGT", "right_seq":"ACGTACGT"})
|
| 311 |
+
assert resp.status_code == 503
|
| 312 |
+
|
| 313 |
+
|
| 314 |
+
def test_heterodimer_parity_and_custom_conditions():
|
| 315 |
+
import primer3
|
| 316 |
+
from app.tools.oligo_qc import DEFAULT_CONDITIONS, oligo_report
|
| 317 |
+
c = {**DEFAULT_CONDITIONS, "mv_conc": 75.0, "temp_c": 25.0}
|
| 318 |
+
a, b = "ATGGAGGAGCCGCAGTCAGA", "TCAGTCTGAGTCAGGCCCTT"
|
| 319 |
+
assert dimer_analysis(a, b, conditions=c)["dg"] == pytest.approx(primer3.calc_heterodimer(a, b, **c).dg / 1000, abs=0.001)
|
| 320 |
+
report = oligo_report(a, c)
|
| 321 |
+
assert report["conditions"] == c
|
| 322 |
+
assert report["tm_50mM"] is None
|
| 323 |
+
assert report["tm_c"] != oligo_report(a)["tm_c"]
|
tests/test_restriction_map.py
ADDED
|
@@ -0,0 +1,73 @@
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
| 1 |
+
import pytest
|
| 2 |
+
from app.tools.restriction_map import restriction_map
|
| 3 |
+
from app.tools.sequence_utilities import SequenceUtilitiesError
|
| 4 |
+
|
| 5 |
+
|
| 6 |
+
def test_ecori_linear_cuts_and_conservation():
|
| 7 |
+
r = restriction_map("AAAAGAATTCAAAA", ["EcoRI"])
|
| 8 |
+
assert r["sites"][0]["start"] == 5
|
| 9 |
+
assert r["cut_boundaries"] == [5]
|
| 10 |
+
assert r["fragment_lengths"] == [5, 9]
|
| 11 |
+
assert r["fragment_total"] == r["length"]
|
| 12 |
+
|
| 13 |
+
|
| 14 |
+
def test_origin_spanning_site_and_one_cut_linearizes():
|
| 15 |
+
r = restriction_map("AATTCAAAAG", ["EcoRI"], True)
|
| 16 |
+
assert r["sites"][0]["start"] == 10
|
| 17 |
+
assert r["sites"][0]["end"] == 5
|
| 18 |
+
assert r["sites"][0]["spans_origin"]
|
| 19 |
+
assert r["cut_boundaries"] == [0]
|
| 20 |
+
assert r["fragment_lengths"] == [10]
|
| 21 |
+
assert r["product_topology"] == "linear"
|
| 22 |
+
|
| 23 |
+
|
| 24 |
+
def test_both_nonpalindromic_strands_and_linear_end_exclusion():
|
| 25 |
+
r = restriction_map("AAAAAAGGTCTCAAAAAAGAGACCAAAAAA", ["BsaI"])
|
| 26 |
+
assert {hit["strand"] for hit in r["sites"]} == {"+", "-"}
|
| 27 |
+
assert sum(r["fragment_lengths"]) == r["length"]
|
| 28 |
+
r = restriction_map("GGTCTC", ["BsaI"])
|
| 29 |
+
assert r["sites"][0]["complete_cut"] is False
|
| 30 |
+
assert r["cut_boundaries"] == []
|
| 31 |
+
|
| 32 |
+
|
| 33 |
+
def test_overlapping_recognition_sites():
|
| 34 |
+
r = restriction_map("CGTCTCGTCTC", ["BsmBI"], True)
|
| 35 |
+
assert [h["start"] for h in r["sites"]] == [1, 6]
|
| 36 |
+
assert sum(r["fragment_lengths"]) == 11
|
| 37 |
+
|
| 38 |
+
|
| 39 |
+
def test_circular_uncut_remains_circular_and_multi_enzyme_conservation():
|
| 40 |
+
assert restriction_map("AAAAAA", ["EcoRI"], True)["product_topology"] == "circular"
|
| 41 |
+
for circular in (False, True):
|
| 42 |
+
r = restriction_map("GAATTCAAAAGGATCCAAAAGAATTC", ["EcoRI", "BamHI", "EcoRI"], circular)
|
| 43 |
+
assert r["fragment_total"] == r["length"]
|
| 44 |
+
assert len(r["enzymes"]) == 2
|
| 45 |
+
|
| 46 |
+
|
| 47 |
+
def test_ambiguity_and_unknown_enzyme_rejected():
|
| 48 |
+
with pytest.raises(SequenceUtilitiesError):
|
| 49 |
+
restriction_map("GAANTC", ["EcoRI"])
|
| 50 |
+
with pytest.raises(SequenceUtilitiesError):
|
| 51 |
+
restriction_map("GAATTC", ["Unknown"])
|
| 52 |
+
|
| 53 |
+
|
| 54 |
+
def test_forward_and_reverse_type_iis_cuts_match_biopython():
|
| 55 |
+
from Bio import Restriction
|
| 56 |
+
from Bio.Seq import Seq
|
| 57 |
+
for enzyme, sequence in [(Restriction.BsaI, "AAAAAAGGTCTCAAAAAAGAGACCAAAAAA"), (Restriction.BsmBI, "AAAAAACGTCTCAAAAAAGAGACGAAAAAA")]:
|
| 58 |
+
for circular in (False, True):
|
| 59 |
+
mapped = restriction_map(sequence, [str(enzyme)], circular)
|
| 60 |
+
reference = sorted({(cut - 1) % len(sequence) if circular else cut - 1 for cut in enzyme.search(Seq(sequence), linear=not circular)})
|
| 61 |
+
assert mapped["cut_boundaries"] == reference
|
| 62 |
+
|
| 63 |
+
|
| 64 |
+
@pytest.mark.asyncio
|
| 65 |
+
async def test_endpoint_preserves_scientific_evidence_and_rejects_ambiguity():
|
| 66 |
+
from fastapi import HTTPException
|
| 67 |
+
from app.routers.seq_tools import RestrictionMapRequest, restriction_map_endpoint
|
| 68 |
+
result = await restriction_map_endpoint(RestrictionMapRequest(sequence="GAATTCAAAA", enzymes=["EcoRI"], circular=True))
|
| 69 |
+
assert result["results"]["fragment_total"] == 10
|
| 70 |
+
assert result["validation"]["experimental_digestion_validated"] is False
|
| 71 |
+
with pytest.raises(HTTPException) as exc:
|
| 72 |
+
await restriction_map_endpoint(RestrictionMapRequest(sequence="NGAATTC"))
|
| 73 |
+
assert exc.value.status_code == 400
|
tests/test_seq_tools.py
CHANGED
|
@@ -73,13 +73,14 @@ class TestAnalyzeSequence:
|
|
| 73 |
tr = res["translation"]
|
| 74 |
assert tr is not None
|
| 75 |
assert tr["frames"]["1"] == "MH*RA*"
|
| 76 |
-
assert set(tr["frames"].keys()) == {"1", "2", "3", "
|
| 77 |
-
# Best ORF is now genome-wide (six frames): the reverse strand's frame
|
| 78 |
# carries MHVNA (5 aa), longer than the forward MH (which stops early).
|
| 79 |
-
assert tr["best"]["frame"] ==
|
| 80 |
-
assert tr["best"]["strand"] == "
|
| 81 |
assert tr["best"]["protein"] == "MHVNA"
|
| 82 |
-
assert tr["best"]["start"] ==
|
|
|
|
| 83 |
|
| 84 |
def test_protein_mw_and_composition(self):
|
| 85 |
res = analyze_sequence("MEEPQSDPSVEP", seq_type="protein")
|
|
|
|
| 73 |
tr = res["translation"]
|
| 74 |
assert tr is not None
|
| 75 |
assert tr["frames"]["1"] == "MH*RA*"
|
| 76 |
+
assert set(tr["frames"].keys()) == {"1", "2", "3", "-1", "-2", "-3"}
|
| 77 |
+
# Best ORF is now genome-wide (six frames): the reverse strand's frame -1
|
| 78 |
# carries MHVNA (5 aa), longer than the forward MH (which stops early).
|
| 79 |
+
assert tr["best"]["frame"] == -1
|
| 80 |
+
assert tr["best"]["strand"] == "-"
|
| 81 |
assert tr["best"]["protein"] == "MHVNA"
|
| 82 |
+
assert tr["best"]["start"] == 2
|
| 83 |
+
assert tr["best"]["end"] == 16
|
| 84 |
|
| 85 |
def test_protein_mw_and_composition(self):
|
| 86 |
res = analyze_sequence("MEEPQSDPSVEP", seq_type="protein")
|
tests/test_workflow_manifest.py
ADDED
|
@@ -0,0 +1,52 @@
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
| 1 |
+
import asyncio
|
| 2 |
+
import hashlib
|
| 3 |
+
import json
|
| 4 |
+
|
| 5 |
+
import pytest
|
| 6 |
+
from fastapi import HTTPException
|
| 7 |
+
|
| 8 |
+
from app.services.workflow_manifest import build_workflow_manifest
|
| 9 |
+
|
| 10 |
+
|
| 11 |
+
def test_states_versions_dependencies_and_canonical_integrity():
|
| 12 |
+
exp = {"experiment_id": "e", "job_id": "j", "input_hash": "abc", "parameters": {"seed": 7}}
|
| 13 |
+
provenance = [{"node_id": "blast", "status": "complete", "tool": "BLAST", "tool_version": "2.16.0",
|
| 14 |
+
"database": "nr", "deps": ["input"], "params": {"evalue": 1e-5}},
|
| 15 |
+
{"node_id": "msa", "status": "skipped"},
|
| 16 |
+
{"node_id": "phylo", "status": "failed"},
|
| 17 |
+
{"node_id": "ai", "status": "running"}, {"node_id": "unknown"}]
|
| 18 |
+
manifest = build_workflow_manifest(exp, provenance)
|
| 19 |
+
assert [n["state"] for n in manifest["steps"]] == ["COMPLETED", "SKIPPED", "FAILED", "RUNNING", "UNKNOWN"]
|
| 20 |
+
assert manifest["steps"][0]["tool_version"] == "2.16.0"
|
| 21 |
+
assert manifest["steps"][0]["database_version"] is None
|
| 22 |
+
assert manifest["edges"] == [{"from": "input", "to": "blast"}]
|
| 23 |
+
assert any(g["reason"] == "UNRECORDED_DEPENDENCY" for g in manifest["reproducibility_gaps"])
|
| 24 |
+
assert manifest["execution"]["executable"] is False
|
| 25 |
+
assert manifest == build_workflow_manifest(exp, provenance)
|
| 26 |
+
digest = manifest.pop("manifest_sha256")
|
| 27 |
+
assert digest == hashlib.sha256(json.dumps(manifest, sort_keys=True, separators=(",", ":"), ensure_ascii=False).encode()).hexdigest()
|
| 28 |
+
|
| 29 |
+
|
| 30 |
+
def test_empty_evidence_cannot_imply_execution():
|
| 31 |
+
result = build_workflow_manifest({"status": "complete"}, [])
|
| 32 |
+
assert result["steps"] == []
|
| 33 |
+
assert any(g["field"] == "steps" for g in result["reproducibility_gaps"])
|
| 34 |
+
assert result["input_sha256"] is None
|
| 35 |
+
|
| 36 |
+
|
| 37 |
+
def test_export_rejects_other_tenant_before_reading_experiment(monkeypatch):
|
| 38 |
+
from app.routers import experiments
|
| 39 |
+
monkeypatch.setattr(experiments, "owns_job", lambda job, user: False)
|
| 40 |
+
monkeypatch.setattr(experiments, "get_experiment", lambda _: pytest.fail("Unauthorized experiment read"))
|
| 41 |
+
with pytest.raises(HTTPException) as error:
|
| 42 |
+
asyncio.run(experiments.export_workflow_manifest("foreign-job", "user"))
|
| 43 |
+
assert error.value.status_code == 404
|
| 44 |
+
|
| 45 |
+
|
| 46 |
+
def test_export_owned_job_preserves_explicit_empty_provenance(monkeypatch):
|
| 47 |
+
from app.routers import experiments
|
| 48 |
+
monkeypatch.setattr(experiments, "owns_job", lambda job, user: True)
|
| 49 |
+
monkeypatch.setattr(experiments, "get_experiment", lambda _: {"experiment_id": "e", "provenance": []})
|
| 50 |
+
monkeypatch.setattr(experiments, "provenance_for_experiment", lambda _: pytest.fail("Unexpected fallback"))
|
| 51 |
+
result = asyncio.run(experiments.export_workflow_manifest("j", "user"))
|
| 52 |
+
assert result["steps"] == []
|