Skills › Scientific Agent Skills › geniml › bed_validator.py
bed_validator.py
363 lines · skills/geniml/scripts/bed_validator.py at 5f2bc03
1#!/usr/bin/env python3
2"""Validate one local BED file and emit a non-mutating normalization plan."""
4from __future__ import annotations
6import argparse
7import re
8import sys
9from collections import Counter
10from pathlib import Path
12from _common import (
13 HARD_MAX_BYTES,
14 HARD_MAX_RECORDS,
15 MAX_COORDINATE,
16 SafetyError,
17 add_path_mode_argument,
18 display_path,
19 fail_json,
20 int_type,
21 iter_text_lines,
22 load_chrom_sizes,
23 local_path,
24 print_json,
25 sha256_file,
26)
29TOOL = "bed-validator"
30_DECIMAL = re.compile(r"^-?[0-9]+$")
31_VALID_STRANDS = {"+", "-", "."}
34def build_parser() -> argparse.ArgumentParser:
35 parser = argparse.ArgumentParser(
36 description=(
37 "Validate a local BED file as 0-based half-open intervals and emit "
38 "a bounded normalization plan. The input is never rewritten."
39 )
40 )
41 parser.add_argument("--input", required=True, help="Local BED or BED.GZ file.")
42 parser.add_argument(
43 "--assembly",
44 required=True,
45 help="Declared assembly/accession recorded in the report.",
46 )
47 parser.add_argument(
48 "--chrom-sizes",
49 help="Optional local two-column chromosome-sizes file for contig/bounds checks.",
50 )
51 parser.add_argument(
52 "--max-bytes",
53 type=int_type(minimum=1, maximum=HARD_MAX_BYTES, label="max-bytes"),
54 default=512 * 1024**2,
55 help="Maximum compressed and expanded input bytes (default: 512 MiB).",
56 )
57 parser.add_argument(
58 "--max-records",
59 type=int_type(minimum=1, maximum=HARD_MAX_RECORDS, label="max-records"),
60 default=1_000_000,
61 help="Maximum text records, including headers (default: 1,000,000).",
62 )
63 parser.add_argument(
64 "--max-examples",
65 type=int_type(minimum=0, maximum=100, label="max-examples"),
66 default=10,
67 help="Maximum redacted issue examples (default: 10).",
68 )
69 parser.add_argument(
70 "--allow-empty",
71 action="store_true",
72 help="Do not treat a BED with zero data records as an error.",
73 )
74 add_path_mode_argument(parser)
75 return parser
78def _record_issue(
79 counter: Counter[str],
80 examples: list[dict[str, int | str]],
81 code: str,
82 line_number: int,
83 maximum: int,
84) -> None:
85 counter[code] += 1
86 if len(examples) < maximum:
87 examples.append({"code": code, "line": line_number})
90def validate(args: argparse.Namespace) -> tuple[dict, int]:
91 if not args.assembly.strip() or len(args.assembly) > 200:
92 raise SafetyError("assembly must be a nonempty value of at most 200 characters")
94 bed_path = local_path(args.input, kind="file")
95 chrom_sizes: dict[str, int] | None = None
96 chrom_order: dict[str, int] | None = None
97 chrom_sizes_digest: str | None = None
98 if args.chrom_sizes:
99 chrom_path = local_path(args.chrom_sizes, kind="file")
100 sizes, order = load_chrom_sizes(
101 chrom_path,
102 max_bytes=min(args.max_bytes, 128 * 1024**2),
103 max_records=min(args.max_records, 1_000_000),
104 )
105 chrom_sizes = sizes
106 chrom_order = {chrom: index for index, chrom in enumerate(order)}
107 chrom_sizes_digest, _ = sha256_file(
108 chrom_path,
109 max_bytes=min(args.max_bytes, 128 * 1024**2),
110 )
112 errors: Counter[str] = Counter()
113 warnings: Counter[str] = Counter()
114 examples: list[dict[str, int | str]] = []
115 records = 0
116 skipped = 0
117 min_columns: int | None = None
118 max_columns = 0
119 contigs: set[str] = set()
120 seen: set[tuple[str, int, int]] = set()
121 duplicate_count = 0
122 overlap_count = 0
123 max_end_by_contig: dict[str, int] = {}
124 previous_sort_key: tuple[int | str, int, int] | None = None
125 unsorted_count = 0
126 strand_rows = 0
128 for line_number, line in iter_text_lines(
129 bed_path,
130 max_bytes=args.max_bytes,
131 max_records=args.max_records,
132 ):
133 if not line:
134 skipped += 1
135 warnings["blank_line"] += 1
136 continue
137 if line.startswith("#") or line.startswith("track") or line.startswith("browser"):
138 skipped += 1
139 continue
141 fields = line.split("\t")
142 if len(fields) < 3:
143 _record_issue(
144 errors,
145 examples,
146 "fewer_than_three_columns",
147 line_number,
148 args.max_examples,
149 )
150 continue
151 records += 1
152 min_columns = len(fields) if min_columns is None else min(min_columns, len(fields))
153 max_columns = max(max_columns, len(fields))
155 chrom = fields[0]
156 if not chrom or any(character.isspace() or ord(character) < 32 for character in chrom):
157 _record_issue(
158 errors,
159 examples,
160 "invalid_contig",
161 line_number,
162 args.max_examples,
163 )
164 continue
166 start_text, end_text = fields[1], fields[2]
167 if not _DECIMAL.fullmatch(start_text) or not _DECIMAL.fullmatch(end_text):
168 _record_issue(
169 errors,
170 examples,
171 "non_decimal_coordinate",
172 line_number,
173 args.max_examples,
174 )
175 continue
176 start, end = int(start_text), int(end_text)
177 if start < 0 or end < 0:
178 _record_issue(
179 errors,
180 examples,
181 "negative_coordinate",
182 line_number,
183 args.max_examples,
184 )
185 continue
186 if start > MAX_COORDINATE or end > MAX_COORDINATE:
187 _record_issue(
188 errors,
189 examples,
190 "coordinate_integer_overflow",
191 line_number,
192 args.max_examples,
193 )
194 continue
195 if end == start:
196 _record_issue(
197 errors,
198 examples,
199 "zero_length_interval",
200 line_number,
201 args.max_examples,
202 )
203 continue
204 if end < start:
205 _record_issue(
206 errors,
207 examples,
208 "end_before_start",
209 line_number,
210 args.max_examples,
211 )
212 continue
214 if chrom_sizes is not None:
215 if chrom not in chrom_sizes:
216 _record_issue(
217 errors,
218 examples,
219 "unknown_contig",
220 line_number,
221 args.max_examples,
222 )
223 continue
224 if end > chrom_sizes[chrom]:
225 _record_issue(
226 errors,
227 examples,
228 "end_beyond_contig",
229 line_number,
230 args.max_examples,
231 )
232 continue
234 if len(fields) >= 6:
235 strand_rows += 1
236 if fields[5] not in _VALID_STRANDS:
237 _record_issue(
238 errors,
239 examples,
240 "invalid_bed6_strand",
241 line_number,
242 args.max_examples,
243 )
245 contigs.add(chrom)
246 key = (chrom, start, end)
247 if key in seen:
248 duplicate_count += 1
249 else:
250 seen.add(key)
252 prior_end = max_end_by_contig.get(chrom)
253 if prior_end is not None and start < prior_end:
254 overlap_count += 1
255 max_end_by_contig[chrom] = max(end, prior_end or end)
257 if chrom_order is None:
258 sort_key: tuple[int | str, int, int] = (chrom, start, end)
259 else:
260 sort_key = (chrom_order[chrom], start, end)
261 if previous_sort_key is not None and sort_key < previous_sort_key:
262 unsorted_count += 1
263 previous_sort_key = sort_key
265 if records == 0 and not args.allow_empty:
266 errors["no_data_records"] += 1
267 if chrom_sizes is None:
268 warnings["bounds_not_checked_without_chrom_sizes"] += 1
269 if unsorted_count:
270 warnings["out_of_order_records"] = unsorted_count
271 if duplicate_count:
272 warnings["duplicate_intervals"] = duplicate_count
273 if overlap_count:
274 warnings["overlapping_intervals"] = overlap_count
276 digest, size = sha256_file(bed_path, max_bytes=args.max_bytes)
277 actions: list[dict[str, str | int]] = []
278 if errors:
279 actions.append(
280 {
281 "action": "reject_or_quarantine_invalid_rows",
282 "count": sum(errors.values()),
283 }
284 )
285 if unsorted_count:
286 actions.append(
287 {
288 "action": "stable_sort",
289 "detail": (
290 "use chromosome-sizes order, then numeric start/end"
291 if chrom_order is not None
292 else "supply chromosome sizes before choosing contig order"
293 ),
294 }
295 )
296 if duplicate_count:
297 actions.append(
298 {
299 "action": "decide_duplicate_policy",
300 "count": duplicate_count,
301 }
302 )
303 if chrom_sizes is None:
304 actions.append(
305 {
306 "action": "supply_chromosome_sizes",
307 "detail": "required before overflow/contig normalization",
308 }
309 )
310 actions.append(
311 {
312 "action": "preserve_coordinate_semantics",
313 "detail": "BED 0-based half-open; no implicit liftover or chr renaming",
314 }
315 )
317 report = {
318 "ok": not errors,
319 "tool": TOOL,
320 "contract": {
321 "assembly": args.assembly,
322 "coordinate_system": "0-based-half-open",
323 "input_mutated": False,
324 "network_used": False,
325 },
326 "input": {
327 "path": display_path(bed_path, 1, args.path_mode),
328 "sha256": digest,
329 "size_bytes": size,
330 },
331 "chromosome_sizes_sha256": chrom_sizes_digest,
332 "summary": {
333 "records": records,
334 "skipped_header_or_blank_lines": skipped,
335 "contig_count": len(contigs),
336 "minimum_columns": min_columns,
337 "maximum_columns": max_columns if records else None,
338 "bed6_or_wider_rows": strand_rows,
339 "duplicate_intervals": duplicate_count,
340 "overlap_observations": overlap_count,
341 "out_of_order_records": unsorted_count,
342 },
343 "errors": dict(sorted(errors.items())),
344 "warnings": dict(sorted(warnings.items())),
345 "issue_examples": examples,
346 "normalization_plan": actions,
347 }
348 return report, 0 if not errors else 2
351def main() -> int:
352 parser = build_parser()
353 args = parser.parse_args()
354 try:
355 report, exit_code = validate(args)
356 print_json(report)
357 return exit_code
358 except (OSError, SafetyError, UnicodeError) as exc:
359 return fail_json(TOOL, exc)
362if __name__ == "__main__":
363 sys.exit(main())