Skills › Scientific Agent Skills › gtars › refget_digest_plan.py
refget_digest_plan.py
311 lines · skills/gtars/scripts/refget_digest_plan.py at 5f2bc03
1#!/usr/bin/env python3
2"""Validate local refget metadata and FASTA sequence digests without networking."""
4from __future__ import annotations
6import argparse
7import base64
8import hashlib
9import re
10import sys
11from pathlib import Path
13from _common import (
14 HARD_MAX_BYTES,
15 HARD_MAX_RECORDS,
16 SafetyError,
17 add_path_mode_argument,
18 display_path,
19 fail_json,
20 int_type,
21 iter_text_lines,
22 load_json,
23 local_path,
24 print_json,
25 sha256_file,
26)
29TOOL = "gtars-refget-digest-plan"
30_SHA512T24U = re.compile(r"^[A-Za-z0-9_-]{32}$")
31_MD5 = re.compile(r"^[0-9a-f]{32}$")
34def build_parser() -> argparse.ArgumentParser:
35 parser = argparse.ArgumentParser(
36 description=(
37 "Validate a local refget metadata JSON document and optionally compare "
38 "it with a local FASTA. No store is opened and no endpoint is contacted."
39 )
40 )
41 parser.add_argument("--metadata", required=True)
42 parser.add_argument("--fasta", help="Optional local FASTA or FASTA.GZ.")
43 parser.add_argument("--assembly", required=True)
44 parser.add_argument(
45 "--max-bytes",
46 type=int_type(minimum=1, maximum=HARD_MAX_BYTES, label="max-bytes"),
47 default=2 * 1024**3,
48 )
49 parser.add_argument(
50 "--max-records",
51 type=int_type(
52 minimum=1,
53 maximum=HARD_MAX_RECORDS,
54 label="max-records",
55 ),
56 default=5_000_000,
57 help="Maximum FASTA text lines (default: 5,000,000).",
58 )
59 add_path_mode_argument(parser)
60 return parser
63def sha512t24u(sequence: bytes) -> str:
64 """Return the unprefixed GA4GH sha512t24u digest used by gtars."""
65 digest = hashlib.sha512(sequence.upper()).digest()[:24]
66 return base64.urlsafe_b64encode(digest).rstrip(b"=").decode("ascii")
69def _normalize_expected_digest(value: object) -> str | None:
70 if not isinstance(value, str):
71 return None
72 if value.lower().startswith("sq."):
73 value = value[3:]
74 return value if _SHA512T24U.fullmatch(value) else None
77def _read_fasta(
78 path: Path,
79 *,
80 max_bytes: int,
81 max_records: int,
82) -> list[dict[str, object]]:
83 records: list[dict[str, object]] = []
84 name: str | None = None
85 chunks: list[bytes] = []
86 names: set[str] = set()
88 def finish() -> None:
89 nonlocal name, chunks
90 if name is None:
91 return
92 sequence = b"".join(chunks).upper()
93 if not sequence:
94 raise SafetyError("FASTA contains an empty sequence")
95 records.append(
96 {
97 "name": name,
98 "length": len(sequence),
99 "sha512t24u": sha512t24u(sequence),
100 "md5": hashlib.md5(sequence).hexdigest(),
101 }
102 )
104 for line_number, line in iter_text_lines(
105 path,
106 max_bytes=max_bytes,
107 max_records=max_records,
108 ):
109 if line.startswith(">"):
110 finish()
111 header = line[1:].strip()
112 candidate = header.split(maxsplit=1)[0] if header else ""
113 if not candidate or len(candidate) > 1_000:
114 raise SafetyError(f"invalid FASTA name at line {line_number}")
115 if candidate in names:
116 raise SafetyError("duplicate FASTA sequence name")
117 names.add(candidate)
118 name = candidate
119 chunks = []
120 continue
121 if name is None:
122 if line.strip():
123 raise SafetyError("FASTA sequence data appears before the first header")
124 continue
125 compact = "".join(line.split())
126 try:
127 encoded = compact.encode("ascii")
128 except UnicodeEncodeError as exc:
129 raise SafetyError(
130 f"non-ASCII FASTA sequence near line {line_number}"
131 ) from exc
132 if any(byte < 33 or byte > 126 for byte in encoded):
133 raise SafetyError(f"invalid FASTA sequence byte near line {line_number}")
134 chunks.append(encoded)
135 finish()
136 if not records:
137 raise SafetyError("FASTA contains no sequences")
138 return records
141def validate(args: argparse.Namespace) -> tuple[dict, int]:
142 assembly = args.assembly.strip()
143 if not assembly or len(assembly) > 200:
144 raise SafetyError("assembly must contain 1-200 characters")
145 metadata_path = local_path(args.metadata, kind="file")
146 document = load_json(
147 metadata_path,
148 max_bytes=min(args.max_bytes, 64 * 1024**2),
149 )
150 if not isinstance(document, dict):
151 raise SafetyError("metadata root must be a JSON object")
153 errors: dict[str, int] = {}
154 warnings: dict[str, int] = {}
155 if document.get("schema_version") != "1.0":
156 errors["metadata:unsupported_schema_version"] = 1
157 if document.get("assembly") != assembly:
158 errors["metadata:assembly_mismatch"] = 1
159 if document.get("coordinate_system") != "0-based-half-open":
160 errors["metadata:coordinate_system_mismatch"] = 1
162 collection_digest = document.get("collection_digest")
163 if not isinstance(collection_digest, str) or not _SHA512T24U.fullmatch(
164 collection_digest
165 ):
166 errors["metadata:invalid_collection_digest_shape"] = 1
167 else:
168 warnings["collection_digest_shape_only_not_recomputed"] = 1
170 expected = document.get("sequences")
171 if not isinstance(expected, list) or not expected:
172 raise SafetyError("metadata.sequences must be a nonempty array")
173 if len(expected) > args.max_records:
174 raise SafetyError("metadata sequence count exceeds max-records")
176 expected_by_name: dict[str, dict] = {}
177 for item in expected:
178 if not isinstance(item, dict):
179 errors["metadata:non_object_sequence_record"] = (
180 errors.get("metadata:non_object_sequence_record", 0) + 1
181 )
182 continue
183 name = item.get("name")
184 length = item.get("length")
185 digest = _normalize_expected_digest(item.get("sha512t24u"))
186 md5 = item.get("md5")
187 if not isinstance(name, str) or not name or len(name) > 1_000:
188 errors["metadata:invalid_sequence_name"] = (
189 errors.get("metadata:invalid_sequence_name", 0) + 1
190 )
191 continue
192 if name in expected_by_name:
193 errors["metadata:duplicate_sequence_name"] = (
194 errors.get("metadata:duplicate_sequence_name", 0) + 1
195 )
196 continue
197 expected_by_name[name] = item
198 if type(length) is not int or length < 1:
199 errors["metadata:invalid_sequence_length"] = (
200 errors.get("metadata:invalid_sequence_length", 0) + 1
201 )
202 if digest is None:
203 errors["metadata:invalid_sha512t24u"] = (
204 errors.get("metadata:invalid_sha512t24u", 0) + 1
205 )
206 if not isinstance(md5, str) or not _MD5.fullmatch(md5.lower()):
207 errors["metadata:invalid_md5"] = errors.get("metadata:invalid_md5", 0) + 1
209 fasta_report = None
210 if args.fasta:
211 fasta_path = local_path(args.fasta, kind="file")
212 actual = _read_fasta(
213 fasta_path,
214 max_bytes=args.max_bytes,
215 max_records=args.max_records,
216 )
217 actual_by_name = {item["name"]: item for item in actual}
218 missing = set(expected_by_name) - set(actual_by_name)
219 unexpected = set(actual_by_name) - set(expected_by_name)
220 length_mismatch = 0
221 sha_mismatch = 0
222 md5_mismatch = 0
223 for name in set(expected_by_name) & set(actual_by_name):
224 expected_item = expected_by_name[name]
225 actual_item = actual_by_name[name]
226 if expected_item.get("length") != actual_item["length"]:
227 length_mismatch += 1
228 expected_sha = _normalize_expected_digest(
229 expected_item.get("sha512t24u")
230 )
231 if expected_sha != actual_item["sha512t24u"]:
232 sha_mismatch += 1
233 expected_md5 = expected_item.get("md5")
234 if (
235 not isinstance(expected_md5, str)
236 or expected_md5.lower() != actual_item["md5"]
237 ):
238 md5_mismatch += 1
239 for code, count in (
240 ("fasta:missing_sequences", len(missing)),
241 ("fasta:unexpected_sequences", len(unexpected)),
242 ("fasta:length_mismatches", length_mismatch),
243 ("fasta:sha512t24u_mismatches", sha_mismatch),
244 ("fasta:md5_mismatches", md5_mismatch),
245 ):
246 if count:
247 errors[code] = count
248 fasta_digest, fasta_bytes = sha256_file(
249 fasta_path,
250 max_bytes=args.max_bytes,
251 )
252 fasta_report = {
253 "path": display_path(fasta_path, 2, args.path_mode),
254 "sha256": fasta_digest,
255 "size_bytes": fasta_bytes,
256 "sequence_count": len(actual),
257 "all_sequence_digests_match": not any(
258 (missing, unexpected, length_mismatch, sha_mismatch, md5_mismatch)
259 ),
260 }
261 else:
262 warnings["fasta_not_supplied_digest_values_not_recomputed"] = 1
264 metadata_digest, metadata_bytes = sha256_file(
265 metadata_path,
266 max_bytes=min(args.max_bytes, 64 * 1024**2),
267 )
268 report = {
269 "ok": not errors,
270 "tool": TOOL,
271 "contract": {
272 "assembly": assembly,
273 "coordinate_system": "0-based-half-open",
274 "network_used": False,
275 "store_opened": False,
276 "cache_written": False,
277 "packages_imported": False,
278 "symlinks_allowed": False,
279 },
280 "metadata": {
281 "path": display_path(metadata_path, 1, args.path_mode),
282 "sha256": metadata_digest,
283 "size_bytes": metadata_bytes,
284 "sequence_count": len(expected_by_name),
285 "collection_digest_present": isinstance(collection_digest, str),
286 },
287 "fasta": fasta_report,
288 "errors": dict(sorted(errors.items())),
289 "warnings": dict(sorted(warnings.items())),
290 "approved_execution_plan": [
291 "prefer RefgetStore.open_local for a reviewed local store",
292 "for a new store use RefgetStore.in_memory, verify, then persist explicitly",
293 "before open_remote record an HTTPS allowlist, immutable revision, expected digests, cache path, byte quota, and approval",
294 "treat get_substring and stream_sequence on a remote store as network reads",
295 ],
296 }
297 return report, 0 if not errors else 2
300def main() -> int:
301 args = build_parser().parse_args()
302 try:
303 report, status = validate(args)
304 print_json(report)
305 return status
306 except (OSError, SafetyError, UnicodeError) as exc:
307 return fail_json(TOOL, exc)
310if __name__ == "__main__":
311 sys.exit(main())