From 2e912add8100726469a1ffef8bb11a6a3820d3da Mon Sep 17 00:00:00 2001 From: Divide-By-0 Date: Thu, 27 Aug 2026 23:19:20 -0700 Subject: [PATCH 1/3] Add cloning sequence inventory and source indexes --- src/lab_bench_2/README.md | 4 + src/lab_bench_2/cloning_inventory.py | 482 ++++++++++++++++++++ tests/lab_bench_2/test_cloning_inventory.py | 109 +++++ tools/inventory_cloning_sequences.py | 66 +++ 4 files changed, 661 insertions(+) create mode 100644 src/lab_bench_2/cloning_inventory.py create mode 100644 tests/lab_bench_2/test_cloning_inventory.py create mode 100644 tools/inventory_cloning_sequences.py diff --git a/src/lab_bench_2/README.md b/src/lab_bench_2/README.md index d776d02..ee5f797 100644 --- a/src/lab_bench_2/README.md +++ b/src/lab_bench_2/README.md @@ -2,6 +2,10 @@ LABBench2 evaluates language models and research agents on life-science tasks spanning literature reasoning, database access, figures, tables, protocols, source quality, sequence analysis, cloning, patents, and clinical trials. +For cloning attachment acquisition, feature/primer/enzyme inventories, the +current data audit, and the proposed question-difficulty graph, see +[Cloning sequence data](CLONING_SEQUENCE_DATA.md). + Contributed by [@iphan](https://github.com/iphan), [@ItsTania](https://github.com/ItsTania), [@lewtun](https://github.com/lewtun), [@Sunishchal](https://github.com/Sunishchal) diff --git a/src/lab_bench_2/cloning_inventory.py b/src/lab_bench_2/cloning_inventory.py new file mode 100644 index 0000000..165e9d3 --- /dev/null +++ b/src/lab_bench_2/cloning_inventory.py @@ -0,0 +1,482 @@ +"""Build provenance-preserving inventories for cloning sequence attachments.""" + +from __future__ import annotations + +import hashlib +import re +import warnings +from collections import defaultdict +from pathlib import Path +from typing import Any, Iterable, Mapping, Sequence + +SEQUENCE_FORMATS = { + ".fa": "fasta", + ".fasta": "fasta", + ".fna": "fasta", + ".gb": "genbank", + ".gbff": "genbank", + ".gbk": "genbank", +} +LABEL_QUALIFIERS = ("label", "gene", "locus_tag", "product", "standard_name") +TEXT_QUALIFIERS = (*LABEL_QUALIFIERS, "note", "function", "regulatory_class") +PRIMER_TYPES = {"primer", "primer_bind"} +SOURCE_TYPES = {"source"} +SHORT_TERM_MAX_LENGTH = 4 + +ROLE_DESCRIPTIONS = { + "affinity_tag": "affinity or purification tag", + "coding_sequence": "protein-coding sequence", + "enhancer": "transcriptional enhancer", + "epitope_tag": "epitope or detection tag", + "genome_editor": "genome-editing effector", + "intron_or_splicing": "intron or RNA-splicing element", + "localization_signal": "subcellular localization signal", + "origin_of_replication": "replication origin", + "primer_binding_site": "primer binding site", + "promoter": "transcriptional promoter", + "protein_binding_site": "protein binding site", + "recombination_site": "recombination site", + "reporter": "reporter or detectable output", + "selection_marker": "selectable marker", + "terminator_or_polya": "transcription terminator or polyadenylation signal", + "translation_control": "translation-control element", + "viral_packaging_or_expression": "viral packaging or expression element", +} + +ROLE_PATTERNS: tuple[tuple[str, tuple[str, ...]], ...] = ( + ( + "localization_signal", + ("nls", "nuclear localization", "nes", "localization signal", "signal peptide"), + ), + ("epitope_tag", ("flag", "ha tag", "myc tag", "v5 tag", "epitope tag")), + ( + "affinity_tag", + ( + "his tag", + "histidine tag", + "his6", + "6xhis", + "strep tag", + "maltose binding protein", + "mbp tag", + "gst tag", + ), + ), + ( + "selection_marker", + ( + "ampr", + "kanr", + "camr", + "chloramphenicol resistance", + "ampicillin resistance", + "kanamycin resistance", + "puromycin resistance", + "hygromycin resistance", + "blasticidin resistance", + "selectable marker", + ), + ), + ( + "reporter", + ("gfp", "mcherry", "tdtomato", "luciferase", "fluorescent protein", "reporter"), + ), + ("genome_editor", ("cas9", "cas12", "cpf1", "base editor", "prime editor")), + ( + "origin_of_replication", + ( + "origin of replication", + "replication origin", + "cole1", + "pmb1 ori", + "f1 ori", + "sv40 ori", + ), + ), + ( + "terminator_or_polya", + ("terminator", "polyadenylation", "poly(a)", "polya signal", "poly a signal"), + ), + ("promoter", ("promoter",)), + ("enhancer", ("enhancer",)), + ("translation_control", ("ires", "ribosome binding site", "rbs", "kozak")), + ( + "viral_packaging_or_expression", + ( + "woodchuck", + "wpre", + "ltr", + "packaging signal", + "rev response", + "posttranscriptional regulatory", + ), + ), + ("recombination_site", ("loxp", "frt", "attb", "attp", "recombination site")), + ("intron_or_splicing", ("intron", "splice donor", "splice acceptor")), +) + + +def discover_sequence_files(inputs: Sequence[Path]) -> list[Path]: + """Return supported sequence files under explicit files or directories.""" + paths: set[Path] = set() + for input_path in inputs: + path = input_path.expanduser().resolve() + if path.is_file() and path.suffix.lower() in SEQUENCE_FORMATS: + paths.add(path) + elif path.is_dir(): + paths.update( + candidate.resolve() + for candidate in path.rglob("*") + if candidate.is_file() and candidate.suffix.lower() in SEQUENCE_FORMATS + ) + return sorted(paths) + + +def build_cloning_inventory( + inputs: Sequence[Path], + *, + root: Path | None = None, + include_enzymes: bool = True, +) -> dict[str, Any]: + """Parse sequence attachments and build file, feature, primer, and enzyme indexes.""" + files = discover_sequence_files(inputs) + base = (root or _common_parent(files, inputs)).expanduser().resolve() + feature_index: defaultdict[str, set[str]] = defaultdict(set) + role_index: defaultdict[str, set[str]] = defaultdict(set) + primer_index: defaultdict[str, set[str]] = defaultdict(set) + sequence_index: defaultdict[str, set[str]] = defaultdict(set) + enzyme_index: defaultdict[str, dict[str, list[str]]] = defaultdict( + lambda: {"zero_cutters": [], "single_cutters": [], "multi_cutters": []} + ) + entries: list[dict[str, Any]] = [] + errors: list[dict[str, str]] = [] + enzyme_catalog = _neb_enzyme_catalog() if include_enzymes else [] + + for path in files: + display_path = _display_path(path, base) + try: + entry = _inventory_file(path, display_path, include_enzymes) + except Exception as exc: # noqa: BLE001 - preserve per-file parse failures + errors.append( + {"path": display_path, "error": f"{type(exc).__name__}: {exc}"} + ) + continue + entries.append(entry) + for record in entry["records"]: + sequence_index[record["sequence_sha256"]].add(display_path) + for feature in record["features"]: + if feature["normalized_label"]: + feature_index[feature["normalized_label"]].add(display_path) + for role in feature["functional_roles"]: + role_index[role].add(display_path) + for primer in record["primers"]: + if primer["normalized_label"]: + primer_index[primer["normalized_label"]].add(display_path) + for name, cut_count in record.get("neb_restriction_sites", {}).items(): + category = ( + "zero_cutters" + if cut_count == 0 + else "single_cutters" + if cut_count == 1 + else "multi_cutters" + ) + enzyme_index[name][category].append(display_path) + + no_features = [ + entry["path"] for entry in entries if entry["part_feature_count"] == 0 + ] + no_primers = [entry["path"] for entry in entries if entry["primer_count"] == 0] + duplicate_sequences = { + digest: sorted(paths) + for digest, paths in sequence_index.items() + if len(paths) > 1 + } + files_with_parse_warnings = [ + entry["path"] for entry in entries if entry["parse_warnings"] + ] + return { + "schema_version": 1, + "root": str(base), + "summary": { + "discovered_file_count": len(files), + "parsed_file_count": len(entries), + "parse_error_count": len(errors), + "parse_warning_count": sum( + len(entry["parse_warnings"]) for entry in entries + ), + "files_with_parse_warnings": files_with_parse_warnings, + "files_with_no_part_features": no_features, + "files_with_no_primers": no_primers, + "duplicate_sequence_group_count": len(duplicate_sequences), + }, + "files": entries, + "indexes": { + "feature_label_to_files": _sorted_set_mapping(feature_index), + "functional_role_to_files": _sorted_set_mapping(role_index), + "primer_label_to_files": _sorted_set_mapping(primer_index), + "sequence_sha256_to_files": _sorted_set_mapping(sequence_index), + "duplicate_sequences": duplicate_sequences, + "neb_enzyme_to_files": dict(sorted(enzyme_index.items())), + }, + "neb_enzyme_catalog": enzyme_catalog, + "errors": errors, + "provenance": { + "feature_annotations": "source GenBank qualifiers plus deterministic role rules", + "primer_sequences": "5-prime-to-3-prime feature.extract(record.seq)", + "neb_enzyme_catalog": ( + "Biopython Restriction_Dictionary supplier code N; verify live NEB availability" + if include_enzymes + else "not generated" + ), + }, + } + + +def _inventory_file( + path: Path, display_path: str, include_enzymes: bool +) -> dict[str, Any]: + from Bio import BiopythonParserWarning, SeqIO + + file_format = SEQUENCE_FORMATS[path.suffix.lower()] + raw = path.read_bytes() + with warnings.catch_warnings(record=True) as caught_warnings: + warnings.simplefilter("always", BiopythonParserWarning) + with path.open(encoding="utf-8") as handle: + parsed_records = SeqIO.parse( # type: ignore[no-untyped-call] + handle, file_format + ) + records = list(parsed_records) + if not records: + raise ValueError("no sequence records") + record_entries = [ + _inventory_record(record, include_enzymes=include_enzymes) for record in records + ] + return { + "path": display_path, + "format": file_format, + "file_sha256": hashlib.sha256(raw).hexdigest(), + "record_count": len(record_entries), + "sequence_length": sum(record["sequence_length"] for record in record_entries), + "feature_count": sum(record["feature_count"] for record in record_entries), + "annotation_feature_count": sum( + record["annotation_feature_count"] for record in record_entries + ), + "part_feature_count": sum( + record["part_feature_count"] for record in record_entries + ), + "primer_count": sum(record["primer_count"] for record in record_entries), + "parse_warnings": [ + { + "category": warning.category.__name__, + "message": str(warning.message), + } + for warning in caught_warnings + ], + "records": record_entries, + } + + +def _inventory_record(record: Any, *, include_enzymes: bool) -> dict[str, Any]: + sequence = str(record.seq).upper() + feature_entries = [ + _inventory_feature(feature) + for feature in record.features + if feature.type.casefold() not in SOURCE_TYPES | PRIMER_TYPES + ] + primers = [ + _inventory_primer(record, feature) + for feature in record.features + if feature.type.lower() in PRIMER_TYPES + ] + result: dict[str, Any] = { + "id": str(record.id), + "name": str(record.name), + "description": str(record.description), + "sequence_length": len(sequence), + "sequence_sha256": hashlib.sha256(sequence.encode("ascii")).hexdigest(), + "topology": str(record.annotations.get("topology") or "unknown"), + "feature_count": len(record.features), + "annotation_feature_count": sum( + feature.type.casefold() not in SOURCE_TYPES for feature in record.features + ), + "part_feature_count": len(feature_entries), + "primer_count": len(primers), + "features": feature_entries, + "primers": primers, + } + if include_enzymes: + result["neb_restriction_sites"] = _neb_restriction_sites(record) + return result + + +def _inventory_feature(feature: Any) -> dict[str, Any]: + qualifiers = _clean_qualifiers(feature.qualifiers) + label = _first_qualifier(qualifiers, LABEL_QUALIFIERS) + roles = classify_feature_roles(feature.type, qualifiers) + return { + "feature_type": str(feature.type), + "location": str(feature.location), + "start_0_based": int(feature.location.start), + "end_0_based_exclusive": int(feature.location.end), + "strand": feature.location.strand, + "label": label, + "normalized_label": normalize_label(label), + "functional_roles": roles, + "functional_description": describe_feature(label, roles), + "qualifiers": qualifiers, + } + + +def _inventory_primer(record: Any, feature: Any) -> dict[str, Any]: + qualifiers = _clean_qualifiers(feature.qualifiers) + label = _first_qualifier(qualifiers, LABEL_QUALIFIERS) or "unnamed primer" + return { + "label": label, + "normalized_label": normalize_label(label), + "location": str(feature.location), + "start_0_based": int(feature.location.start), + "end_0_based_exclusive": int(feature.location.end), + "strand": feature.location.strand, + "binding_sequence_5to3": str(feature.extract(record.seq)).upper(), + "length": len(feature), + "qualifiers": qualifiers, + } + + +def classify_feature_roles( + feature_type: str, qualifiers: Mapping[str, list[str]] +) -> list[str]: + """Map a structural GenBank feature to zero or more functional roles.""" + feature_type_lower = feature_type.casefold() + text = " ".join( + value for key in TEXT_QUALIFIERS for value in qualifiers.get(key, []) + ).casefold() + roles: set[str] = set() + if feature_type_lower == "cds": + roles.add("coding_sequence") + if feature_type_lower in PRIMER_TYPES: + roles.add("primer_binding_site") + structural_roles = { + "promoter": "promoter", + "enhancer": "enhancer", + "terminator": "terminator_or_polya", + "polya_signal": "terminator_or_polya", + "rep_origin": "origin_of_replication", + "protein_bind": "protein_binding_site", + "intron": "intron_or_splicing", + } + role = structural_roles.get(feature_type_lower) + if role: + roles.add(role) + for candidate, patterns in ROLE_PATTERNS: + if any(_contains_term(text, pattern) for pattern in patterns): + roles.add(candidate) + return sorted(roles) + + +def describe_feature(label: str, roles: Sequence[str]) -> str: + """Create a conservative prompt-ready function description.""" + descriptions = [ + ROLE_DESCRIPTIONS[role] for role in roles if role in ROLE_DESCRIPTIONS + ] + if label and descriptions: + return f"{label}: {', '.join(descriptions)}" + if label: + return label + return ", ".join(descriptions) or "unclassified sequence annotation" + + +def normalize_label(label: str) -> str: + """Normalize a feature/primer label for cross-file indexing.""" + return re.sub(r"[^a-z0-9]+", " ", label.casefold()).strip() + + +def _contains_term(text: str, term: str) -> bool: + if len(term) <= SHORT_TERM_MAX_LENGTH and term.replace(" ", "").isalnum(): + return bool(re.search(rf"(? dict[str, list[str]]: + cleaned: dict[str, list[str]] = {} + for key, raw_values in qualifiers.items(): + values = raw_values if isinstance(raw_values, list) else [raw_values] + cleaned[str(key)] = [str(value) for value in values] + return cleaned + + +def _first_qualifier(qualifiers: Mapping[str, list[str]], keys: Iterable[str]) -> str: + for key in keys: + values = qualifiers.get(key, []) + if values and values[0].strip(): + return values[0].strip() + return "" + + +def _neb_enzyme_names() -> list[str]: + from Bio.Restriction import Restriction_Dictionary + + supplier = Restriction_Dictionary.suppliers["N"] + return sorted(str(name) for name in supplier[1]) + + +def _neb_enzyme_catalog() -> list[dict[str, Any]]: + from Bio import __version__ as biopython_version + from Bio.Restriction import Restriction_Dictionary + + catalog = [] + for name in _neb_enzyme_names(): + data = Restriction_Dictionary.rest_dict[name] + catalog.append( + { + "name": name, + "recognition_site": data.get("site"), + "cut_5prime": data.get("fst5"), + "cut_3prime": data.get("fst3"), + "second_cut_5prime": data.get("scd5"), + "second_cut_3prime": data.get("scd3"), + "overhang_length": data.get("ovhg"), + "overhang_sequence": data.get("ovhgseq"), + "optimal_temperature_c": data.get("opt_temp"), + "inactivation_temperature_c": data.get("inact_temp"), + "rebase_uri": data.get("uri"), + "supplier_code": "N", + "biopython_version": biopython_version, + } + ) + return catalog + + +def _neb_restriction_sites(record: Any) -> dict[str, int]: + from Bio.Restriction.Restriction import RestrictionBatch + + batch = RestrictionBatch(_neb_enzyme_names()) # type: ignore[no-untyped-call] + linear = str(record.annotations.get("topology") or "").casefold() != "circular" + matches = batch.search( # type: ignore[no-untyped-call] + record.seq, linear=linear + ) + return dict( + sorted((str(enzyme), len(positions)) for enzyme, positions in matches.items()) + ) + + +def _display_path(path: Path, root: Path) -> str: + try: + return path.relative_to(root).as_posix() + except ValueError: + return str(path) + + +def _common_parent(files: Sequence[Path], inputs: Sequence[Path]) -> Path: + if files: + import os + + return Path(os.path.commonpath([str(path.parent) for path in files])) + if inputs: + first = inputs[0].expanduser().resolve() + return first if first.is_dir() else first.parent + return Path.cwd() + + +def _sorted_set_mapping(mapping: Mapping[str, set[str]]) -> dict[str, list[str]]: + return {key: sorted(values) for key, values in sorted(mapping.items())} diff --git a/tests/lab_bench_2/test_cloning_inventory.py b/tests/lab_bench_2/test_cloning_inventory.py new file mode 100644 index 0000000..6e4c77a --- /dev/null +++ b/tests/lab_bench_2/test_cloning_inventory.py @@ -0,0 +1,109 @@ +from pathlib import Path + +from Bio import SeqIO +from Bio.Seq import Seq +from Bio.SeqFeature import FeatureLocation, SeqFeature +from Bio.SeqRecord import SeqRecord + +from lab_bench_2.cloning_inventory import ( + build_cloning_inventory, + classify_feature_roles, +) + + +def _write_genbank(path: Path) -> None: + record = SeqRecord( + Seq("AAAAGGTCTCAAAAAAAAAAAAAAAAAAAAAA"), + id="plasmid", + name="plasmid", + description="synthetic fixture", + ) + record.annotations = {"molecule_type": "DNA", "topology": "linear"} + record.features = [ + SeqFeature(FeatureLocation(0, 34), type="source"), + SeqFeature( + FeatureLocation(10, 19), + type="CDS", + qualifiers={"label": ["SV40 NLS"]}, + ), + SeqFeature( + FeatureLocation(19, 28), + type="CDS", + qualifiers={"label": ["AmpR"]}, + ), + SeqFeature( + FeatureLocation(0, 4, strand=-1), + type="primer_bind", + qualifiers={"label": ["Reverse primer"]}, + ), + SeqFeature( + FeatureLocation(4, 10), + type="promoter", + qualifiers={"label": ["CMV promoter"]}, + ), + ] + SeqIO.write(record, path, "genbank") + + +def test_builds_feature_primer_duplicate_and_missing_indexes(tmp_path: Path) -> None: + gbk = tmp_path / "plasmid.gbk" + duplicate = tmp_path / "duplicate.gb" + fasta = tmp_path / "bare.fa" + _write_genbank(gbk) + duplicate.write_bytes(gbk.read_bytes()) + SeqIO.write(SeqRecord(Seq("ACGT"), id="bare", description=""), fasta, "fasta") + + inventory = build_cloning_inventory( + [tmp_path], root=tmp_path, include_enzymes=False + ) + + summary = inventory["summary"] + assert summary["discovered_file_count"] == 3 + assert summary["files_with_no_part_features"] == ["bare.fa"] + assert summary["files_with_no_primers"] == ["bare.fa"] + assert inventory["indexes"]["feature_label_to_files"]["sv40 nls"] == [ + "duplicate.gb", + "plasmid.gbk", + ] + assert inventory["indexes"]["functional_role_to_files"]["localization_signal"] == [ + "duplicate.gb", + "plasmid.gbk", + ] + assert inventory["indexes"]["primer_label_to_files"]["reverse primer"] == [ + "duplicate.gb", + "plasmid.gbk", + ] + assert summary["duplicate_sequence_group_count"] == 1 + + plasmid = next( + entry for entry in inventory["files"] if entry["path"] == "plasmid.gbk" + ) + assert plasmid["records"][0]["primers"][0]["binding_sequence_5to3"] == "TTTT" + nls = next( + feature + for feature in plasmid["records"][0]["features"] + if feature["label"] == "SV40 NLS" + ) + assert nls["functional_roles"] == ["coding_sequence", "localization_signal"] + assert "subcellular localization signal" in nls["functional_description"] + + +def test_inventory_includes_neb_catalog_and_cut_counts(tmp_path: Path) -> None: + gbk = tmp_path / "plasmid.gbk" + _write_genbank(gbk) + + inventory = build_cloning_inventory([gbk], root=tmp_path) + + bsa_i = next( + enzyme for enzyme in inventory["neb_enzyme_catalog"] if enzyme["name"] == "BsaI" + ) + assert bsa_i["recognition_site"] == "GGTCTC" + assert inventory["files"][0]["records"][0]["neb_restriction_sites"]["BsaI"] == 1 + assert inventory["indexes"]["neb_enzyme_to_files"]["BsaI"]["single_cutters"] == [ + "plasmid.gbk" + ] + + +def test_feature_type_and_function_are_separate() -> None: + roles = classify_feature_roles("CDS", {"label": ["FLAG-NLS"]}) + assert roles == ["coding_sequence", "epitope_tag", "localization_signal"] diff --git a/tools/inventory_cloning_sequences.py b/tools/inventory_cloning_sequences.py new file mode 100644 index 0000000..93bdb50 --- /dev/null +++ b/tools/inventory_cloning_sequences.py @@ -0,0 +1,66 @@ +#!/usr/bin/env python3 +"""Generate JSON indexes for cloning GenBank and FASTA attachments.""" + +from __future__ import annotations + +import argparse +import json +import sys +from pathlib import Path + +from lab_bench_2.cloning_inventory import build_cloning_inventory + + +def _arguments() -> argparse.Namespace: + parser = argparse.ArgumentParser( + description=( + "Inventory sequence files, GenBank features/primers, duplicates, and " + "New England Biolabs restriction-enzyme sites." + ) + ) + parser.add_argument("inputs", nargs="+", type=Path) + parser.add_argument( + "--root", type=Path, help="Root used to make paths in the JSON relative" + ) + parser.add_argument("--output", type=Path, required=True) + parser.add_argument( + "--no-enzymes", + action="store_true", + help="Skip the Biopython supplier-N enzyme catalog and site index", + ) + return parser.parse_args() + + +def main() -> int: + args = _arguments() + try: + inventory = build_cloning_inventory( + args.inputs, + root=args.root, + include_enzymes=not args.no_enzymes, + ) + args.output.parent.mkdir(parents=True, exist_ok=True) + args.output.write_text( + json.dumps(inventory, indent=2, sort_keys=True) + "\n", + encoding="utf-8", + ) + except (OSError, ValueError) as exc: + print(f"error: {exc}", file=sys.stderr) + return 2 + + summary = inventory["summary"] + print( + f"Inventoried {summary['parsed_file_count']}/" + f"{summary['discovered_file_count']} sequence files; " + f"{summary['parse_error_count']} parse error(s). Output: {args.output}" + ) + print( + f"No part features: " + f"{len(summary['files_with_no_part_features'])}; " + f"no primers: {len(summary['files_with_no_primers'])}" + ) + return 1 if summary["parse_error_count"] else 0 + + +if __name__ == "__main__": + raise SystemExit(main()) From dba650bde63d7f0e6a886a8ebb1caa0e61e2acd5 Mon Sep 17 00:00:00 2001 From: Divide-By-0 Date: Fri, 28 Aug 2026 00:42:39 -0700 Subject: [PATCH 2/3] Add provenance-aware cloning enrichment sources --- src/lab_bench_2/CLONING_SEQUENCE_DATA.md | 262 +++++ src/lab_bench_2/cloning_external_sources.py | 998 ++++++++++++++++++ src/lab_bench_2/cloning_inventory.py | 218 +++- src/lab_bench_2/plannotate_runner.py | 410 +++++++ .../test_cloning_external_sources.py | 361 +++++++ tests/lab_bench_2/test_cloning_inventory.py | 87 ++ tests/lab_bench_2/test_plannotate_runner.py | 118 +++ tools/annotate_cloning_sequences.py | 67 ++ tools/fetch_cloning_source_data.py | 93 ++ tools/inventory_cloning_sequences.py | 17 +- 10 files changed, 2625 insertions(+), 6 deletions(-) create mode 100644 src/lab_bench_2/CLONING_SEQUENCE_DATA.md create mode 100644 src/lab_bench_2/cloning_external_sources.py create mode 100644 src/lab_bench_2/plannotate_runner.py create mode 100644 tests/lab_bench_2/test_cloning_external_sources.py create mode 100644 tests/lab_bench_2/test_plannotate_runner.py create mode 100644 tools/annotate_cloning_sequences.py create mode 100644 tools/fetch_cloning_source_data.py diff --git a/src/lab_bench_2/CLONING_SEQUENCE_DATA.md b/src/lab_bench_2/CLONING_SEQUENCE_DATA.md new file mode 100644 index 0000000..fd58ffb --- /dev/null +++ b/src/lab_bench_2/CLONING_SEQUENCE_DATA.md @@ -0,0 +1,262 @@ +# Cloning sequence data and enrichment + +This document describes how the LAB-Bench 2 cloning attachments reach a model, +what is present in those files, and how to build separate, provenance-preserving +inventories from Biopython, REBASE, the iGEM Registry, and pLannotate. The audit +and live-source counts below were collected on 2026-08-27. + +## What cloning models receive + +LAB-Bench 2 has 14 cloning questions. At the dataset revision pinned by this +repository, their official attachment set contains 30 file instances: 29 +GenBank files and one FASTA file. Depending on evaluation mode, the same raw +source contents are injected into the prompt, attached as files, or copied into +the retrieval sandbox. + +The model does not currently receive a normalized parts list, a functional +feature inventory, a feature-to-file index, a primer index, an enzyme catalog, +or external annotations. The tools here generate those structures for analysis +and future dataset design; they do not silently change the existing evaluation. + +The local `~/.cache/labbench2/labbench2-data-public/cloning` snapshot contains +31 supported sequence files because one official GenBank file has a `(1)` copy. +It also contains one `.fa.orig` repair artifact, which this inventory ignores. +Those local artifacts are not additional official inputs. + +## Data layers + +Keep structural annotations, inferred function, and external matches separate: + +| Layer | Data | Trust and use | +| --- | --- | --- | +| Source GBK/FASTA | Sequence, coordinates, feature type, qualifiers, annotated primers | Ground truth for what the benchmark supplied | +| Biopython | Parsed records, feature/primer indexes, sequence hashes, restriction sites | Deterministic computation from the source file | +| REBASE | Current restriction-enzyme records with supplier codes | Live catalog evidence; supplier `N` is New England Biolabs | +| iGEM Registry | Specific published parts, descriptions, roles, sequences, and vocabulary | External cross-reference with explicit match evidence | +| pLannotate | SnapGene, FPbase, Swiss-Prot, and Rfam sequence-search hits | Sequence-derived enrichment candidate, kept in separate outputs | +| LLM | Prompt-ready functional description | Inference only; retain model, prompt, evidence, and confidence | + +For example, `CDS` says that a region encodes a protein or peptide. It does not +say whether that product is Cas9, a localization signal, an epitope tag, or a +selection marker. The inventory therefore preserves both dimensions: + +```json +{ + "feature_type": "CDS", + "label": "SV40 NLS", + "functional_roles": ["coding_sequence", "localization_signal"], + "functional_description": "SV40 NLS: protein-coding sequence, subcellular localization signal" +} +``` + +## Fetch the independent source data + +No Addgene login or API token is involved in this workflow. Fetch the public +vocabularies and current enzyme catalog once; subsequent runs use a validated +local cache: + +```bash +uv run --extra lab_bench_2 python tools/fetch_cloning_source_data.py \ + ~/.cache/labbench2/labbench2-data-public/cloning \ + --output cloning-external-sources.json +``` + +The fetcher is serial, rate-aware, retried, HTTPS-only, host-allowlisted, and +atomically cached. Each response records its URL, retrieval time, SHA-256, and +selected response headers. + +The live result used for this audit contains: + +- REBASE release 608, dated July 31, 2026: 234 restriction enzymes whose + supplier codes include `N` (New England Biolabs). +- Biopython 1.87 supplier `N`: 238 enzymes. It directly supports 233 of the 234 + live REBASE names. The spelling difference `CviKI-1` versus `CviKI_1` accounts + for one mismatch; five names remain only in Biopython's supplier set. +- iGEM Registry API 1.0: 88,418 published parts, 283 nondeprecated Sequence + Ontology roles, 329 categories, and 53 `//function/...` categories. + +The [REBASE monthly EMBOSS files](https://rebase.neb.com/rebase/rebase.files.html) +provide recognition and cut geometry separately from supplier records. This is +better evidence for current NEB availability than treating Biopython's bundled +supplier list as a purchasing catalog. Biopython remains the computation engine +for finding sites in each sequence. + +The [iGEM Registry API](https://api.registry.igem.org/docs) exposes published +part searches, individual JSON/FASTA/GenBank/SBOL records, roles, and categories +without authentication. Its one-shot `/v1/parts/export.fasta` endpoint currently +requires login, so this workflow fetches the complete small vocabularies and +uses targeted part lookups rather than pretending an anonymous bulk sequence +mirror exists. + +## Generate the source inventory + +```bash +uv run --extra lab_bench_2 python tools/inventory_cloning_sequences.py \ + ~/.cache/labbench2/labbench2-data-public/cloning \ + --root ~/.cache/labbench2/labbench2-data-public/cloning \ + --external-sources cloning-external-sources.json \ + --output cloning-sequence-inventory.json +``` + +The JSON contains: + +- A per-file and per-record inventory with file and sequence hashes, format, + length, topology, raw feature count, biological-part count, primer count, + parser warnings, and errors. +- A separate parts list for each non-source, non-primer feature: structural + type, coordinates, strand, raw qualifiers, normalized label, deterministic + functional roles, and a conservative prompt-ready description. +- A primer list with labels, coordinates, strand, and binding sequence in + 5-prime-to-3-prime orientation. +- Reverse indexes from feature label, functional role, primer label, sequence + hash, selected iGEM part/role, and restriction enzyme to files. +- Duplicate-sequence groups and explicit lists of files with no biological + parts or no annotated primers. +- Biopython's restriction metadata plus per-record zero-, single-, and + multi-cutter results. + +iGEM enrichment keeps three different concepts separate: + +- `feature_type` is GenBank syntax such as `CDS`, `promoter`, or + `protein_bind`. It is copied from the source and is never used as an iGEM + role. +- `functional_roles` are conservative, deterministic prompt terms such as + `epitope_tag`, `localization_signal`, or `selection_marker` inferred from the + source label and qualifiers. +- `specific_parts[].role` is the role attached by iGEM to one particular + published part record. It is shown only with that part's title, useful + description, URL, and nucleotide/protein match evidence. + +For each unique source feature signature, the fetcher searches the anonymous +published-part API by source label and exact source length. It retrieves the +best candidate details, checks the candidate's actual iGEM role against the +source-derived functional role, and auto-selects at most one part only when the +nucleotide sequence or translated peptide is exact. A small auditable alias +table resolves especially ambiguous canonical names such as `FLAG`, but the +same role and sequence-evidence gate still applies. Unverified search results +remain review candidates rather than selected cross-references. Exact source +label lookup against iGEM's role/category vocabulary is retained separately as +a low-confidence aid; GenBank structural types are not fed into that lookup. + +On `addgene-plasmid-181752-sequence-353936.gbk`, 12 of 20 unique feature +signatures have a selected specific iGEM part; because some signatures occur +more than once, this covers 13 of 21 source feature instances. For example, +source `FLAG` is a GenBank `CDS`, has local function `epitope_tag`, and maps to +[BBa_K4587111](https://registry.igem.org/parts/bba-k4587111), whose iGEM role is +`Tag` (`SO:0000324`). The source DNA and iGEM DNA are only 83.333% identical, +but both translate exactly to `DYKDDDDK`, so the match is protein-verified and +not falsely presented as a nucleotide-identical part. + +## Run pLannotate separately + +```bash +uv run --extra lab_bench_2 python tools/annotate_cloning_sequences.py \ + ~/.cache/labbench2/labbench2-data-public/cloning \ + --output-dir cloning-plannotate +``` + +If pLannotate is absent, the tool automatically creates a cached conda +environment pinned to pLannotate 2.0.0, downloads its checksum-validated +database bundle, and records the complete database manifest. Fast mode uses +SnapGene and FPbase; `--full` also searches Swiss-Prot and Rfam. Original +attachments are never modified. Generated GBK/CSV files and their hashes live +in a separate directory with a source-file manifest. Automatic setup requires +an installed conda, mamba, or micromamba executable; this machine already has +conda, so the validated run required no manual setup. + +The installed pLannotate v2 database used here was built July 5, 2026 and +records SnapGene, FPbase, UniProtKB/Swiss-Prot 2026_02, and Rfam 15.1 sources. +See the [pLannotate repository](https://github.com/mmcguffi/pLannotate) for its +search and nested-feature policies. + +Short peptide tags expose an important limitation. The selected Addgene +example contains a source-annotated FLAG feature with DNA +`GATTACAAAGACGATGACGATAAG`. pLannotate's SnapGene database contains a FLAG +reference with DNA `GATTACAAGGACGACGATGACAAG`. These 24 bp encodings are only +83.333% identical, below the configured SnapGene BLAST threshold of 95% +nucleotide identity, while both translate exactly to `DYKDDDDK`. Fast mode +therefore omits FLAG even though the peptide function is the same. The specific +iGEM lookup retains both DNA and translated-peptide evidence so this synonymous +codon case is visible rather than silently lost or falsely called +nucleotide-identical. + +Real validation on the sole featureless FASTA found one 16 bp SnapGene fragment +at 100% identity but only 3.67% reference coverage. pLannotate correctly marked +it as a fragment. That is not enough evidence to promote it into source truth, +and it illustrates why the CSV evidence, fragment flag, identity, and coverage +must remain visible to models and dataset authors. + +On a pET-28b GenBank attachment, the same fast pass found 20 candidate +annotations: 14 non-fragments and ten exact full-length matches. Exact calls +included the bacterial origin, f1 origin, rop, lacI promoter, T7 promoter and +terminator, T7 tag, lac operator, RBS, and thrombin site. Lower-coverage TcR, +AmpR-promoter, SET1, polyhistidine, and capTEV calls remained visibly marked as +fragments. This is the useful case for pLannotate: filling or checking a plasmid +parts inventory while preserving the match evidence. + +## Current attachment audit + +All 31 supported files in the local cache parse with Biopython. Ignoring the +local duplicate leaves the 30 official attachment instances. The sole file +without biological part annotations is: + +- `61e4b666-1ee5-4046-b304-d57e183c8593/Homo_sapiens_ENST00000612748_1_sequence.fa` + +The official attachment instances with no annotated primers are: + +- `4fb34135-1416-4b17-8bb9-c3a4dd3a2da7/NM_002478.5.gb` +- `4fb34135-1416-4b17-8bb9-c3a4dd3a2da7/pivt-fluc-bxb.gb` +- `5e7bf2b5-5d04-48c6-83d4-38d97ea8542c/pMBP-bdSUMO-dArc1.gb` +- `61e4b666-1ee5-4046-b304-d57e183c8593/Homo_sapiens_ENST00000612748_1_sequence.fa` +- `a4bf037c-2477-4cca-9ca3-12c5ee63c44f/npas4-201-enst00000311034.gb` +- `ad7daee9-2b65-43dc-906b-42ecb07f58b1/pet-28b.gb` +- `ae62bcdb-197b-4815-991f-cb7a9c151ff6/pet-28b.gb` +- `dff28bd4-89ea-401e-a3c1-ce8733512b79/pMBP-bdSUMO-dArc1.gb` +- `fb8fc27d-592a-40e8-a65f-9e1a60b7a708/sorcs2-201-ensmust00000037370.gb` + +That is nine attachment instances but seven unique filenames. Across the 29 +official GenBank instances, 21 contain `primer_bind` features. The source +records contain 313 primer annotations and 50 distinct original labels before +normalized aliases are collapsed. + +Biopython warns that six local files encode origin-wrapping locations in a +nonstandard form and repairs them while parsing. The inventory captures all ten +warnings; it does not rewrite the benchmark inputs. + +## Model-facing inventory and question difficulty + +For a cloning prompt, provide only the slice needed for the question: + +1. Candidate source files and exact source features. +2. Functional-role candidates with their provenance and review status. +3. Relevant primers or the explicit fact that no primers are annotated. +4. Relevant enzyme sites and whether the enzyme is in the live REBASE/NEB set. +5. External sequence-search evidence only when it materially resolves an + annotation gap. + +Represent difficulty as an evidence graph. Nodes include requested functions, +files, features, fragments, primers, enzymes, edits, junctions, and the final +construct. Typed edges record the necessary connection: prompt function to +role, role to feature, feature to boundary, fragment to junction, junction to +primer or edit, and enzyme to compatible recognition site/overhang. + +A useful raw vector is: + +```json +{ + "inventory_gaps": 0, + "semantic_resolution_edges": 2, + "candidate_branches": 1, + "fragment_selection_edges": 2, + "junctions": 1, + "primers_to_design": 2, + "sequence_edits": 0, + "enzyme_constraints": 3, + "external_provenance_hops": 0, + "required_connection_count": 10 +} +``` + +The connection count should describe the smallest validated reasoning path, +not every possible option. Keep the raw vector even if it is later mapped to a +coarse D0-D4 label, so difficulty changes remain explainable. diff --git a/src/lab_bench_2/cloning_external_sources.py b/src/lab_bench_2/cloning_external_sources.py new file mode 100644 index 0000000..a692e9f --- /dev/null +++ b/src/lab_bench_2/cloning_external_sources.py @@ -0,0 +1,998 @@ +"""Fetch versioned public vocabularies used to enrich cloning inventories.""" + +from __future__ import annotations + +import hashlib +import json +import re +import shutil +import subprocess +import tempfile +import time +from dataclasses import dataclass +from datetime import UTC, datetime +from difflib import SequenceMatcher +from pathlib import Path +from typing import Any, Callable, Mapping, Sequence +from urllib.error import HTTPError, URLError +from urllib.parse import urlencode, urlsplit +from urllib.request import Request, urlopen + +from lab_bench_2.cloning_inventory import normalize_label + +REBASE_BASE_URL = "https://rebase.neb.com/rebase" +IGEM_API_BASE_URL = "https://api.registry.igem.org" +IGEM_REGISTRY_PART_URL = "https://registry.igem.org/parts" +RETRYABLE_HTTP_STATUSES = {429, 500, 502, 503, 504} +REBASE_ENZYME_FIELD_COUNT = 9 +REBASE_REFERENCE_HEADER_FIELD_COUNT = 7 +MAX_IGEM_PART_DETAILS_PER_FEATURE = 8 + +# A tiny, auditable override table is appropriate for canonical aliases whose +# Registry search has many synonymous entries. The evidence block still records +# whether the source nucleotide or translated peptide actually matches the part. +IGEM_CANONICAL_PART_SLUGS = {"flag": "bba-k4587111"} + +EXPECTED_IGEM_ROLE_ACCESSIONS = { + "affinity_tag": ("SO:0000324", "SO:0000807"), + "coding_sequence": ("SO:0000316",), + "enhancer": ("SO:0000165", "IGEM:0000006"), + "epitope_tag": ("SO:0000324", "SO:0000807"), + "genome_editor": ("SO:0000316",), + "intron_or_splicing": ("SO:0000188", "IGEM:0000036"), + "localization_signal": ("SO:0001527", "IGEM:0000010", "SO:0000417"), + "origin_of_replication": ("SO:0000296",), + "promoter": ("SO:0000167", "IGEM:0000006"), + "protein_binding_site": ("SO:0000057", "SO:0000410", "IGEM:0000006"), + "selection_marker": ("SO:0002232", "SO:0000316"), + "terminator_or_polya": ("SO:0000141", "SO:0000551"), +} + + +class ExternalSourceError(RuntimeError): + """Raised when a cloning-data source cannot be fetched or parsed.""" + + +@dataclass(frozen=True) +class CachedResponse: + """A cached or newly fetched HTTP response with immutable provenance.""" + + body: bytes + url: str + fetched_at: str + sha256: str + headers: dict[str, str] + cache_status: str + + def provenance(self) -> dict[str, Any]: + """Return the serializable provenance fields for this response.""" + return { + "url": self.url, + "fetched_at": self.fetched_at, + "sha256": self.sha256, + "headers": self.headers, + "cache_status": self.cache_status, + } + + +OpenUrl = Callable[..., Any] +SleepFunction = Callable[[float], None] +MonotonicFunction = Callable[[], float] + + +class CachedHttpClient: + """Serial HTTPS client with validation, retry, and a content-addressed cache.""" + + def __init__( + self, + cache_dir: Path, + *, + allowed_hosts: Sequence[str], + min_delay: float = 0.2, + retries: int = 3, + timeout: float = 60.0, + open_url: OpenUrl | None = None, + sleep: SleepFunction = time.sleep, + monotonic: MonotonicFunction = time.monotonic, + ) -> None: + if min_delay < 0 or retries < 0 or timeout <= 0: + raise ValueError("Require nonnegative delay/retries and positive timeout") + self.cache_dir = cache_dir.expanduser().resolve() + self.allowed_hosts = {host.casefold() for host in allowed_hosts} + self.min_delay = min_delay + self.retries = retries + self.timeout = timeout + self._open_url = open_url or _default_open_url + self._sleep = sleep + self._monotonic = monotonic + self._last_request_at: float | None = None + + def get(self, url: str, *, cache_key: str, refresh: bool = False) -> CachedResponse: + """Fetch one allowlisted HTTPS URL, reusing a validated cache by default.""" + self._validate_url(url) + body_path, metadata_path = self._cache_paths(cache_key) + if not refresh: + cached = self._read_cache(body_path, metadata_path, expected_url=url) + if cached is not None: + return cached + + last_error: Exception | None = None + for attempt in range(self.retries + 1): + self._wait_for_rate_limit() + request = Request( + url, + headers={ + "Accept": "application/json, text/plain, */*", + "User-Agent": "lab-bench-cloning-source-inventory/1.0", + }, + ) + self._last_request_at = self._monotonic() + try: + with self._open_url(request, timeout=self.timeout) as response: + body = response.read() + headers = _selected_headers(response.headers) + final_url = str(response.geturl()) + self._validate_url(final_url) + if not body: + raise ExternalSourceError(f"Empty response from {url}") + result = CachedResponse( + body=body, + url=final_url, + fetched_at=datetime.now(UTC).isoformat(), + sha256=hashlib.sha256(body).hexdigest(), + headers=headers, + cache_status="fetched", + ) + self._write_cache(body_path, metadata_path, result, request_url=url) + return result + except HTTPError as exc: + last_error = exc + if exc.code not in RETRYABLE_HTTP_STATUSES or attempt >= self.retries: + break + self._sleep(_retry_delay(exc.headers, attempt)) + except (URLError, TimeoutError, OSError) as exc: + last_error = exc + if attempt >= self.retries: + break + self._sleep(float(2**attempt)) + raise ExternalSourceError( + f"Could not fetch {url}: {last_error}" + ) from last_error + + def get_json( + self, url: str, *, cache_key: str, refresh: bool = False + ) -> tuple[Any, CachedResponse]: + """Fetch and decode a JSON response.""" + response = self.get(url, cache_key=cache_key, refresh=refresh) + try: + payload = json.loads(response.body) + except (UnicodeDecodeError, json.JSONDecodeError) as exc: + raise ExternalSourceError(f"Invalid JSON from {url}") from exc + return payload, response + + def _validate_url(self, url: str) -> None: + parsed = urlsplit(url) + if ( + parsed.scheme != "https" + or (parsed.hostname or "").casefold() not in self.allowed_hosts + ): + raise ExternalSourceError(f"Refusing non-allowlisted source URL: {url}") + + def _cache_paths(self, cache_key: str) -> tuple[Path, Path]: + safe_key = re.sub(r"[^A-Za-z0-9_.-]+", "_", cache_key).strip("_") + if not safe_key: + raise ValueError("cache_key must contain a filename-safe character") + source_hash = hashlib.sha256(cache_key.encode()).hexdigest()[:12] + stem = f"{safe_key[:80]}-{source_hash}" + return self.cache_dir / f"{stem}.body", self.cache_dir / f"{stem}.json" + + def _read_cache( + self, body_path: Path, metadata_path: Path, *, expected_url: str + ) -> CachedResponse | None: + try: + body = body_path.read_bytes() + metadata = json.loads(metadata_path.read_text(encoding="utf-8")) + digest = hashlib.sha256(body).hexdigest() + if digest != metadata["sha256"] or not body: + return None + url = str(metadata["url"]) + request_url = str(metadata.get("request_url") or url) + if request_url != expected_url: + return None + self._validate_url(url) + headers = { + str(key): str(value) + for key, value in dict(metadata.get("headers") or {}).items() + } + return CachedResponse( + body=body, + url=url, + fetched_at=str(metadata["fetched_at"]), + sha256=digest, + headers=headers, + cache_status="cache-hit", + ) + except (OSError, KeyError, TypeError, ValueError, json.JSONDecodeError): + return None + + def _write_cache( + self, + body_path: Path, + metadata_path: Path, + result: CachedResponse, + *, + request_url: str, + ) -> None: + self.cache_dir.mkdir(parents=True, exist_ok=True) + body_temp = body_path.with_suffix(body_path.suffix + ".tmp") + metadata_temp = metadata_path.with_suffix(metadata_path.suffix + ".tmp") + body_temp.write_bytes(result.body) + metadata = {**result.provenance(), "request_url": request_url} + metadata_temp.write_text( + json.dumps(metadata, indent=2, sort_keys=True) + "\n", encoding="utf-8" + ) + body_temp.replace(body_path) + metadata_temp.replace(metadata_path) + + def _wait_for_rate_limit(self) -> None: + if self._last_request_at is None: + return + remaining = self.min_delay - (self._monotonic() - self._last_request_at) + if remaining > 0: + self._sleep(remaining) + + +def fetch_rebase_neb_catalog( + cache_dir: Path, + *, + refresh: bool = False, + http: CachedHttpClient | None = None, +) -> dict[str, Any]: + """Fetch the current REBASE restriction catalog sold by NEB (supplier N).""" + client = http or CachedHttpClient( + cache_dir / "rebase", + allowed_hosts=["rebase.neb.com"], + min_delay=0.2, + ) + source_responses = { + name: client.get( + f"{REBASE_BASE_URL}/{link}", cache_key=f"rebase-{name}", refresh=refresh + ) + for name, link in ( + ("enzymes", "link_emboss_e"), + ("references", "link_emboss_r"), + ("suppliers", "link_emboss_s"), + ) + } + texts = { + name: response.body.decode("utf-8") + for name, response in source_responses.items() + } + releases = {_parse_rebase_release(text) for text in texts.values()} + if len(releases) != 1: + raise ExternalSourceError( + f"REBASE source files disagree on release: {releases}" + ) + release = releases.pop() + suppliers = parse_rebase_suppliers(texts["suppliers"]) + if suppliers.get("N") != "New England Biolabs": + raise ExternalSourceError("REBASE supplier N is not New England Biolabs") + enzymes = parse_rebase_enzymes(texts["enzymes"]) + references = parse_rebase_references(texts["references"]) + biopython_names, biopython_version = _biopython_neb_names() + + catalog = [] + for name, reference in references.items(): + supplier_codes = reference["supplier_codes"] + enzyme = enzymes.get(name) + if "N" not in supplier_codes or enzyme is None: + continue + catalog.append( + { + **enzyme, + "organism": reference["organism"], + "isoschizomers": reference["isoschizomers"], + "methylation": reference["methylation"], + "source": reference["source"], + "supplier_codes": supplier_codes, + "supplier_names": [ + suppliers[code] for code in supplier_codes if code in suppliers + ], + "reference_count": reference["reference_count"], + "analysis_supported_by_biopython": name in biopython_names, + } + ) + catalog.sort(key=lambda item: str(item["name"])) + rebase_names = {str(item["name"]) for item in catalog} + return { + "schema_version": 1, + "source": { + "name": "REBASE", + "release": release, + "release_date": _parse_rebase_release_date(texts["enzymes"]), + "supplier_code": "N", + "supplier_name": suppliers["N"], + "files": { + name: response.provenance() + for name, response in source_responses.items() + }, + }, + "enzyme_count": len(catalog), + "enzymes": catalog, + "biopython_comparison": { + "biopython_version": biopython_version, + "supplier_n_count": len(biopython_names), + "overlap_count": len(rebase_names & biopython_names), + "only_in_live_rebase": sorted(rebase_names - biopython_names), + "only_in_biopython": sorted(biopython_names - rebase_names), + }, + } + + +class IgemRegistryClient: + """Read the anonymous, published subset of the iGEM Registry API.""" + + def __init__( + self, + cache_dir: Path, + *, + refresh: bool = False, + http: CachedHttpClient | None = None, + ) -> None: + self.refresh = refresh + self.http = http or CachedHttpClient( + cache_dir / "igem", + allowed_hosts=["api.registry.igem.org"], + min_delay=2.1, + ) + + def fetch_vocabulary( + self, feature_queries: Sequence[Mapping[str, Any]] = () + ) -> dict[str, Any]: + """Fetch current vocabularies and evidence-ranked public part matches.""" + openapi, openapi_response = self._get_json("/docs-json", "openapi") + roles, role_responses = self._get_paginated( + "/v1/parts/roles", {"deprecated": "false"}, "roles" + ) + categories, category_responses = self._get_paginated( + "/v1/parts/categories", {}, "categories" + ) + part_page, part_response = self._get_json( + "/v1/parts", "parts-count", {"page": 1, "pageSize": 1} + ) + if not isinstance(openapi, dict) or not isinstance(part_page, dict): + raise ExternalSourceError("iGEM returned an unexpected schema") + role_aliases = _igem_role_alias_index(roles) + category_aliases = _igem_category_alias_index(categories) + part_matches, part_match_responses, part_match_errors = ( + self._fetch_specific_part_matches(feature_queries, roles) + ) + info = openapi.get("info") or {} + return { + "schema_version": 1, + "source": { + "name": "iGEM Registry", + "api_base_url": IGEM_API_BASE_URL, + "api_title": str(info.get("title") or ""), + "api_version": str(info.get("version") or ""), + "access": "anonymous published records", + "openapi": openapi_response.provenance(), + "pages": [ + response.provenance() + for response in [*role_responses, *category_responses] + ], + "parts_count_request": part_response.provenance(), + "part_match_requests": [ + response.provenance() for response in part_match_responses + ], + "part_match_errors": part_match_errors, + }, + "summary": { + "published_part_count": int(part_page.get("total") or 0), + "role_count": len(roles), + "category_count": len(categories), + "functional_category_count": sum( + _category_path(item).startswith("//function/") + for item in categories + ), + "feature_query_count": len(feature_queries), + "specific_part_match_count": sum( + any(bool(part.get("selected")) for part in matches) + for matches in part_matches.values() + ), + "part_match_error_count": len(part_match_errors), + }, + "roles": sorted(roles, key=lambda item: str(item.get("accession") or "")), + "categories": sorted( + categories, key=lambda item: _category_path(item).casefold() + ), + "indexes": { + "role_alias_to_terms": role_aliases, + "category_alias_to_terms": category_aliases, + "feature_signature_to_part_matches": part_matches, + }, + } + + def _fetch_specific_part_matches( + self, + feature_queries: Sequence[Mapping[str, Any]], + roles: Sequence[Mapping[str, Any]], + ) -> tuple[ + dict[str, list[dict[str, Any]]], list[CachedResponse], list[dict[str, str]] + ]: + del roles + matches: dict[str, list[dict[str, Any]]] = {} + responses: list[CachedResponse] = [] + errors: list[dict[str, str]] = [] + seen_signatures: set[str] = set() + for query in feature_queries: + signature = str(query.get("feature_signature") or "") + label = str(query.get("normalized_label") or "") + if not signature or not label or signature in seen_signatures: + continue + seen_signatures.add(signature) + canonical_slug = IGEM_CANONICAL_PART_SLUGS.get(label) + if canonical_slug: + try: + payload, response = self._get_json( + f"/v1/parts/slugs/{canonical_slug}", + f"part-{canonical_slug}", + ) + except ExternalSourceError as exc: + errors.append({"feature_signature": signature, "error": str(exc)}) + continue + responses.append(response) + if isinstance(payload, Mapping): + matches[signature] = [ + _igem_part_candidate( + payload, + query, + canonical=True, + expected_accessions=_expected_role_accessions(query), + provenance=response.provenance(), + ) + ] + continue + + expected_accessions = _expected_role_accessions(query) + length = int(query.get("nucleotide_length") or 0) + parameters: dict[str, Any] = { + "page": 1, + "pageSize": 100, + "search": str(query.get("label") or label), + } + if length: + parameters.update({"minLength": length, "maxLength": length}) + try: + payload, search_response = self._get_json( + "/v1/parts", + "part-search-" + f"{hashlib.sha256(signature.encode()).hexdigest()[:24]}", + parameters, + ) + except ExternalSourceError as exc: + errors.append({"feature_signature": signature, "error": str(exc)}) + continue + responses.append(search_response) + raw_parts = payload.get("data", []) if isinstance(payload, Mapping) else [] + candidates = [part for part in raw_parts if isinstance(part, Mapping)] + candidates.sort( + key=lambda part: _igem_part_preliminary_score( + part, query, expected_accessions + ), + reverse=True, + ) + evaluated: list[dict[str, Any]] = [] + for part in candidates[:MAX_IGEM_PART_DETAILS_PER_FEATURE]: + slug = str(part.get("slug") or "") + if not slug: + continue + try: + detail, detail_response = self._get_json( + f"/v1/parts/slugs/{slug}", f"part-{slug}" + ) + except ExternalSourceError as exc: + errors.append( + { + "feature_signature": signature, + "part_slug": slug, + "error": str(exc), + } + ) + continue + responses.append(detail_response) + if isinstance(detail, Mapping): + evaluated.append( + _igem_part_candidate( + detail, + query, + canonical=False, + expected_accessions=expected_accessions, + provenance=detail_response.provenance(), + ) + ) + evaluated.sort(key=_igem_part_sort_key, reverse=True) + selected = next( + ( + part + for part in evaluated + if part["evidence"]["expected_functional_role"] + and ( + part["evidence"]["nucleotide_exact"] + or part["evidence"]["translated_peptide_exact"] + ) + ), + None, + ) + if selected is not None: + selected["selected"] = True + if evaluated: + matches[signature] = evaluated[:3] + return matches, responses, errors + + def _get_paginated( + self, route: str, parameters: Mapping[str, Any], cache_prefix: str + ) -> tuple[list[dict[str, Any]], list[CachedResponse]]: + page = 1 + total: int | None = None + items: list[dict[str, Any]] = [] + responses: list[CachedResponse] = [] + while total is None or len(items) < total: + payload, response = self._get_json( + route, + f"{cache_prefix}-page-{page}", + {**parameters, "page": page, "pageSize": 100}, + ) + if not isinstance(payload, dict) or not isinstance( + payload.get("data"), list + ): + raise ExternalSourceError(f"Unexpected paginated iGEM result: {route}") + page_items = payload["data"] + typed_items = [item for item in page_items if isinstance(item, dict)] + if len(typed_items) != len(page_items): + raise ExternalSourceError(f"Non-object item in iGEM result: {route}") + items.extend(typed_items) + responses.append(response) + total = int(payload.get("total") or 0) + if not page_items: + break + page += 1 + return items, responses + + def _get_json( + self, + route: str, + cache_key: str, + parameters: Mapping[str, Any] | None = None, + ) -> tuple[Any, CachedResponse]: + query = f"?{urlencode(parameters, doseq=True)}" if parameters else "" + return self.http.get_json( + f"{IGEM_API_BASE_URL}{route}{query}", + cache_key=f"igem-{cache_key}", + refresh=self.refresh, + ) + + +def fetch_cloning_external_sources( + cache_dir: Path, + *, + refresh: bool = False, + include_rebase: bool = True, + include_igem: bool = True, + feature_queries: Sequence[Mapping[str, Any]] = (), +) -> dict[str, Any]: + """Fetch the independent external catalogs used by cloning inventories.""" + result: dict[str, Any] = { + "schema_version": 1, + "generated_at": datetime.now(UTC).isoformat(), + "sources": {}, + } + sources = result["sources"] + if not isinstance(sources, dict): + raise AssertionError("sources must be a mapping") + if include_rebase: + sources["rebase_neb"] = fetch_rebase_neb_catalog(cache_dir, refresh=refresh) + if include_igem: + sources["igem_registry"] = IgemRegistryClient( + cache_dir, refresh=refresh + ).fetch_vocabulary(feature_queries) + return result + + +def _expected_role_accessions(query: Mapping[str, Any]) -> tuple[str, ...]: + accessions: set[str] = set() + functional_roles = query.get("functional_roles") or [] + if isinstance(functional_roles, Sequence) and not isinstance( + functional_roles, (str, bytes) + ): + for role in functional_roles: + accessions.update(EXPECTED_IGEM_ROLE_ACCESSIONS.get(str(role), ())) + return tuple(sorted(accessions)) + + +def _igem_part_preliminary_score( + part: Mapping[str, Any], + query: Mapping[str, Any], + expected_accessions: Sequence[str], +) -> float: + label = normalize_label(str(query.get("label") or "")) + title = normalize_label(str(part.get("title") or "")) + name = normalize_label(str(part.get("name") or "")) + source_length = int(query.get("nucleotide_length") or 0) + candidate_length = int(part.get("sequenceLength") or 0) + role = part.get("role") + role_accession = ( + str(role.get("accession") or "") if isinstance(role, Mapping) else "" + ) + score = SequenceMatcher(None, label, title).ratio() * 40 + if label and label in {title, name}: + score += 80 + elif label and all(token in title.split() for token in label.split()): + score += 35 + if source_length and source_length == candidate_length: + score += 40 + if role_accession in expected_accessions: + score += 30 + score += min(float(part.get("usageCount") or 0), 20.0) / 4 + return score + + +def _igem_part_candidate( + part: Mapping[str, Any], + query: Mapping[str, Any], + *, + canonical: bool, + expected_accessions: Sequence[str], + provenance: Mapping[str, Any], +) -> dict[str, Any]: + from Bio.Seq import Seq + + source_sequence = str(query.get("nucleotide_sequence") or "").upper() + source_translation = _normalize_translation(str(query.get("translation") or "")) + candidate_sequence = str(part.get("sequence") or "").upper() + candidate_translation = "" + if candidate_sequence and len(candidate_sequence) % 3 == 0: + candidate_translation = _normalize_translation( + str(Seq(candidate_sequence).translate()) # type: ignore[no-untyped-call] + ) + nucleotide_exact = bool(source_sequence) and candidate_sequence in { + source_sequence, + str(Seq(source_sequence).reverse_complement()), # type: ignore[no-untyped-call] + } + translated_exact = bool(source_translation) and ( + candidate_translation == source_translation + ) + identity = _same_length_identity(source_sequence, candidate_sequence) + role = part.get("role") + role_mapping = role if isinstance(role, Mapping) else {} + role_accession = str(role_mapping.get("accession") or "") + slug = str(part.get("slug") or "") + evidence = { + "canonical_alias": canonical, + "nucleotide_exact": nucleotide_exact, + "same_length_nucleotide_identity_percent": identity, + "source_translation": source_translation, + "candidate_translation": candidate_translation, + "translated_peptide_exact": translated_exact, + "expected_functional_role": role_accession in expected_accessions, + "source_feature_type_used_as_role": False, + } + selected = ( + canonical + and role_accession in expected_accessions + and (nucleotide_exact or translated_exact) + ) + return { + "uuid": str(part.get("uuid") or ""), + "name": str(part.get("name") or ""), + "slug": slug, + "url": f"{IGEM_REGISTRY_PART_URL}/{slug}", + "title": str(part.get("title") or ""), + "description": str(part.get("description") or ""), + "role": { + "uuid": str(role_mapping.get("uuid") or ""), + "accession": role_accession, + "label": str(role_mapping.get("label") or ""), + "definition": str(role_mapping.get("definition") or ""), + "source": str(role_mapping.get("source") or ""), + }, + "sequence_length": len(candidate_sequence), + "sequence_sha256": ( + hashlib.sha256(candidate_sequence.encode("ascii")).hexdigest() + if candidate_sequence + else "" + ), + "usage_count": int(part.get("usageCount") or 0), + "selected": selected, + "review_required": not (nucleotide_exact or translated_exact), + "selection_method": ( + "canonical source-label alias, expected iGEM role, and " + "sequence/translation evidence" + if canonical + else "ranked public search with expected iGEM role and " + "sequence/translation evidence" + ), + "evidence": evidence, + "source": {str(key): value for key, value in provenance.items()}, + } + + +def _normalize_translation(sequence: str) -> str: + return re.sub(r"[^A-Z*]", "", sequence.upper()).rstrip("*") + + +def _same_length_identity(left: str, right: str) -> float | None: + if not left or len(left) != len(right): + return None + return round( + 100 * sum(a == b for a, b in zip(left, right, strict=True)) / len(left), 3 + ) + + +def _igem_part_sort_key( + part: Mapping[str, Any], +) -> tuple[int, int, int, int, float, int]: + evidence = part.get("evidence") + evidence_mapping = evidence if isinstance(evidence, Mapping) else {} + identity = evidence_mapping.get("same_length_nucleotide_identity_percent") + return ( + int(bool(evidence_mapping.get("canonical_alias"))), + int(bool(evidence_mapping.get("expected_functional_role"))), + int(bool(evidence_mapping.get("nucleotide_exact"))), + int(bool(evidence_mapping.get("translated_peptide_exact"))), + float(identity) if isinstance(identity, (int, float)) else -1.0, + int(part.get("usage_count") or 0), + ) + + +def parse_rebase_enzymes(text: str) -> dict[str, dict[str, Any]]: + """Parse the REBASE EMBOSS enzyme-pattern file.""" + enzymes: dict[str, dict[str, Any]] = {} + for line in _data_lines(text): + fields = line.split() + if len(fields) != REBASE_ENZYME_FIELD_COUNT: + raise ExternalSourceError(f"Malformed REBASE enzyme row: {line}") + name, recognition_site = fields[:2] + try: + length, cut_count, blunt, cut_5, cut_3, second_5, second_3 = map( + int, fields[2:] + ) + except ValueError as exc: + raise ExternalSourceError(f"Malformed REBASE cut geometry: {line}") from exc + enzymes[name] = { + "name": name, + "recognition_site": recognition_site.upper(), + "recognition_length": length, + "cut_count": cut_count, + "blunt": bool(blunt), + "cut_5prime": cut_5, + "cut_3prime": cut_3, + "second_cut_5prime": second_5, + "second_cut_3prime": second_3, + } + return enzymes + + +def parse_rebase_suppliers(text: str) -> dict[str, str]: + """Parse the REBASE EMBOSS supplier file.""" + suppliers: dict[str, str] = {} + for line in _data_lines(text): + code, separator, name = line.partition(" ") + if len(code) != 1 or not separator or not name.strip(): + raise ExternalSourceError(f"Malformed REBASE supplier row: {line}") + suppliers[code] = name.strip() + return suppliers + + +def parse_rebase_references(text: str) -> dict[str, dict[str, Any]]: + """Parse supplier and source fields from REBASE EMBOSS reference records.""" + records: dict[str, dict[str, Any]] = {} + blocks = re.split(r"^//\s*$", text, flags=re.MULTILINE) + for block in blocks: + lines = [ + line.rstrip() for line in block.splitlines() if not line.startswith("#") + ] + while lines and not lines[0].strip(): + lines.pop(0) + if not lines: + continue + if len(lines) < REBASE_REFERENCE_HEADER_FIELD_COUNT: + raise ExternalSourceError(f"Malformed REBASE reference entry: {lines[0]}") + try: + reference_count = int(lines[6].strip()) + except ValueError as exc: + raise ExternalSourceError( + f"Malformed REBASE reference count for {lines[0]}" + ) from exc + supplier_codes = list(lines[5].strip()) + records[lines[0].strip()] = { + "organism": lines[1].strip(), + "isoschizomers": [ + item.strip() for item in lines[2].split(",") if item.strip() + ], + "methylation": lines[3].strip(), + "source": lines[4].strip(), + "supplier_codes": supplier_codes, + "reference_count": reference_count, + } + return records + + +def _data_lines(text: str) -> list[str]: + return [ + line.strip() + for line in text.splitlines() + if line.strip() and not line.startswith("#") + ] + + +def _parse_rebase_release(text: str) -> str: + match = re.search(r"REBASE version\s+(\d+)", text) + if not match: + raise ExternalSourceError("REBASE file has no release number") + return match.group(1) + + +def _parse_rebase_release_date(text: str) -> str: + match = re.search(r"\b([A-Z][a-z]{2}\s+\d{1,2}\s+\d{4})\b", text) + return match.group(1) if match else "unknown" + + +def _biopython_neb_names() -> tuple[set[str], str]: + try: + from Bio import __version__ as biopython_version + from Bio.Restriction import Restriction_Dictionary + except ModuleNotFoundError: + return set(), "not installed" + names = {str(name) for name in Restriction_Dictionary.suppliers["N"][1]} + return names, biopython_version + + +def _igem_role_alias_index( + roles: Sequence[Mapping[str, Any]], +) -> dict[str, list[dict[str, str]]]: + aliases: dict[str, dict[str, dict[str, str]]] = {} + for role in roles: + term = { + "accession": str(role.get("accession") or ""), + "label": str(role.get("label") or ""), + "definition": str(role.get("definition") or ""), + "source": str(role.get("source") or ""), + } + raw_aliases = [term["label"], *list(role.get("synonyms") or [])] + for raw_alias in raw_aliases: + alias = normalize_label(str(raw_alias)) + if alias: + aliases.setdefault(alias, {})[term["accession"]] = term + return { + alias: sorted(terms.values(), key=lambda term: term["accession"]) + for alias, terms in sorted(aliases.items()) + } + + +def _igem_category_alias_index( + categories: Sequence[Mapping[str, Any]], +) -> dict[str, list[dict[str, str]]]: + aliases: dict[str, dict[str, dict[str, str]]] = {} + for category in categories: + category_path = _category_path(category) + term = { + "uuid": str(category.get("uuid") or ""), + "path": category_path, + "description": str(category.get("description") or ""), + } + raw_aliases = [category_path, category_path.rsplit("/", 1)[-1]] + for raw_alias in raw_aliases: + alias = normalize_label(raw_alias) + if alias: + aliases.setdefault(alias, {})[term["uuid"]] = term + return { + alias: sorted(terms.values(), key=lambda term: term["path"]) + for alias, terms in sorted(aliases.items()) + } + + +def _category_path(category: Mapping[str, Any]) -> str: + return str(category.get("value") or category.get("label") or "") + + +def _selected_headers(headers: Mapping[str, Any]) -> dict[str, str]: + selected = {} + for key in ( + "Content-Type", + "Content-Length", + "ETag", + "Last-Modified", + "X-RateLimit-Limit-Short", + "X-RateLimit-Limit-Medium", + "X-RateLimit-Limit-Large", + ): + value = headers.get(key) + if value is not None: + selected[key.casefold()] = str(value) + return selected + + +class _CurlResponse: + def __init__(self, body: bytes, headers: dict[str, str], final_url: str) -> None: + self._body = body + self.headers = headers + self._final_url = final_url + + def __enter__(self) -> _CurlResponse: + return self + + def __exit__(self, *args: object) -> None: + return None + + def read(self) -> bytes: + return self._body + + def geturl(self) -> str: + return self._final_url + + +def _default_open_url(request: Request, *, timeout: float) -> Any: + curl = shutil.which("curl") + if curl is None: + return urlopen(request, timeout=timeout) + with tempfile.TemporaryDirectory(prefix="lab-bench-source-fetch-") as temp_dir: + temp_path = Path(temp_dir) + body_path = temp_path / "body" + headers_path = temp_path / "headers" + command = [ + curl, + "--fail-with-body", + "--silent", + "--show-error", + "--location", + "--max-time", + str(timeout), + "--dump-header", + str(headers_path), + "--output", + str(body_path), + "--write-out", + "%{url_effective}", + ] + for key, value in request.header_items(): + command.extend(["--header", f"{key}: {value}"]) + command.append(request.full_url) + try: + result = subprocess.run( + command, + check=False, + capture_output=True, + text=True, + timeout=timeout + 5, + ) + except (OSError, subprocess.TimeoutExpired) as exc: + raise URLError(str(exc)) from exc + if result.returncode != 0: + raise URLError((result.stderr or "curl request failed").strip()) + body = body_path.read_bytes() + headers = _parse_curl_headers(headers_path.read_text(errors="replace")) + return _CurlResponse(body, headers, result.stdout.strip() or request.full_url) + + +def _parse_curl_headers(text: str) -> dict[str, str]: + blocks = [block for block in re.split(r"\r?\n\r?\n", text) if block.strip()] + header_block = next( + (block for block in reversed(blocks) if block.startswith("HTTP/")), "" + ) + headers: dict[str, str] = {} + for line in header_block.splitlines()[1:]: + key, separator, value = line.partition(":") + if separator: + headers[key.strip()] = value.strip() + return headers + + +def _retry_delay(headers: Any, attempt: int) -> float: + retry_after = headers.get("Retry-After") + if retry_after is not None: + try: + return max(0.0, min(float(retry_after), 120.0)) + except ValueError: + pass + return float(2**attempt) diff --git a/src/lab_bench_2/cloning_inventory.py b/src/lab_bench_2/cloning_inventory.py index 165e9d3..093ea14 100644 --- a/src/lab_bench_2/cloning_inventory.py +++ b/src/lab_bench_2/cloning_inventory.py @@ -3,6 +3,7 @@ from __future__ import annotations import hashlib +import os import re import warnings from collections import defaultdict @@ -137,6 +138,7 @@ def build_cloning_inventory( *, root: Path | None = None, include_enzymes: bool = True, + external_sources: Mapping[str, Any] | None = None, ) -> dict[str, Any]: """Parse sequence attachments and build file, feature, primer, and enzyme indexes.""" files = discover_sequence_files(inputs) @@ -148,6 +150,9 @@ def build_cloning_inventory( enzyme_index: defaultdict[str, dict[str, list[str]]] = defaultdict( lambda: {"zero_cutters": [], "single_cutters": [], "multi_cutters": []} ) + external_role_index: defaultdict[str, set[str]] = defaultdict(set) + external_category_index: defaultdict[str, set[str]] = defaultdict(set) + external_part_index: defaultdict[str, set[str]] = defaultdict(set) entries: list[dict[str, Any]] = [] errors: list[dict[str, str]] = [] enzyme_catalog = _neb_enzyme_catalog() if include_enzymes else [] @@ -182,6 +187,27 @@ def build_cloning_inventory( ) enzyme_index[name][category].append(display_path) + external_enrichment = _apply_external_enrichment(entries, external_sources) + for entry in entries: + for record in entry["records"]: + for feature in record["features"]: + igem = feature.get("external_function_candidates", {}).get( + "igem_registry", {} + ) + for role in igem.get("roles", []): + external_role_index[str(role["accession"])].add(entry["path"]) + for category in igem.get("categories", []): + external_category_index[str(category["path"])].add(entry["path"]) + for part in igem.get("specific_parts", []): + if not part.get("selected"): + continue + external_part_index[str(part["name"])].add(entry["path"]) + part_role = part.get("role") + if isinstance(part_role, Mapping) and part_role.get("accession"): + external_role_index[str(part_role["accession"])].add( + entry["path"] + ) + no_features = [ entry["path"] for entry in entries if entry["part_feature_count"] == 0 ] @@ -208,6 +234,9 @@ def build_cloning_inventory( "files_with_no_part_features": no_features, "files_with_no_primers": no_primers, "duplicate_sequence_group_count": len(duplicate_sequences), + "features_with_external_candidates": external_enrichment[ + "matched_feature_count" + ], }, "files": entries, "indexes": { @@ -217,6 +246,9 @@ def build_cloning_inventory( "sequence_sha256_to_files": _sorted_set_mapping(sequence_index), "duplicate_sequences": duplicate_sequences, "neb_enzyme_to_files": dict(sorted(enzyme_index.items())), + "igem_role_to_files": _sorted_set_mapping(external_role_index), + "igem_category_to_files": _sorted_set_mapping(external_category_index), + "igem_part_to_files": _sorted_set_mapping(external_part_index), }, "neb_enzyme_catalog": enzyme_catalog, "errors": errors, @@ -228,10 +260,136 @@ def build_cloning_inventory( if include_enzymes else "not generated" ), + "external_enrichment": external_enrichment, }, } +def _apply_external_enrichment( + entries: Sequence[dict[str, Any]], external_sources: Mapping[str, Any] | None +) -> dict[str, Any]: + if not external_sources: + return { + "status": "not provided", + "matched_feature_count": 0, + "method": "none", + } + sources = external_sources.get("sources") + if not isinstance(sources, Mapping): + raise ValueError("External source data has no sources mapping") + igem = sources.get("igem_registry") + if not isinstance(igem, Mapping): + return { + "status": "iGEM vocabulary absent", + "matched_feature_count": 0, + "method": "none", + } + indexes = igem.get("indexes") + if not isinstance(indexes, Mapping): + raise ValueError("iGEM source data has no indexes mapping") + role_index = indexes.get("role_alias_to_terms") + category_index = indexes.get("category_alias_to_terms") + part_index = indexes.get("feature_signature_to_part_matches", {}) + if not isinstance(role_index, Mapping) or not isinstance(category_index, Mapping): + raise ValueError("iGEM source data has invalid alias indexes") + if not isinstance(part_index, Mapping): + raise ValueError("iGEM source data has invalid specific-part index") + + matched_count = 0 + specific_part_candidate_count = 0 + selected_specific_part_count = 0 + for entry in entries: + for record in entry["records"]: + for feature in record["features"]: + match_keys = [str(feature.get("normalized_label") or "")] + match_keys = [key for key in match_keys if key] + roles = _external_terms_for_aliases(role_index, match_keys, "accession") + categories = _external_terms_for_aliases( + category_index, match_keys, "uuid" + ) + signature = str(feature.get("feature_signature") or "") + raw_parts = part_index.get(signature, []) + specific_parts = [ + {str(key): value for key, value in part.items()} + for part in raw_parts + if isinstance(part, Mapping) + ] + selected_parts = [ + part for part in specific_parts if part.get("selected") + ] + if not roles and not categories and not specific_parts: + continue + feature["external_function_candidates"] = { + "igem_registry": { + "match_keys": match_keys, + "match_method": ( + "specific part search by source label, sequence length, " + "functional-role filter, and sequence/translation evidence; " + "vocabulary lookup by source label only" + ), + "nucleotide_sequence_verified": any( + bool(part.get("evidence", {}).get("nucleotide_exact")) + for part in selected_parts + ), + "translated_peptide_verified": any( + bool( + part.get("evidence", {}).get("translated_peptide_exact") + ) + for part in selected_parts + ), + "review_required": not selected_parts + or any( + bool(part.get("review_required")) for part in selected_parts + ), + "roles": roles, + "categories": categories, + "specific_parts": specific_parts, + } + } + matched_count += 1 + specific_part_candidate_count += bool(specific_parts) + selected_specific_part_count += bool(selected_parts) + + source = igem.get("source") + summary = igem.get("summary") + source_mapping = source if isinstance(source, Mapping) else {} + summary_mapping = summary if isinstance(summary, Mapping) else {} + openapi = source_mapping.get("openapi") + openapi_mapping = openapi if isinstance(openapi, Mapping) else {} + return { + "status": "applied", + "matched_feature_count": matched_count, + "specific_part_candidate_count": specific_part_candidate_count, + "specific_part_match_count": selected_specific_part_count, + "selected_specific_part_count": selected_specific_part_count, + "method": ( + "source-label vocabulary lookup plus provenance-preserving specific-part " + "matches; GenBank structural types are never treated as iGEM roles" + ), + "igem_api_version": str(source_mapping.get("api_version") or ""), + "igem_openapi_sha256": str(openapi_mapping.get("sha256") or ""), + "igem_role_count": int(summary_mapping.get("role_count") or 0), + "igem_category_count": int(summary_mapping.get("category_count") or 0), + } + + +def _external_terms_for_aliases( + index: Mapping[str, Any], aliases: Sequence[str], identity_key: str +) -> list[dict[str, Any]]: + terms: dict[str, dict[str, Any]] = {} + for alias in aliases: + candidates = index.get(alias, []) + if not isinstance(candidates, list): + continue + for candidate in candidates: + if not isinstance(candidate, Mapping): + continue + identity = str(candidate.get(identity_key) or "") + if identity: + terms[identity] = {str(key): value for key, value in candidate.items()} + return sorted(terms.values(), key=lambda term: str(term.get(identity_key) or "")) + + def _inventory_file( path: Path, display_path: str, include_enzymes: bool ) -> dict[str, Any]: @@ -279,7 +437,7 @@ def _inventory_file( def _inventory_record(record: Any, *, include_enzymes: bool) -> dict[str, Any]: sequence = str(record.seq).upper() feature_entries = [ - _inventory_feature(feature) + _inventory_feature(record, feature) for feature in record.features if feature.type.casefold() not in SOURCE_TYPES | PRIMER_TYPES ] @@ -309,10 +467,11 @@ def _inventory_record(record: Any, *, include_enzymes: bool) -> dict[str, Any]: return result -def _inventory_feature(feature: Any) -> dict[str, Any]: +def _inventory_feature(record: Any, feature: Any) -> dict[str, Any]: qualifiers = _clean_qualifiers(feature.qualifiers) label = _first_qualifier(qualifiers, LABEL_QUALIFIERS) roles = classify_feature_roles(feature.type, qualifiers) + sequence = str(feature.extract(record.seq)).upper() return { "feature_type": str(feature.type), "location": str(feature.location), @@ -321,12 +480,63 @@ def _inventory_feature(feature: Any) -> dict[str, Any]: "strand": feature.location.strand, "label": label, "normalized_label": normalize_label(label), + "feature_sequence_length": len(sequence), + "feature_sequence_sha256": hashlib.sha256(sequence.encode("ascii")).hexdigest(), + "feature_signature": feature_signature(label, sequence), "functional_roles": roles, "functional_description": describe_feature(label, roles), "qualifiers": qualifiers, } +def collect_feature_queries(inputs: Sequence[Path]) -> list[dict[str, Any]]: + """Extract deduplicated, sequence-bearing feature queries for external search.""" + from Bio import SeqIO + + queries: dict[str, dict[str, Any]] = {} + for path in discover_sequence_files(inputs): + file_format = SEQUENCE_FORMATS[path.suffix.lower()] + with path.open(encoding="utf-8") as handle: + records = SeqIO.parse(handle, file_format) # type: ignore[no-untyped-call] + for record in records: + for feature in record.features: + if feature.type.casefold() in SOURCE_TYPES | PRIMER_TYPES: + continue + qualifiers = _clean_qualifiers(feature.qualifiers) + label = _first_qualifier(qualifiers, LABEL_QUALIFIERS) + if not label: + continue + sequence = str(feature.extract(record.seq)).upper() + signature = feature_signature(label, sequence) + queries.setdefault( + signature, + { + "feature_signature": signature, + "label": label, + "normalized_label": normalize_label(label), + "feature_type": str(feature.type), + "functional_roles": classify_feature_roles( + feature.type, qualifiers + ), + "nucleotide_sequence": sequence, + "nucleotide_length": len(sequence), + "translation": _first_qualifier( + qualifiers, ("translation",) + ), + "source_description": _first_qualifier( + qualifiers, ("function", "product", "note") + ), + }, + ) + return [queries[key] for key in sorted(queries)] + + +def feature_signature(label: str, sequence: str) -> str: + """Return a stable identity for one labeled feature sequence.""" + sequence_digest = hashlib.sha256(sequence.upper().encode("ascii")).hexdigest() + return f"{normalize_label(label)}:{sequence_digest}" + + def _inventory_primer(record: Any, feature: Any) -> dict[str, Any]: qualifiers = _clean_qualifiers(feature.qualifiers) label = _first_qualifier(qualifiers, LABEL_QUALIFIERS) or "unnamed primer" @@ -462,15 +672,13 @@ def _neb_restriction_sites(record: Any) -> dict[str, int]: def _display_path(path: Path, root: Path) -> str: try: - return path.relative_to(root).as_posix() + return os.fspath(path.relative_to(root)) except ValueError: return str(path) def _common_parent(files: Sequence[Path], inputs: Sequence[Path]) -> Path: if files: - import os - return Path(os.path.commonpath([str(path.parent) for path in files])) if inputs: first = inputs[0].expanduser().resolve() diff --git a/src/lab_bench_2/plannotate_runner.py b/src/lab_bench_2/plannotate_runner.py new file mode 100644 index 0000000..2644483 --- /dev/null +++ b/src/lab_bench_2/plannotate_runner.py @@ -0,0 +1,410 @@ +"""Run pLannotate as a separate, provenance-preserving enrichment pass.""" + +from __future__ import annotations + +import csv +import hashlib +import json +import os +import re +import shutil +import subprocess +from dataclasses import dataclass +from datetime import UTC, datetime +from pathlib import Path +from typing import Any, Callable, Sequence + +from lab_bench_2.cloning_inventory import SEQUENCE_FORMATS, discover_sequence_files + +PLANNOTATE_VERSION = "2.0.0" +DEFAULT_CACHE_DIR = Path.home() / ".cache/labbench2/plannotate-v2" +EXACT_IDENTITY_PERCENT = 100.0 +FULL_LENGTH_COVERAGE_PERCENT = 99.9 + + +class PlannotateError(RuntimeError): + """Raised when the managed pLannotate enrichment cannot run safely.""" + + +@dataclass(frozen=True) +class PlannotateEnvironment: + """Resolved executable and database provenance for one annotation run.""" + + command: list[str] + version_report: str + database_manifest: dict[str, Any] + managed_environment: bool + + +RunCommand = Callable[..., subprocess.CompletedProcess[str]] +WhichCommand = Callable[[str], str | None] + + +def prepare_plannotate( + cache_dir: Path = DEFAULT_CACHE_DIR, + *, + auto_setup: bool = True, + run: RunCommand = subprocess.run, + which: WhichCommand = shutil.which, +) -> PlannotateEnvironment: + """Find pLannotate or create a cached conda environment, then install its DBs.""" + cache_dir = cache_dir.expanduser().resolve() + executable = which("plannotate") + managed = False + if executable: + command = [executable] + else: + managed = True + command = [_managed_executable(cache_dir)] + if not Path(command[0]).is_file(): + if not auto_setup: + raise PlannotateError( + "pLannotate is not installed and automatic setup is disabled" + ) + _create_managed_environment(cache_dir, run=run, which=which) + if not Path(command[0]).is_file(): + raise PlannotateError("Managed pLannotate setup did not create its CLI") + + database_result = _run([*command, "databases"], run=run, timeout=60, check=False) + if database_result.returncode != 0: + if not auto_setup: + raise PlannotateError( + "pLannotate databases are absent and automatic setup is disabled" + ) + _run([*command, "setupdb"], run=run, timeout=900) + database_result = _run([*command, "databases"], run=run, timeout=60) + try: + database_manifest = json.loads(database_result.stdout) + except json.JSONDecodeError as exc: + raise PlannotateError( + "pLannotate returned an invalid database manifest" + ) from exc + if not isinstance(database_manifest, dict): + raise PlannotateError("pLannotate database manifest is not an object") + version_result = _run([*command, "--version"], run=run, timeout=60) + version_report = version_result.stdout.strip() + if f"pLannotate {PLANNOTATE_VERSION}" not in version_report: + raise PlannotateError( + f"Expected pLannotate {PLANNOTATE_VERSION}, got: {version_report}" + ) + return PlannotateEnvironment( + command=command, + version_report=version_report, + database_manifest=database_manifest, + managed_environment=managed, + ) + + +def annotate_with_plannotate( + inputs: Sequence[Path], + output_dir: Path, + *, + cache_dir: Path = DEFAULT_CACHE_DIR, + fast: bool = True, + cores: int = 4, + refresh: bool = False, + auto_setup: bool = True, + run: RunCommand = subprocess.run, + which: WhichCommand = shutil.which, +) -> dict[str, Any]: + """Annotate files into a separate directory and return a source-linked manifest.""" + if cores < 1: + raise ValueError("cores must be positive") + files = discover_sequence_files(inputs) + if not files: + raise ValueError("No supported sequence files found") + output_dir = output_dir.expanduser().resolve() + output_dir.mkdir(parents=True, exist_ok=True) + environment = prepare_plannotate( + cache_dir, auto_setup=auto_setup, run=run, which=which + ) + prior = _read_prior_manifest(output_dir / "plannotate-manifest.json") + prior_by_hash = { + str(item.get("source_sha256")): item + for item in prior.get("results", []) + if isinstance(item, dict) + } + + results: list[dict[str, Any]] = [] + errors: list[dict[str, str]] = [] + for source_file in files: + source_hash = _sha256_file(source_file) + cached = prior_by_hash.get(source_hash) + if not refresh and cached and _outputs_exist(output_dir, cached): + results.append({**cached, "status": "cache-hit"}) + continue + output_stem = _output_stem(source_file, source_hash) + command = build_annotation_command( + environment.command, + source_file, + output_dir, + output_stem, + fast=fast, + cores=cores, + ) + try: + _run(command, run=run, timeout=1800) + except PlannotateError as exc: + errors.append({"source_path": str(source_file), "error": str(exc)}) + continue + outputs = sorted( + path.name for path in output_dir.glob(f"{output_stem}*") if path.is_file() + ) + if not outputs: + errors.append( + { + "source_path": str(source_file), + "error": "pLannotate created no output files", + } + ) + continue + results.append( + { + "source_path": str(source_file), + "source_sha256": source_hash, + "source_format": SEQUENCE_FORMATS[source_file.suffix.casefold()], + "status": "annotated", + "mode": "fast" if fast else "full", + "linear": _is_linear(source_file), + "outputs": outputs, + "output_sha256": { + name: _sha256_file(output_dir / name) for name in outputs + }, + "annotation_summary": _summarize_plannotate_outputs( + output_dir, outputs + ), + } + ) + + manifest = { + "schema_version": 1, + "generated_at": datetime.now(UTC).isoformat(), + "tool": { + "name": "pLannotate", + "version": PLANNOTATE_VERSION, + "version_report": environment.version_report, + "managed_environment": environment.managed_environment, + "database_manifest": environment.database_manifest, + }, + "policy": { + "source_files_modified": False, + "outputs_are_enrichment_candidates": True, + "mode": "fast" if fast else "full", + }, + "summary": { + "discovered_file_count": len(files), + "annotated_or_cached_count": len(results), + "error_count": len(errors), + }, + "results": results, + "errors": errors, + } + _write_json_atomic(output_dir / "plannotate-manifest.json", manifest) + return manifest + + +def build_annotation_command( + plannotate_command: Sequence[str], + source_file: Path, + output_dir: Path, + output_stem: str, + *, + fast: bool, + cores: int, +) -> list[str]: + """Build the pLannotate CLI invocation for one source file.""" + command = [ + *plannotate_command, + "batch", + "--input", + str(source_file), + "--output", + str(output_dir), + "--file-name", + output_stem, + "--suffix", + "", + "--csv", + "--cores", + str(cores), + ] + if _is_linear(source_file): + command.append("--linear") + if fast: + command.append("--fast") + return command + + +def _create_managed_environment( + cache_dir: Path, *, run: RunCommand, which: WhichCommand +) -> None: + manager = next( + ( + resolved + for name in ("micromamba", "mamba", "conda") + if (resolved := which(name)) + ), + None, + ) + if manager is None: + raise PlannotateError( + "Automatic pLannotate setup requires micromamba, mamba, or conda" + ) + cache_dir.parent.mkdir(parents=True, exist_ok=True) + _run( + [ + manager, + "create", + "--yes", + "--prefix", + str(cache_dir), + "--channel", + "conda-forge", + "--channel", + "bioconda", + f"plannotate={PLANNOTATE_VERSION}", + ], + run=run, + timeout=1800, + ) + + +def _managed_executable(cache_dir: Path) -> str: + unix = cache_dir / "bin/plannotate" + windows = cache_dir / "Scripts/plannotate.exe" + return str(windows if windows.is_file() else unix) + + +def _is_linear(source_file: Path) -> bool: + if SEQUENCE_FORMATS[source_file.suffix.casefold()] == "fasta": + return True + try: + from Bio import SeqIO + + with source_file.open(encoding="utf-8") as handle: + records = list(SeqIO.parse(handle, "genbank")) # type: ignore[no-untyped-call] + except (ModuleNotFoundError, OSError, ValueError): + return False + return bool(records) and all( + str(record.annotations.get("topology") or "").casefold() == "linear" + for record in records + ) + + +def _output_stem(source_file: Path, source_hash: str) -> str: + safe_name = re.sub(r"[^A-Za-z0-9._-]+", "_", source_file.stem).strip("_") + return f"{safe_name or 'sequence'}-{source_hash[:12]}" + + +def _read_prior_manifest(path: Path) -> dict[str, Any]: + try: + payload = json.loads(path.read_text(encoding="utf-8")) + return payload if isinstance(payload, dict) else {} + except (OSError, json.JSONDecodeError): + return {} + + +def _outputs_exist(output_dir: Path, item: dict[str, Any]) -> bool: + outputs = item.get("outputs") + hashes = item.get("output_sha256") + if not isinstance(outputs, list) or not isinstance(hashes, dict) or not outputs: + return False + return all( + isinstance(name, str) + and (output_dir / name).is_file() + and _sha256_file(output_dir / name) == hashes.get(name) + for name in outputs + ) + + +def _summarize_plannotate_outputs( + output_dir: Path, outputs: Sequence[str] +) -> dict[str, Any]: + annotations: list[dict[str, Any]] = [] + database_counts: dict[str, int] = {} + for name in outputs: + if Path(name).suffix.casefold() != ".csv": + continue + with (output_dir / name).open(encoding="utf-8", newline="") as handle: + for row in csv.DictReader(handle): + database = str(row.get("database") or "") + fragment = str(row.get("fragment") or "").casefold() == "true" + identity = _optional_float(row.get("percent identity")) + match_coverage = _optional_float(row.get("percent match length")) + annotations.append( + { + "feature": str(row.get("Feature") or ""), + "feature_type": str(row.get("Type") or ""), + "description": str(row.get("Description") or ""), + "database": database, + "fragment": fragment, + "percent_identity": identity, + "percent_match_length": match_coverage, + } + ) + database_counts[database] = database_counts.get(database, 0) + 1 + return { + "annotation_count": len(annotations), + "fragment_count": sum(item["fragment"] for item in annotations), + "non_fragment_count": sum(not item["fragment"] for item in annotations), + "exact_full_length_count": sum( + item["percent_identity"] == EXACT_IDENTITY_PERCENT + and item["percent_match_length"] is not None + and item["percent_match_length"] >= FULL_LENGTH_COVERAGE_PERCENT + for item in annotations + ), + "database_counts": dict(sorted(database_counts.items())), + "annotations": annotations, + } + + +def _optional_float(value: object) -> float | None: + try: + return float(str(value)) + except (TypeError, ValueError): + return None + + +def _sha256_file(path: Path) -> str: + digest = hashlib.sha256() + with path.open("rb") as handle: + for chunk in iter(lambda: handle.read(1024 * 1024), b""): + digest.update(chunk) + return digest.hexdigest() + + +def _write_json_atomic(path: Path, payload: dict[str, Any]) -> None: + temporary = path.with_suffix(path.suffix + ".tmp") + temporary.write_text( + json.dumps(payload, indent=2, sort_keys=True) + "\n", encoding="utf-8" + ) + temporary.replace(path) + + +def _run( + command: Sequence[str], + *, + run: RunCommand, + timeout: float, + check: bool = True, +) -> subprocess.CompletedProcess[str]: + executable_directory = str(Path(command[0]).expanduser().resolve().parent) + environment = { + **os.environ, + "PATH": executable_directory + os.pathsep + os.environ.get("PATH", ""), + } + try: + result = run( + list(command), + capture_output=True, + text=True, + timeout=timeout, + check=False, + env=environment, + ) + except (OSError, subprocess.TimeoutExpired) as exc: + raise PlannotateError(f"Command failed to start: {command[0]}: {exc}") from exc + if check and result.returncode != 0: + detail = (result.stderr or result.stdout or "unknown error").strip() + raise PlannotateError(f"pLannotate failed: {detail}") + return result diff --git a/tests/lab_bench_2/test_cloning_external_sources.py b/tests/lab_bench_2/test_cloning_external_sources.py new file mode 100644 index 0000000..3ce2702 --- /dev/null +++ b/tests/lab_bench_2/test_cloning_external_sources.py @@ -0,0 +1,361 @@ +from __future__ import annotations + +import json +from pathlib import Path +from typing import Any +from urllib.parse import parse_qs, urlsplit + +import pytest + +from lab_bench_2.cloning_external_sources import ( + CachedHttpClient, + CachedResponse, + ExternalSourceError, + IgemRegistryClient, + fetch_rebase_neb_catalog, + parse_rebase_enzymes, + parse_rebase_references, +) +from lab_bench_2.cloning_inventory import feature_signature + + +class NetworkResponse: + def __init__(self, url: str) -> None: + self.url = url + self.headers = {"Content-Type": "application/json"} + + def __enter__(self) -> NetworkResponse: + return self + + def __exit__(self, *args: object) -> None: + return None + + def read(self) -> bytes: + return b'{"ok": true}' + + def geturl(self) -> str: + return self.url + + +def _response(body: str, url: str = "https://example.invalid/source") -> CachedResponse: + return CachedResponse( + body=body.encode(), + url=url, + fetched_at="2026-08-27T00:00:00+00:00", + sha256="digest", + headers={}, + cache_status="fetched", + ) + + +class RebaseHttp: + def __init__(self, payloads: dict[str, str]) -> None: + self.payloads = payloads + + def get(self, url: str, *, cache_key: str, refresh: bool) -> CachedResponse: + del cache_key, refresh + return _response(self.payloads[url.rsplit("/", 1)[-1]], url) + + +class IgemHttp: + def get_json( + self, url: str, *, cache_key: str, refresh: bool + ) -> tuple[Any, CachedResponse]: + del cache_key, refresh + parsed = urlsplit(url) + page = int(parse_qs(parsed.query).get("page", ["1"])[0]) + payload: Any + if parsed.path == "/docs-json": + payload = {"info": {"title": "Registry", "version": "1.0"}} + elif parsed.path == "/v1/parts/roles": + payload = { + "data": ( + [ + { + "accession": "SO:0000316", + "label": "CDS", + "definition": "protein coding region", + "source": "SO", + "synonyms": ["coding sequence"], + }, + { + "uuid": "role-tag", + "accession": "SO:0000324", + "label": "Tag", + "definition": "identifying sequence tag", + "source": "SO", + "synonyms": [], + }, + ] + if page == 1 + else [] + ), + "page": page, + "total": 2, + } + elif parsed.path == "/v1/parts/categories": + payload = { + "data": ( + [ + { + "uuid": "category-cas9", + "label": "//function/crispr/cas9", + "value": "//function/crispr/cas9", + "description": "", + } + ] + if page == 1 + else [] + ), + "page": page, + "total": 1, + } + elif parsed.path == "/v1/parts": + payload = {"data": [], "page": 1, "total": 42} + elif parsed.path == "/v1/parts/slugs/bba-k4587111": + payload = { + "uuid": "part-flag", + "name": "BBa_K4587111", + "slug": "bba-k4587111", + "title": "FLAG Tag", + "description": "A detection tag.", + "sequence": "GACTACAAGGACGATGATGATAAA", + "usageCount": 6, + "role": { + "uuid": "role-tag", + "accession": "SO:0000324", + "label": "Tag", + }, + } + else: + raise AssertionError(url) + return payload, _response(json.dumps(payload), url) + + +class IgemWrongRoleHttp(IgemHttp): + def get_json( + self, url: str, *, cache_key: str, refresh: bool + ) -> tuple[Any, CachedResponse]: + parsed = urlsplit(url) + query = parse_qs(parsed.query) + if parsed.path == "/v1/parts" and query.get("search"): + payload = { + "data": [ + { + "name": "BBa_TEST_WRONG_ROLE", + "slug": "bba-test-wrong-role", + "title": "Test promoter", + "sequenceLength": 6, + "usageCount": 1, + "role": {"accession": "SO:0000316", "label": "CDS"}, + } + ], + "page": 1, + "total": 1, + } + return payload, _response(json.dumps(payload), url) + if parsed.path == "/v1/parts/slugs/bba-test-wrong-role": + payload = { + "uuid": "part-wrong-role", + "name": "BBa_TEST_WRONG_ROLE", + "slug": "bba-test-wrong-role", + "title": "Test promoter", + "description": "An intentionally wrong-role test candidate.", + "sequence": "AAAAAA", + "usageCount": 1, + "role": {"accession": "SO:0000316", "label": "CDS"}, + } + return payload, _response(json.dumps(payload), url) + return super().get_json(url, cache_key=cache_key, refresh=refresh) + + +def test_http_cache_reuses_validated_response_and_rejects_other_hosts( + tmp_path: Path, +) -> None: + calls: list[str] = [] + + def open_url(request: Any, *, timeout: float) -> NetworkResponse: + del timeout + calls.append(str(request.full_url)) + return NetworkResponse(str(request.full_url)) + + client = CachedHttpClient( + tmp_path, + allowed_hosts=["api.registry.igem.org"], + min_delay=0, + open_url=open_url, + ) + + first = client.get("https://api.registry.igem.org/test", cache_key="test-response") + second = client.get("https://api.registry.igem.org/test", cache_key="test-response") + third = client.get( + "https://api.registry.igem.org/test?different=1", + cache_key="test-response", + ) + + assert first.cache_status == "fetched" + assert second.cache_status == "cache-hit" + assert third.cache_status == "fetched" + assert calls == [ + "https://api.registry.igem.org/test", + "https://api.registry.igem.org/test?different=1", + ] + with pytest.raises(ExternalSourceError, match="non-allowlisted"): + client.get("https://example.org/test", cache_key="bad-host") + + +def test_parses_rebase_cut_geometry_and_supplier_records() -> None: + enzyme_text = """# REBASE version 608 +AatII GACGTC 6 2 0 5 1 0 0 +""" + reference_text = """# REBASE version 608 +AatII +Acetobacter aceti + +2(6) +IFO 3281 +BINRV +1 +Example reference. +// +""" + + assert parse_rebase_enzymes(enzyme_text)["AatII"] == { + "name": "AatII", + "recognition_site": "GACGTC", + "recognition_length": 6, + "cut_count": 2, + "blunt": False, + "cut_5prime": 5, + "cut_3prime": 1, + "second_cut_5prime": 0, + "second_cut_3prime": 0, + } + assert parse_rebase_references(reference_text)["AatII"]["supplier_codes"] == [ + "B", + "I", + "N", + "R", + "V", + ] + + +def test_rebase_catalog_filters_to_neb_and_records_release(tmp_path: Path) -> None: + common = "# REBASE version 608\n# Rich Roberts Jul 31 2026\n" + client = RebaseHttp( + { + "link_emboss_e": common + + "AatII GACGTC 6 2 0 5 1 0 0\nEcoRV GATATC 6 2 1 3 3 0 0\n", + "link_emboss_s": common + + "B Thermo Fisher Scientific\nN New England Biolabs\n", + "link_emboss_r": common + + "AatII\norganism\n\n\nsource\nN\n0\n//\n" + + "EcoRV\norganism\n\n\nsource\nB\n0\n//\n", + } + ) + + catalog = fetch_rebase_neb_catalog(tmp_path, http=client) # type: ignore[arg-type] + + assert catalog["source"]["release"] == "608" + assert catalog["source"]["release_date"] == "Jul 31 2026" + assert [item["name"] for item in catalog["enzymes"]] == ["AatII"] + + +def test_rebase_rejects_mismatched_releases(tmp_path: Path) -> None: + client = RebaseHttp( + { + "link_emboss_e": "# REBASE version 608\nAatII GACGTC 6 2 0 5 1 0 0\n", + "link_emboss_s": "# REBASE version 607\nN New England Biolabs\n", + "link_emboss_r": "# REBASE version 608\n", + } + ) + + with pytest.raises(ExternalSourceError, match="disagree"): + fetch_rebase_neb_catalog(tmp_path, http=client) # type: ignore[arg-type] + + +def test_igem_vocabulary_builds_role_and_function_aliases(tmp_path: Path) -> None: + vocabulary = IgemRegistryClient( + tmp_path, + http=IgemHttp(), # type: ignore[arg-type] + ).fetch_vocabulary() + + assert vocabulary["summary"] == { + "published_part_count": 42, + "role_count": 2, + "category_count": 1, + "functional_category_count": 1, + "feature_query_count": 0, + "specific_part_match_count": 0, + "part_match_error_count": 0, + } + assert ( + vocabulary["indexes"]["role_alias_to_terms"]["coding sequence"][0]["accession"] + == "SO:0000316" + ) + assert ( + vocabulary["indexes"]["category_alias_to_terms"]["cas9"][0]["path"] + == "//function/crispr/cas9" + ) + + +def test_igem_flag_match_is_specific_and_translation_verified(tmp_path: Path) -> None: + source_sequence = "GATTACAAAGACGATGACGATAAG" + signature = feature_signature("FLAG", source_sequence) + vocabulary = IgemRegistryClient( + tmp_path, + http=IgemHttp(), # type: ignore[arg-type] + ).fetch_vocabulary( + [ + { + "feature_signature": signature, + "label": "FLAG", + "normalized_label": "flag", + "functional_roles": ["coding_sequence", "epitope_tag"], + "nucleotide_sequence": source_sequence, + "nucleotide_length": len(source_sequence), + "translation": "DYKDDDDK", + } + ] + ) + + match = vocabulary["indexes"]["feature_signature_to_part_matches"][signature][0] + assert match["name"] == "BBa_K4587111" + assert match["url"] == "https://registry.igem.org/parts/bba-k4587111" + assert match["description"] == "A detection tag." + assert match["role"]["accession"] == "SO:0000324" + assert match["role"]["label"] == "Tag" + assert match["evidence"]["nucleotide_exact"] is False + assert match["evidence"]["translated_peptide_exact"] is True + assert match["evidence"]["expected_functional_role"] is True + assert match["evidence"]["source_feature_type_used_as_role"] is False + assert match["selected"] is True + assert match["review_required"] is False + + +def test_igem_does_not_select_exact_sequence_with_wrong_role(tmp_path: Path) -> None: + source_sequence = "AAAAAA" + signature = feature_signature("Test promoter", source_sequence) + vocabulary = IgemRegistryClient( + tmp_path, + http=IgemWrongRoleHttp(), # type: ignore[arg-type] + ).fetch_vocabulary( + [ + { + "feature_signature": signature, + "label": "Test promoter", + "normalized_label": "test promoter", + "functional_roles": ["promoter"], + "nucleotide_sequence": source_sequence, + "nucleotide_length": len(source_sequence), + "translation": "", + } + ] + ) + + match = vocabulary["indexes"]["feature_signature_to_part_matches"][signature][0] + assert match["evidence"]["nucleotide_exact"] is True + assert match["evidence"]["expected_functional_role"] is False + assert match["evidence"]["source_feature_type_used_as_role"] is False + assert match["selected"] is False + assert vocabulary["summary"]["specific_part_match_count"] == 0 diff --git a/tests/lab_bench_2/test_cloning_inventory.py b/tests/lab_bench_2/test_cloning_inventory.py index 6e4c77a..7df61de 100644 --- a/tests/lab_bench_2/test_cloning_inventory.py +++ b/tests/lab_bench_2/test_cloning_inventory.py @@ -8,6 +8,7 @@ from lab_bench_2.cloning_inventory import ( build_cloning_inventory, classify_feature_roles, + feature_signature, ) @@ -107,3 +108,89 @@ def test_inventory_includes_neb_catalog_and_cut_counts(tmp_path: Path) -> None: def test_feature_type_and_function_are_separate() -> None: roles = classify_feature_roles("CDS", {"label": ["FLAG-NLS"]}) assert roles == ["coding_sequence", "epitope_tag", "localization_signal"] + + +def test_external_igem_specific_match_preserves_evidence_and_indexes( + tmp_path: Path, +) -> None: + record = SeqRecord(Seq("ATG" * 30), id="record", name="record") + record.annotations["molecule_type"] = "DNA" + record.features = [ + SeqFeature( + FeatureLocation(0, 30), + type="CDS", + qualifiers={"label": ["Cas9"]}, + ) + ] + gbk = tmp_path / "cas9.gbk" + SeqIO.write(record, gbk, "genbank") + signature = feature_signature("Cas9", "ATG" * 10) + source_data = { + "sources": { + "igem_registry": { + "source": { + "api_version": "1.0", + "openapi": {"sha256": "openapi-digest"}, + }, + "summary": {"role_count": 1, "category_count": 1}, + "indexes": { + "role_alias_to_terms": { + "cds": [ + { + "accession": "SO:0000316", + "label": "CDS", + "definition": "coding region", + "source": "SO", + } + ] + }, + "category_alias_to_terms": { + "cas9": [ + { + "uuid": "category-cas9", + "path": "//function/crispr/cas9", + "description": "", + } + ] + }, + "feature_signature_to_part_matches": { + signature: [ + { + "name": "BBa_TEST_CAS9", + "url": "https://registry.igem.org/parts/bba-test-cas9", + "title": "Cas9", + "description": "RNA-guided nuclease.", + "role": { + "accession": "SO:0000316", + "label": "CDS", + }, + "selected": True, + "evidence": {"nucleotide_exact": True}, + } + ] + }, + }, + } + } + } + + inventory = build_cloning_inventory( + [gbk], include_enzymes=False, external_sources=source_data + ) + candidates = inventory["files"][0]["records"][0]["features"][0][ + "external_function_candidates" + ]["igem_registry"] + + assert candidates["nucleotide_sequence_verified"] is True + assert candidates["translated_peptide_verified"] is False + assert candidates["review_required"] is False + assert candidates["roles"] == [] + assert candidates["categories"][0]["path"] == "//function/crispr/cas9" + assert candidates["specific_parts"][0]["name"] == "BBa_TEST_CAS9" + assert inventory["indexes"]["igem_category_to_files"] == { + "//function/crispr/cas9": ["cas9.gbk"] + } + assert inventory["indexes"]["igem_part_to_files"] == {"BBa_TEST_CAS9": ["cas9.gbk"]} + enrichment = inventory["provenance"]["external_enrichment"] + assert enrichment["specific_part_candidate_count"] == 1 + assert enrichment["selected_specific_part_count"] == 1 diff --git a/tests/lab_bench_2/test_plannotate_runner.py b/tests/lab_bench_2/test_plannotate_runner.py new file mode 100644 index 0000000..fa9af3e --- /dev/null +++ b/tests/lab_bench_2/test_plannotate_runner.py @@ -0,0 +1,118 @@ +from __future__ import annotations + +import json +import subprocess +from pathlib import Path + +from lab_bench_2.plannotate_runner import ( + annotate_with_plannotate, + build_annotation_command, + prepare_plannotate, +) + + +class FakePlannotate: + def __init__(self) -> None: + self.commands: list[list[str]] = [] + self.environments: list[dict[str, str]] = [] + + def __call__( + self, command: list[str], **kwargs: object + ) -> subprocess.CompletedProcess[str]: + environment = kwargs.get("env") + assert isinstance(environment, dict) + self.environments.append(environment) + self.commands.append(command) + if command[-1] == "databases": + return subprocess.CompletedProcess( + command, + 0, + stdout=json.dumps( + { + "schema_version": 1, + "bundle": "plannotate-databases-v2", + "files": {"BLAST_dbs/snapgene.db": {"sha256": "db"}}, + } + ), + stderr="", + ) + if command[-1] == "--version": + return subprocess.CompletedProcess( + command, + 0, + stdout="pLannotate 2.0.0\ndatabase: plannotate-databases-v2\n", + stderr="", + ) + if "batch" in command: + output_dir = Path(command[command.index("--output") + 1]) + stem = command[command.index("--file-name") + 1] + (output_dir / f"{stem}.gbk").write_text("annotated genbank") + (output_dir / f"{stem}.csv").write_text( + "Feature,Type,Description,database,fragment,percent identity,percent match length\n" + "Cas9,CDS,editor,snapgene,False,100.0,100.0\n" + ) + return subprocess.CompletedProcess(command, 0, stdout="", stderr="") + raise AssertionError(command) + + +def test_prepare_reads_version_and_database_provenance(tmp_path: Path) -> None: + fake = FakePlannotate() + + environment = prepare_plannotate( + tmp_path, + run=fake, + which=lambda name: "/usr/local/bin/plannotate" + if name == "plannotate" + else None, + ) + + assert environment.command == ["/usr/local/bin/plannotate"] + assert environment.database_manifest["bundle"] == "plannotate-databases-v2" + assert environment.managed_environment is False + assert fake.environments[0]["PATH"].split(":", 1)[0] == "/usr/local/bin" + + +def test_fasta_annotation_command_is_linear_and_fast(tmp_path: Path) -> None: + fasta = tmp_path / "sequence.fa" + fasta.write_text(">sequence\nACGT\n") + + command = build_annotation_command( + ["plannotate"], fasta, tmp_path / "out", "sequence-hash", fast=True, cores=3 + ) + + assert command[-2:] == ["--linear", "--fast"] + assert command[command.index("--cores") + 1] == "3" + + +def test_annotations_are_separate_provenance_linked_and_cached(tmp_path: Path) -> None: + fasta = tmp_path / "sequence.fa" + fasta.write_text(">sequence\nACGTACGT\n") + output = tmp_path / "annotations" + fake = FakePlannotate() + + def which(name: str) -> str | None: + return "/usr/local/bin/plannotate" if name == "plannotate" else None + + first = annotate_with_plannotate( + [fasta], output, cache_dir=tmp_path / "cache", run=fake, which=which + ) + second = annotate_with_plannotate( + [fasta], output, cache_dir=tmp_path / "cache", run=fake, which=which + ) + refreshed = annotate_with_plannotate( + [fasta], + output, + cache_dir=tmp_path / "cache", + refresh=True, + run=fake, + which=which, + ) + + assert fasta.read_text() == ">sequence\nACGTACGT\n" + assert first["policy"]["source_files_modified"] is False + assert first["policy"]["outputs_are_enrichment_candidates"] is True + assert first["results"][0]["status"] == "annotated" + assert first["results"][0]["annotation_summary"]["exact_full_length_count"] == 1 + assert second["results"][0]["status"] == "cache-hit" + assert refreshed["results"][0]["status"] == "annotated" + assert sum("batch" in command for command in fake.commands) == 2 diff --git a/tools/annotate_cloning_sequences.py b/tools/annotate_cloning_sequences.py new file mode 100644 index 0000000..e02abff --- /dev/null +++ b/tools/annotate_cloning_sequences.py @@ -0,0 +1,67 @@ +#!/usr/bin/env python3 +"""Annotate cloning files with a managed, provenance-pinned pLannotate setup.""" + +from __future__ import annotations + +import argparse +import sys +from pathlib import Path + +from lab_bench_2.plannotate_runner import ( + DEFAULT_CACHE_DIR, + PlannotateError, + annotate_with_plannotate, +) + + +def _arguments() -> argparse.Namespace: + parser = argparse.ArgumentParser( + description=( + "Run pLannotate into a separate output directory. Missing conda " + "environment and databases are installed automatically and cached." + ) + ) + parser.add_argument("inputs", nargs="+", type=Path) + parser.add_argument("--output-dir", type=Path, required=True) + parser.add_argument("--cache-dir", type=Path, default=DEFAULT_CACHE_DIR) + parser.add_argument("--cores", type=int, default=4) + parser.add_argument( + "--full", + action="store_true", + help="Also search Swiss-Prot and Rfam; default fast mode uses SnapGene/FPbase", + ) + parser.add_argument("--refresh", action="store_true") + parser.add_argument( + "--no-auto-setup", + action="store_true", + help="Fail instead of creating the cached conda environment/databases", + ) + return parser.parse_args() + + +def main() -> int: + args = _arguments() + try: + manifest = annotate_with_plannotate( + args.inputs, + args.output_dir, + cache_dir=args.cache_dir, + fast=not args.full, + cores=args.cores, + refresh=args.refresh, + auto_setup=not args.no_auto_setup, + ) + except (OSError, PlannotateError, ValueError) as exc: + print(f"error: {exc}", file=sys.stderr) + return 2 + summary = manifest["summary"] + print( + f"pLannotate enriched {summary['annotated_or_cached_count']}/" + f"{summary['discovered_file_count']} file(s); " + f"{summary['error_count']} error(s). Output: {args.output_dir}" + ) + return 1 if summary["error_count"] else 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/tools/fetch_cloning_source_data.py b/tools/fetch_cloning_source_data.py new file mode 100644 index 0000000..1ca70ca --- /dev/null +++ b/tools/fetch_cloning_source_data.py @@ -0,0 +1,93 @@ +#!/usr/bin/env python3 +"""Fetch public iGEM and REBASE data for cloning-sequence enrichment.""" + +from __future__ import annotations + +import argparse +import json +import sys +from pathlib import Path + +from lab_bench_2.cloning_external_sources import ( + ExternalSourceError, + fetch_cloning_external_sources, +) +from lab_bench_2.cloning_inventory import collect_feature_queries + + +def _arguments() -> argparse.Namespace: + parser = argparse.ArgumentParser( + description=( + "Fetch versioned public feature vocabularies from iGEM and the " + "current supplier-N restriction catalog from REBASE." + ) + ) + parser.add_argument("--output", type=Path, required=True) + parser.add_argument( + "inputs", + nargs="*", + type=Path, + help=( + "Optional GenBank/FASTA files or directories. GenBank features are " + "used for specific, evidence-ranked iGEM part searches." + ), + ) + parser.add_argument( + "--cache-dir", + type=Path, + default=Path.home() / ".cache/labbench2/cloning-sources", + ) + parser.add_argument("--refresh", action="store_true") + parser.add_argument("--no-igem", action="store_true") + parser.add_argument("--no-rebase", action="store_true") + return parser.parse_args() + + +def main() -> int: + args = _arguments() + if args.no_igem and args.no_rebase: + print("error: at least one source must be enabled", file=sys.stderr) + return 2 + try: + feature_queries = collect_feature_queries(args.inputs) if args.inputs else [] + data = fetch_cloning_external_sources( + args.cache_dir, + refresh=args.refresh, + include_rebase=not args.no_rebase, + include_igem=not args.no_igem, + feature_queries=feature_queries, + ) + args.output.parent.mkdir(parents=True, exist_ok=True) + args.output.write_text( + json.dumps(data, indent=2, sort_keys=True) + "\n", encoding="utf-8" + ) + except (ExternalSourceError, OSError, ValueError) as exc: + print(f"error: {exc}", file=sys.stderr) + return 2 + + sources = data["sources"] + rebase = sources.get("rebase_neb") + if rebase: + print( + f"REBASE {rebase['source']['release']}: " + f"{rebase['enzyme_count']} NEB-supplied restriction enzymes" + ) + igem = sources.get("igem_registry") + if igem: + summary = igem["summary"] + print( + f"iGEM: {summary['published_part_count']} published parts, " + f"{summary['role_count']} roles, {summary['category_count']} categories" + ) + if feature_queries: + print( + f"Specific iGEM part matches: " + f"{summary['specific_part_match_count']}/" + f"{summary['feature_query_count']} feature queries" + ) + print(f"Output: {args.output}") + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/tools/inventory_cloning_sequences.py b/tools/inventory_cloning_sequences.py index 93bdb50..2eb3aef 100644 --- a/tools/inventory_cloning_sequences.py +++ b/tools/inventory_cloning_sequences.py @@ -28,23 +28,38 @@ def _arguments() -> argparse.Namespace: action="store_true", help="Skip the Biopython supplier-N enzyme catalog and site index", ) + parser.add_argument( + "--external-sources", + type=Path, + help=( + "JSON created by fetch_cloning_source_data.py; adds provenance-marked " + "specific iGEM parts, descriptions, roles, and label-vocabulary " + "candidates without replacing source annotations" + ), + ) return parser.parse_args() def main() -> int: args = _arguments() try: + external_sources = ( + json.loads(args.external_sources.read_text(encoding="utf-8")) + if args.external_sources + else None + ) inventory = build_cloning_inventory( args.inputs, root=args.root, include_enzymes=not args.no_enzymes, + external_sources=external_sources, ) args.output.parent.mkdir(parents=True, exist_ok=True) args.output.write_text( json.dumps(inventory, indent=2, sort_keys=True) + "\n", encoding="utf-8", ) - except (OSError, ValueError) as exc: + except (json.JSONDecodeError, OSError, ValueError) as exc: print(f"error: {exc}", file=sys.stderr) return 2 From bb5771a413adb75cc94f40968cbfcc6ae6c9f976 Mon Sep 17 00:00:00 2001 From: Divide-By-0 Date: Fri, 28 Aug 2026 01:01:29 -0700 Subject: [PATCH 3/3] Add one-command plasmid evidence reports --- src/lab_bench_2/CLONING_SEQUENCE_DATA.md | 56 ++- src/lab_bench_2/cloning_inventory.py | 63 ++- src/lab_bench_2/cloning_report.py | 506 ++++++++++++++++++++ tests/lab_bench_2/test_cloning_inventory.py | 17 + tests/lab_bench_2/test_cloning_report.py | 164 +++++++ tools/report_cloning_sequence.py | 213 ++++++++ 6 files changed, 999 insertions(+), 20 deletions(-) create mode 100644 src/lab_bench_2/cloning_report.py create mode 100644 tests/lab_bench_2/test_cloning_report.py create mode 100644 tools/report_cloning_sequence.py diff --git a/src/lab_bench_2/CLONING_SEQUENCE_DATA.md b/src/lab_bench_2/CLONING_SEQUENCE_DATA.md index fd58ffb..4a5e691 100644 --- a/src/lab_bench_2/CLONING_SEQUENCE_DATA.md +++ b/src/lab_bench_2/CLONING_SEQUENCE_DATA.md @@ -27,14 +27,18 @@ Those local artifacts are not additional official inputs. Keep structural annotations, inferred function, and external matches separate: -| Layer | Data | Trust and use | -| --- | --- | --- | -| Source GBK/FASTA | Sequence, coordinates, feature type, qualifiers, annotated primers | Ground truth for what the benchmark supplied | -| Biopython | Parsed records, feature/primer indexes, sequence hashes, restriction sites | Deterministic computation from the source file | -| REBASE | Current restriction-enzyme records with supplier codes | Live catalog evidence; supplier `N` is New England Biolabs | -| iGEM Registry | Specific published parts, descriptions, roles, sequences, and vocabulary | External cross-reference with explicit match evidence | -| pLannotate | SnapGene, FPbase, Swiss-Prot, and Rfam sequence-search hits | Sequence-derived enrichment candidate, kept in separate outputs | -| LLM | Prompt-ready functional description | Inference only; retain model, prompt, evidence, and confidence | +- **Source GBK/FASTA:** Sequence, coordinates, feature type, qualifiers, and + annotated primers. This is ground truth for what the benchmark supplied. +- **Biopython:** Parsed records, feature/primer indexes, sequence hashes, and + restriction sites. These are deterministic computations from the source. +- **REBASE:** Current restriction-enzyme records with supplier codes. This is + live catalog evidence; supplier `N` is New England Biolabs. +- **iGEM Registry:** Specific published parts, descriptions, roles, sequences, + and vocabulary. These are external cross-references with explicit evidence. +- **pLannotate:** SnapGene, FPbase, Swiss-Prot, and Rfam sequence-search hits. + These remain separate sequence-derived enrichment candidates. +- **LLM:** Prompt-ready descriptions generated by a model would be inference + only and must retain the model, prompt, evidence, and confidence. For example, `CDS` says that a region encodes a protein or peptide. It does not say whether that product is Cas9, a localization signal, an epitope tag, or a @@ -49,6 +53,39 @@ selection marker. The inventory therefore preserves both dimensions: } ``` +The generated inventory also retains evidence for every rule-derived role. A +role records whether it came from the GenBank feature type or from a particular +qualifier value and matched term. Regulatory features such as an `AmpR +promoter` are not themselves classified as selectable markers merely because +their label contains `AmpR`. + +## Generate everything for one plasmid + +The complete workflow is available as one command: + +```bash +uv run --extra lab_bench_2 python tools/report_cloning_sequence.py \ + path/to/plasmid.gbk \ + --output-dir plasmid-report +``` + +It creates: + +- `inventory.json`: source features, all qualifiers, rule-derived functional + summaries with per-role evidence, primers, enzyme sites, and reverse indexes. +- `external-sources.json`: cached iGEM and REBASE data with request URLs, + timestamps, response hashes, and selected specific-part evidence. +- `plannotate/`: separate pLannotate GBK/CSV output and its provenance manifest. +- `report.html`: a self-contained, color-coded report that can be opened + directly in a browser. +- `report-manifest.json`: source and artifact hashes plus the provenance policy. + +The first run automatically prepares the pinned pLannotate environment and +databases; later runs reuse validated caches. Use `--skip-plannotate`, +`--skip-igem`, or `--skip-rebase` only when deliberately requesting a partial +report. `--full-plannotate` adds Swiss-Prot and Rfam to the default fast +SnapGene/FPbase pass. + ## Fetch the independent source data No Addgene login or API token is involved in this workflow. Fetch the public @@ -105,7 +142,8 @@ The JSON contains: parser warnings, and errors. - A separate parts list for each non-source, non-primer feature: structural type, coordinates, strand, raw qualifiers, normalized label, deterministic - functional roles, and a conservative prompt-ready description. + functional roles, exact rule evidence, and a conservative prompt-ready + description. - A primer list with labels, coordinates, strand, and binding sequence in 5-prime-to-3-prime orientation. - Reverse indexes from feature label, functional role, primer label, sequence diff --git a/src/lab_bench_2/cloning_inventory.py b/src/lab_bench_2/cloning_inventory.py index 093ea14..479f820 100644 --- a/src/lab_bench_2/cloning_inventory.py +++ b/src/lab_bench_2/cloning_inventory.py @@ -470,7 +470,8 @@ def _inventory_record(record: Any, *, include_enzymes: bool) -> dict[str, Any]: def _inventory_feature(record: Any, feature: Any) -> dict[str, Any]: qualifiers = _clean_qualifiers(feature.qualifiers) label = _first_qualifier(qualifiers, LABEL_QUALIFIERS) - roles = classify_feature_roles(feature.type, qualifiers) + role_evidence = classify_feature_role_evidence(feature.type, qualifiers) + roles = [item["role"] for item in role_evidence] sequence = str(feature.extract(record.seq)).upper() return { "feature_type": str(feature.type), @@ -484,6 +485,7 @@ def _inventory_feature(record: Any, feature: Any) -> dict[str, Any]: "feature_sequence_sha256": hashlib.sha256(sequence.encode("ascii")).hexdigest(), "feature_signature": feature_signature(label, sequence), "functional_roles": roles, + "functional_role_evidence": role_evidence, "functional_description": describe_feature(label, roles), "qualifiers": qualifiers, } @@ -557,15 +559,31 @@ def classify_feature_roles( feature_type: str, qualifiers: Mapping[str, list[str]] ) -> list[str]: """Map a structural GenBank feature to zero or more functional roles.""" + return [ + item["role"] + for item in classify_feature_role_evidence(feature_type, qualifiers) + ] + + +def classify_feature_role_evidence( + feature_type: str, qualifiers: Mapping[str, list[str]] +) -> list[dict[str, Any]]: + """Map a GenBank feature to roles and retain the exact rule evidence.""" feature_type_lower = feature_type.casefold() - text = " ".join( - value for key in TEXT_QUALIFIERS for value in qualifiers.get(key, []) - ).casefold() - roles: set[str] = set() + evidence_by_role: defaultdict[str, list[dict[str, str]]] = defaultdict(list) + + def add_structural(role: str) -> None: + evidence_by_role[role].append( + { + "method": "genbank_feature_type", + "feature_type": feature_type, + } + ) + if feature_type_lower == "cds": - roles.add("coding_sequence") + add_structural("coding_sequence") if feature_type_lower in PRIMER_TYPES: - roles.add("primer_binding_site") + add_structural("primer_binding_site") structural_roles = { "promoter": "promoter", "enhancer": "enhancer", @@ -577,11 +595,34 @@ def classify_feature_roles( } role = structural_roles.get(feature_type_lower) if role: - roles.add(role) + add_structural(role) for candidate, patterns in ROLE_PATTERNS: - if any(_contains_term(text, pattern) for pattern in patterns): - roles.add(candidate) - return sorted(roles) + if candidate == "selection_marker" and feature_type_lower in { + "enhancer", + "intron", + "polya_signal", + "promoter", + "protein_bind", + "rep_origin", + "terminator", + }: + continue + for qualifier in TEXT_QUALIFIERS: + for source_value in qualifiers.get(qualifier, []): + for pattern in patterns: + if _contains_term(source_value.casefold(), pattern): + evidence_by_role[candidate].append( + { + "method": "curated_qualifier_pattern", + "qualifier": qualifier, + "source_value": source_value, + "matched_term": pattern, + } + ) + return [ + {"role": role_name, "evidence": evidence_by_role[role_name]} + for role_name in sorted(evidence_by_role) + ] def describe_feature(label: str, roles: Sequence[str]) -> str: diff --git a/src/lab_bench_2/cloning_report.py b/src/lab_bench_2/cloning_report.py new file mode 100644 index 0000000..b710787 --- /dev/null +++ b/src/lab_bench_2/cloning_report.py @@ -0,0 +1,506 @@ +"""Render one provenance-separated cloning inventory as standalone HTML.""" + +from __future__ import annotations + +import csv +import html +import json +from pathlib import Path +from typing import Any, Mapping, Sequence + + +class CloningReportError(ValueError): + """Raised when inventory data cannot describe one plasmid report.""" + + +def render_cloning_report( + source_path: Path, + inventory: Mapping[str, Any], + external_sources: Mapping[str, Any], + *, + plannotate_manifest: Mapping[str, Any] | None = None, + plannotate_output_dir: Path | None = None, +) -> str: + """Return a self-contained report for a single-record plasmid file.""" + file_entry, record = _single_record(inventory) + features = _mapping_list(record.get("features")) + primers = _mapping_list(record.get("primers")) + feature_colors = [_series_color(index) for index in range(len(features))] + primer_colors = [ + _series_color(index + len(features)) for index in range(len(primers)) + ] + sites = _string_int_mapping(record.get("neb_restriction_sites")) + zero_cutters = sorted(name for name, count in sites.items() if count == 0) + single_cutters = sorted(name for name, count in sites.items() if count == 1) + multi_cutters = sorted(name for name, count in sites.items() if count > 1) + selected_parts = [_selected_part(feature) for feature in features] + selected_count = sum(part is not None for part in selected_parts) + plannotate_rows = _plannotate_rows( + source_path, + plannotate_manifest, + plannotate_output_dir, + ) + + feature_rows = "".join( + _feature_row(index, feature, feature_colors[index]) + for index, feature in enumerate(features) + ) + primer_rows = "".join( + _primer_row(index, primer, primer_colors[index]) + for index, primer in enumerate(primers) + ) + igem_rows = "".join( + _igem_row(index, feature, selected_parts[index], feature_colors[index]) + for index, feature in enumerate(features) + ) + plannotate_rows_html = "".join( + _plannotate_row(index, row) for index, row in enumerate(plannotate_rows) + ) + feature_data = [ + { + "label": str(feature.get("label") or "unnamed feature"), + "type": str(feature.get("feature_type") or ""), + "start": int(feature.get("start_0_based") or 0), + "end": int(feature.get("end_0_based_exclusive") or 0), + "strand": feature.get("strand"), + "color": feature_colors[index], + } + for index, feature in enumerate(features) + ] + primer_data = [ + { + "label": str(primer.get("label") or "unnamed primer"), + "start": int(primer.get("start_0_based") or 0), + "end": int(primer.get("end_0_based_exclusive") or 0), + "strand": primer.get("strand"), + "color": primer_colors[index], + } + for index, primer in enumerate(primers) + ] + source_provenance = _source_provenance(external_sources) + plannotate_status = _plannotate_status(plannotate_manifest) + file_name = source_path.name + file_sha256 = str(file_entry.get("file_sha256") or "") + sequence_length = int(record.get("sequence_length") or 0) + topology = str(record.get("topology") or "unknown") + + body = f"""
+

Plasmid evidence report

+

{_esc(file_name)} · SHA-256 {_esc(file_sha256)}

+
+
{sequence_length:,} bp{_esc(topology)} DNA
+
{len(features)}source features
+
{len(primers)}source primers
+
{selected_count}/{len(features)}specific iGEM parts selected
+
+
+

Downloaded file · source truth

Sequence, topology, feature types, coordinates, labels, qualifiers, translations, and annotated primers.

+ +

Deterministic computation

Biopython parsing, rule-derived functional summaries with rule evidence, primer sequences, restriction sites, and hashes.

+ +

External evidence · separate

Specific iGEM records and roles, pLannotate sequence-search candidates, and the current REBASE supplier-N catalog.

+
+

Every feature and primer instance has a stable color shared by the map and its table row. Coordinates are 0-based and end-exclusive.

+
+ +
+

Feature type, location, label, and qualifiers are copied from the source file. “Rule-derived functional summary” is generated from explicit local rules; each row exposes the rule evidence.

+
{feature_rows or _empty_row(8, "No source features")}
#LabelGenBank typeLocationSource descriptionRule-derived functional summaryRule evidenceAll source qualifiers
+
+ + + + + + +
""" + return _standalone_html(body) + + +def write_cloning_report( + output_path: Path, + source_path: Path, + inventory: Mapping[str, Any], + external_sources: Mapping[str, Any], + *, + plannotate_manifest: Mapping[str, Any] | None = None, + plannotate_output_dir: Path | None = None, +) -> None: + """Write a standalone report, creating its parent directory.""" + report = render_cloning_report( + source_path, + inventory, + external_sources, + plannotate_manifest=plannotate_manifest, + plannotate_output_dir=plannotate_output_dir, + ) + output_path.parent.mkdir(parents=True, exist_ok=True) + temporary = output_path.with_suffix(output_path.suffix + ".tmp") + temporary.write_text(report, encoding="utf-8") + temporary.replace(output_path) + + +def _single_record( + inventory: Mapping[str, Any], +) -> tuple[Mapping[str, Any], Mapping[str, Any]]: + files = _mapping_list(inventory.get("files")) + if len(files) != 1: + raise CloningReportError("Report requires exactly one parsed plasmid file") + records = _mapping_list(files[0].get("records")) + if len(records) != 1: + raise CloningReportError("Report requires exactly one sequence record") + return files[0], records[0] + + +def _feature_row(index: int, feature: Mapping[str, Any], color: str) -> str: + qualifiers = feature.get("qualifiers") + qualifiers_mapping = qualifiers if isinstance(qualifiers, Mapping) else {} + evidence = _role_evidence_html(feature) + return ( + "" + f"{_swatch(color)}{index + 1}" + f"{_esc(feature.get('label') or 'unnamed feature')}" + f"{_esc(feature.get('feature_type') or '')}" + f"{_esc(feature.get('location') or '')}" + f"{_esc(_source_description(qualifiers_mapping))}" + f"{_esc(feature.get('functional_description') or '')}" + f"{evidence}" + f"{_json_details(qualifiers_mapping)}" + "" + ) + + +def _primer_row(index: int, primer: Mapping[str, Any], color: str) -> str: + return ( + "" + f"{_swatch(color)}{index + 1}" + f"{_esc(primer.get('label') or 'unnamed primer')}" + f"{_esc(primer.get('location') or '')}" + f"{_esc(primer.get('binding_sequence_5to3') or '')}" + f"{_json_details(primer.get('qualifiers') or {})}" + "" + ) + + +def _igem_row( + index: int, + feature: Mapping[str, Any], + part: Mapping[str, Any] | None, + color: str, +) -> str: + if part is None: + part_cell = "No validated specific part" + role_cell = "—" + description_cell = "—" + evidence_cell = "not selected" + else: + url = str(part.get("url") or "") + part_cell = ( + f'' + f"{_esc(part.get('name') or '')}: {_esc(part.get('title') or '')}" + ) + role = part.get("role") + role_mapping = role if isinstance(role, Mapping) else {} + role_cell = ( + f"{_esc(role_mapping.get('label') or '')} " + f"({_esc(role_mapping.get('accession') or '')})" + ) + description_cell = _esc( + html.unescape(str(part.get("description") or "No description supplied")) + ) + evidence_cell = _esc(_part_evidence(part)) + return ( + "" + f"{_swatch(color)}{index + 1}" + f"{_esc(feature.get('label') or 'unnamed feature')}" + f"{_esc(feature.get('feature_type') or '')}" + f"{_esc(feature.get('functional_description') or '')}" + f"{part_cell}{role_cell}" + f"{description_cell}{evidence_cell}" + "" + ) + + +def _plannotate_row(index: int, row: Mapping[str, Any]) -> str: + return ( + "" + f"{index + 1}" + f"{_esc(row.get('Feature') or row.get('feature') or '')}" + f"{_esc(row.get('Type') or row.get('feature_type') or '')}" + f"{_esc(_plannotate_location(row))}" + f"{_percent(row.get('percent identity', row.get('percent_identity')))}" + f"{_percent(row.get('percent match length', row.get('percent_match_length')))}" + f"{_esc(row.get('fragment') or '')}" + f"{_esc(row.get('database') or '')}" + f"{_esc(row.get('Description') or row.get('description') or '')}" + "" + ) + + +def _selected_part(feature: Mapping[str, Any]) -> Mapping[str, Any] | None: + external = feature.get("external_function_candidates") + external_mapping = external if isinstance(external, Mapping) else {} + igem = external_mapping.get("igem_registry") + igem_mapping = igem if isinstance(igem, Mapping) else {} + for part in _mapping_list(igem_mapping.get("specific_parts")): + if part.get("selected"): + return part + return None + + +def _role_evidence_html(feature: Mapping[str, Any]) -> str: + entries = _mapping_list(feature.get("functional_role_evidence")) + if not entries: + return "No role assigned" + items: list[str] = [] + for entry in entries: + role = _esc(entry.get("role") or "") + evidence_list = _mapping_list(entry.get("evidence")) + fragments: list[str] = [] + for evidence in evidence_list: + if evidence.get("method") == "genbank_feature_type": + fragments.append( + f"GBK type {_esc(evidence.get('feature_type') or '')}" + ) + else: + fragments.append( + f"/{_esc(evidence.get('qualifier') or '')} " + f"matched {_esc(evidence.get('matched_term') or '')}" + ) + items.append(f"
  • {role} ← {'; '.join(fragments)}
    • ") + return f"
      {''.join(items)}
    " + + +def _plannotate_rows( + source_path: Path, + manifest: Mapping[str, Any] | None, + output_dir: Path | None, +) -> list[Mapping[str, Any]]: + if manifest is None: + return [] + results = _mapping_list(manifest.get("results")) + result = next( + ( + item + for item in results + if Path(str(item.get("source_path") or "")).resolve() + == source_path.resolve() + ), + results[0] if len(results) == 1 else None, + ) + if result is None: + return [] + if output_dir is not None: + for output in result.get("outputs") or []: + name = str(output) + if Path(name).suffix.casefold() != ".csv": + continue + csv_path = output_dir / name + if csv_path.is_file(): + with csv_path.open(encoding="utf-8", newline="") as handle: + return [dict(row) for row in csv.DictReader(handle)] + summary = result.get("annotation_summary") + summary_mapping = summary if isinstance(summary, Mapping) else {} + return _mapping_list(summary_mapping.get("annotations")) + + +def _plannotate_status(manifest: Mapping[str, Any] | None) -> str: + if manifest is None: + return "not run" + tool = manifest.get("tool") + tool_mapping = tool if isinstance(tool, Mapping) else {} + summary = manifest.get("summary") + summary_mapping = summary if isinstance(summary, Mapping) else {} + return ( + f"pLannotate {tool_mapping.get('version') or 'unknown'}; " + f"{summary_mapping.get('annotated_or_cached_count') or 0} annotated/cached; " + f"{summary_mapping.get('error_count') or 0} errors" + ) + + +def _source_provenance(external_sources: Mapping[str, Any]) -> dict[str, str]: + sources = external_sources.get("sources") + sources_mapping = sources if isinstance(sources, Mapping) else {} + igem = sources_mapping.get("igem_registry") + rebase = sources_mapping.get("rebase_neb") + igem_mapping = igem if isinstance(igem, Mapping) else {} + rebase_mapping = rebase if isinstance(rebase, Mapping) else {} + igem_source = igem_mapping.get("source") + rebase_source = rebase_mapping.get("source") + igem_source_mapping = igem_source if isinstance(igem_source, Mapping) else {} + rebase_source_mapping = rebase_source if isinstance(rebase_source, Mapping) else {} + return { + "igem": ( + f"API {igem_source_mapping.get('api_version') or 'not fetched'}; " + f"anonymous published records; response provenance retained in external-sources.json" + ), + "rebase": ( + f"release {rebase_source_mapping.get('release') or 'not fetched'} " + f"({rebase_source_mapping.get('release_date') or 'unknown date'}); " + "response provenance retained in external-sources.json" + ), + } + + +def _part_evidence(part: Mapping[str, Any]) -> str: + evidence = part.get("evidence") + evidence_mapping = evidence if isinstance(evidence, Mapping) else {} + values: list[str] = [] + if evidence_mapping.get("nucleotide_exact"): + values.append("exact DNA") + else: + identity = evidence_mapping.get("same_length_nucleotide_identity_percent") + if isinstance(identity, (int, float)): + values.append(f"{identity:.3f}% same-length DNA identity") + if evidence_mapping.get("translated_peptide_exact"): + values.append("exact translated peptide") + return " · ".join(values) or "not sequence-verified" + + +def _source_description(qualifiers: Mapping[str, Any]) -> str: + for key in ("note", "product", "bound_moiety", "function"): + values = qualifiers.get(key) + if isinstance(values, Sequence) and not isinstance(values, (str, bytes)): + if values: + return str(values[0]) + return "—" + + +def _series_color(index: int) -> str: + base = [f"var(--series-{number})" for number in range(1, 7)] + if index < len(base): + return base[index] + first = index % 6 + second = (index * 5 + 3) % 6 + if first == second: + second = (second + 1) % 6 + ratio = 35 + (index * 13) % 46 + return f"color-mix(in srgb, {base[first]} {ratio}%, {base[second]})" + + +def _standalone_html(body: str) -> str: + return f""" + +Plasmid evidence report +{body}""" + + +def _mapping_list(value: object) -> list[Mapping[str, Any]]: + if not isinstance(value, Sequence) or isinstance(value, (str, bytes)): + return [] + return [item for item in value if isinstance(item, Mapping)] + + +def _string_int_mapping(value: object) -> dict[str, int]: + if not isinstance(value, Mapping): + return {} + return { + str(key): int(count) for key, count in value.items() if isinstance(count, int) + } + + +def _script_json(value: object) -> str: + return ( + json.dumps(value, ensure_ascii=False, separators=(",", ":")) + .replace("<", "\\u003c") + .replace(">", "\\u003e") + .replace("&", "\\u0026") + ) + + +def _json_details(value: object) -> str: + payload = _esc(json.dumps(value, indent=2, sort_keys=True, ensure_ascii=False)) + return f"
    View
    {payload}
    " + + +def _plannotate_location(row: Mapping[str, Any]) -> str: + start = row.get("start location") + end = row.get("end location") + return f"{start}–{end}" if start is not None and end is not None else "—" + + +def _percent(value: object) -> str: + try: + return f"{float(str(value)):.3g}%" + except (TypeError, ValueError): + return "—" + + +def _swatch(color: str) -> str: + return f'' + + +def _chip_list(values: Sequence[str]) -> str: + return ( + " ".join(f'{_esc(value)}' for value in values) or "—" + ) + + +def _empty_row(columns: int, text: str) -> str: + return f'{_esc(text)}' + + +def _esc(value: object) -> str: + return html.escape(str(value), quote=True) diff --git a/tests/lab_bench_2/test_cloning_inventory.py b/tests/lab_bench_2/test_cloning_inventory.py index 7df61de..6ab67bf 100644 --- a/tests/lab_bench_2/test_cloning_inventory.py +++ b/tests/lab_bench_2/test_cloning_inventory.py @@ -7,6 +7,7 @@ from lab_bench_2.cloning_inventory import ( build_cloning_inventory, + classify_feature_role_evidence, classify_feature_roles, feature_signature, ) @@ -110,6 +111,22 @@ def test_feature_type_and_function_are_separate() -> None: assert roles == ["coding_sequence", "epitope_tag", "localization_signal"] +def test_role_evidence_names_exact_source_rule_and_avoids_promoter_marker() -> None: + evidence = classify_feature_role_evidence( + "CDS", {"label": ["AmpR"], "note": ["confers ampicillin resistance"]} + ) + selection = next(item for item in evidence if item["role"] == "selection_marker") + assert selection["evidence"][0] == { + "method": "curated_qualifier_pattern", + "qualifier": "label", + "source_value": "AmpR", + "matched_term": "ampr", + } + assert classify_feature_roles("promoter", {"label": ["AmpR promoter"]}) == [ + "promoter" + ] + + def test_external_igem_specific_match_preserves_evidence_and_indexes( tmp_path: Path, ) -> None: diff --git a/tests/lab_bench_2/test_cloning_report.py b/tests/lab_bench_2/test_cloning_report.py new file mode 100644 index 0000000..427af8a --- /dev/null +++ b/tests/lab_bench_2/test_cloning_report.py @@ -0,0 +1,164 @@ +from __future__ import annotations + +from pathlib import Path + +import pytest + +from lab_bench_2.cloning_report import CloningReportError, render_cloning_report + + +def _inventory() -> dict[str, object]: + return { + "files": [ + { + "path": "plasmid.gbk", + "file_sha256": "file-digest", + "parse_warnings": [], + "records": [ + { + "sequence_length": 100, + "sequence_sha256": "sequence-digest", + "topology": "circular", + "features": [ + { + "label": "FLAG", + "feature_type": "CDS", + "location": "[10:34](+)", + "start_0_based": 10, + "end_0_based_exclusive": 34, + "strand": 1, + "functional_description": ( + "FLAG: protein-coding sequence, " + "epitope or detection tag" + ), + "functional_role_evidence": [ + { + "role": "coding_sequence", + "evidence": [ + { + "method": "genbank_feature_type", + "feature_type": "CDS", + } + ], + }, + { + "role": "epitope_tag", + "evidence": [ + { + "method": "curated_qualifier_pattern", + "qualifier": "label", + "source_value": "FLAG", + "matched_term": "flag", + } + ], + }, + ], + "qualifiers": { + "label": ["FLAG"], + "product": ["FLAG epitope tag"], + }, + "external_function_candidates": { + "igem_registry": { + "specific_parts": [ + { + "name": "BBa_K4587111", + "title": "FLAG Tag", + "url": ( + "https://registry.igem.org/parts/" + "bba-k4587111" + ), + "description": "A detection tag.", + "selected": True, + "role": { + "label": "Tag", + "accession": "SO:0000324", + }, + "evidence": { + "nucleotide_exact": False, + "same_length_nucleotide_identity_percent": 83.333, + "translated_peptide_exact": True, + }, + } + ] + } + }, + } + ], + "primers": [ + { + "label": "Primer A", + "location": "[1:9](+)", + "start_0_based": 1, + "end_0_based_exclusive": 9, + "strand": 1, + "binding_sequence_5to3": "ATGCGTAA", + "qualifiers": {"note": ["source primer"]}, + } + ], + "neb_restriction_sites": {"EcoRI": 1, "BamHI": 0}, + } + ], + } + ] + } + + +def test_report_separates_source_rules_and_specific_igem_data(tmp_path: Path) -> None: + source = tmp_path / "plasmid.gbk" + source.write_text("source is represented by inventory fixture") + manifest = { + "tool": {"version": "2.0.0"}, + "summary": {"annotated_or_cached_count": 1, "error_count": 0}, + "results": [ + { + "source_path": str(source), + "annotation_summary": { + "annotations": [ + { + "feature": "FLAG", + "feature_type": "CDS", + "description": "database candidate", + "database": "snapgene", + "fragment": False, + "percent_identity": 100.0, + "percent_match_length": 100.0, + } + ] + }, + } + ], + } + + report = render_cloning_report( + source, + _inventory(), + {"sources": {}}, + plannotate_manifest=manifest, + ) + + assert "Rule-derived functional summary" in report + assert "epitope_tag" in report + assert "matched flag" in report + assert "BBa_K4587111: FLAG Tag" in report + assert "Tag (SO:0000324)" in report + assert "83.333% same-length DNA identity · exact translated peptide" in report + assert "A detection tag." in report + assert "Primer A" in report + assert "EcoRI" in report + assert "pLannotate 2.0.0; 1 annotated/cached; 0 errors" in report + assert "background:var(--series-1)" in report + assert "background:var(--series-2)" in report + + +def test_report_rejects_more_than_one_record(tmp_path: Path) -> None: + inventory = _inventory() + files = inventory["files"] + assert isinstance(files, list) + file_entry = files[0] + assert isinstance(file_entry, dict) + records = file_entry["records"] + assert isinstance(records, list) + records.append(dict(records[0])) + + with pytest.raises(CloningReportError, match="exactly one sequence record"): + render_cloning_report(tmp_path / "multi.gbk", inventory, {"sources": {}}) diff --git a/tools/report_cloning_sequence.py b/tools/report_cloning_sequence.py new file mode 100644 index 0000000..2a3caf4 --- /dev/null +++ b/tools/report_cloning_sequence.py @@ -0,0 +1,213 @@ +#!/usr/bin/env python3 +"""Build a complete provenance-separated report for one plasmid file.""" + +from __future__ import annotations + +import argparse +import hashlib +import json +import sys +from datetime import UTC, datetime +from pathlib import Path +from typing import Any + +from lab_bench_2.cloning_external_sources import ( + ExternalSourceError, + fetch_cloning_external_sources, +) +from lab_bench_2.cloning_inventory import ( + SEQUENCE_FORMATS, + build_cloning_inventory, + collect_feature_queries, +) +from lab_bench_2.cloning_report import CloningReportError, write_cloning_report +from lab_bench_2.plannotate_runner import ( + DEFAULT_CACHE_DIR, + PlannotateError, + annotate_with_plannotate, +) + + +def _arguments() -> argparse.Namespace: + parser = argparse.ArgumentParser( + description=( + "Turn one GenBank/FASTA plasmid into inventory JSON, cached public " + "iGEM/REBASE evidence, optional pLannotate outputs, and a standalone " + "color-coded HTML report." + ) + ) + parser.add_argument("input", type=Path, help="One single-record plasmid file") + parser.add_argument("--output-dir", type=Path, required=True) + parser.add_argument( + "--source-cache-dir", + type=Path, + default=Path.home() / ".cache/labbench2/cloning-sources", + ) + parser.add_argument("--plannotate-cache-dir", type=Path, default=DEFAULT_CACHE_DIR) + parser.add_argument("--cores", type=int, default=4) + parser.add_argument("--refresh", action="store_true") + parser.add_argument("--full-plannotate", action="store_true") + parser.add_argument("--skip-plannotate", action="store_true") + parser.add_argument("--skip-igem", action="store_true") + parser.add_argument("--skip-rebase", action="store_true") + parser.add_argument( + "--no-auto-setup", + action="store_true", + help="Fail instead of automatically installing pLannotate and its databases", + ) + return parser.parse_args() + + +def main() -> int: + args = _arguments() + source = args.input.expanduser().resolve() + output_dir = args.output_dir.expanduser().resolve() + try: + _validate_source(source) + if args.cores < 1: + raise ValueError("cores must be positive") + output_dir.mkdir(parents=True, exist_ok=True) + external_sources = fetch_cloning_external_sources( + args.source_cache_dir, + refresh=args.refresh, + include_rebase=not args.skip_rebase, + include_igem=not args.skip_igem, + feature_queries=( + collect_feature_queries([source]) if not args.skip_igem else [] + ), + ) + external_path = output_dir / "external-sources.json" + _write_json(external_path, external_sources) + inventory = build_cloning_inventory( + [source], + root=source.parent, + external_sources=external_sources, + ) + inventory_path = output_dir / "inventory.json" + _write_json(inventory_path, inventory) + + plannotate_manifest: dict[str, Any] | None = None + plannotate_dir: Path | None = None + if not args.skip_plannotate: + plannotate_dir = output_dir / "plannotate" + plannotate_manifest = annotate_with_plannotate( + [source], + plannotate_dir, + cache_dir=args.plannotate_cache_dir, + fast=not args.full_plannotate, + cores=args.cores, + refresh=args.refresh, + auto_setup=not args.no_auto_setup, + ) + + report_path = output_dir / "report.html" + write_cloning_report( + report_path, + source, + inventory, + external_sources, + plannotate_manifest=plannotate_manifest, + plannotate_output_dir=plannotate_dir, + ) + pipeline_manifest = _pipeline_manifest( + source, + output_dir, + inventory_path, + external_path, + report_path, + plannotate_manifest, + ) + _write_json(output_dir / "report-manifest.json", pipeline_manifest) + except ( + CloningReportError, + ExternalSourceError, + OSError, + PlannotateError, + ValueError, + ) as exc: + print(f"error: {exc}", file=sys.stderr) + return 2 + + summary = inventory["summary"] + plannotate_errors = ( + int(plannotate_manifest["summary"]["error_count"]) + if plannotate_manifest is not None + else 0 + ) + print( + f"Reported {summary['parsed_file_count']}/1 plasmid file; " + f"{plannotate_errors} pLannotate error(s). Report: {report_path}" + ) + return 1 if plannotate_errors else 0 + + +def _validate_source(source: Path) -> None: + if not source.is_file(): + raise ValueError(f"Input is not a file: {source}") + if source.suffix.casefold() not in SEQUENCE_FORMATS: + supported = ", ".join(sorted(SEQUENCE_FORMATS)) + raise ValueError(f"Unsupported input suffix; expected one of: {supported}") + + +def _pipeline_manifest( + source: Path, + output_dir: Path, + inventory_path: Path, + external_path: Path, + report_path: Path, + plannotate_manifest: dict[str, Any] | None, +) -> dict[str, Any]: + artifacts = { + "inventory": _artifact(output_dir, inventory_path), + "external_sources": _artifact(output_dir, external_path), + "report": _artifact(output_dir, report_path), + } + if plannotate_manifest is not None: + plannotate_manifest_path = output_dir / "plannotate/plannotate-manifest.json" + artifacts["plannotate_manifest"] = _artifact( + output_dir, plannotate_manifest_path + ) + return { + "schema_version": 1, + "generated_at": datetime.now(UTC).isoformat(), + "source": { + "path": str(source), + "sha256": _sha256_file(source), + }, + "artifacts": artifacts, + "provenance_policy": { + "source_file_modified": False, + "source_annotations_are_ground_truth": True, + "functional_summaries_are_rule_derived": True, + "external_matches_remain_separate": True, + }, + } + + +def _artifact(root: Path, path: Path) -> dict[str, Any]: + return { + "path": str(path.relative_to(root)), + "sha256": _sha256_file(path), + "size_bytes": path.stat().st_size, + } + + +def _sha256_file(path: Path) -> str: + digest = hashlib.sha256() + with path.open("rb") as handle: + for chunk in iter(lambda: handle.read(1024 * 1024), b""): + digest.update(chunk) + return digest.hexdigest() + + +def _write_json(path: Path, payload: dict[str, Any]) -> None: + path.parent.mkdir(parents=True, exist_ok=True) + temporary = path.with_suffix(path.suffix + ".tmp") + temporary.write_text( + json.dumps(payload, indent=2, sort_keys=True) + "\n", encoding="utf-8" + ) + temporary.replace(path) + + +if __name__ == "__main__": + raise SystemExit(main())