Skills › Scientific Agent Skills › genomic-coordinates › audit_intervals.py
audit_intervals.py
539 lines · skills/genomic-coordinates/scripts/audit_intervals.py at 5f2bc03
1#!/usr/bin/env python3
2"""Check an interval or variant file against its own format's conventions.
4A BED file holding 1-based coordinates parses cleanly, sorts cleanly, and
5intersects cleanly. Nothing downstream complains; every result is shifted by one
6base. These checks look for the evidence that survives that kind of mistake --
7coordinates a format cannot legally hold, features whose length contradicts
8their intent, positions past the end of the contig.
10 python3 audit_intervals.py peaks.bed
11 python3 audit_intervals.py gencode.gtf --genome hg38.chrom.sizes
12 python3 audit_intervals.py cohort.vcf --genome GRCh38.fa.fai
14Exit codes: 0 clean or warnings only, 1 at least one fatal finding,
152 usage error.
16"""
18from __future__ import annotations
20import argparse
21import re
22import sys
23from collections import Counter
24from pathlib import Path
26sys.path.insert(0, str(Path(__file__).resolve().parent))
28from _common import canonical_contig, emit, naming_style # noqa: E402
29from normalize_variant import special_alt # noqa: E402
31COLUMNS = ["severity", "rule", "line", "detail"]
32VALID_BASES = set("ACGTNacgtn")
35class Findings:
36 """Collected findings, capped per rule so one broken file cannot flood stdout."""
38 def __init__(self, max_examples: int) -> None:
39 self.rows: list[dict] = []
40 self.counts: Counter = Counter()
41 self.by_severity: Counter = Counter()
42 self.max_examples = max_examples
44 def add(self, severity: str, rule: str, line: int | str, detail: str) -> None:
45 self.counts[rule] += 1
46 self.by_severity[severity] += 1
47 if self.counts[rule] <= self.max_examples:
48 self.rows.append(
49 {"severity": severity, "rule": rule, "line": line, "detail": detail}
50 )
51 elif self.counts[rule] == self.max_examples + 1:
52 self.rows.append(
53 {
54 "severity": severity,
55 "rule": rule,
56 "line": "...",
57 "detail": "further occurrences suppressed (--max-examples)",
58 }
59 )
61 @property
62 def fatal(self) -> int:
63 """Every fatal occurrence, including those suppressed from the table."""
64 return self.by_severity["fatal"]
67def detect_format(path: Path) -> str:
68 suffixes = [s.lower() for s in path.suffixes]
69 for suffix, fmt in (
70 (".bed", "bed"), (".narrowpeak", "narrowpeak"), (".broadpeak", "broadpeak"),
71 (".bedgraph", "bedgraph"), (".gtf", "gtf"), (".gff3", "gff3"), (".gff", "gff3"),
72 (".vcf", "vcf"),
73 ):
74 if suffix in suffixes:
75 return fmt
76 raise SystemExit(
77 f"error: cannot infer the format of {path.name}; pass --format "
78 "bed|gtf|gff3|vcf"
79 )
82def load_genome(path: Path | None) -> dict[str, int]:
83 if path is None:
84 return {}
85 sizes = {}
86 for line in path.read_text().splitlines():
87 if not line.strip() or line.startswith("#"):
88 continue
89 parts = line.split("\t") if "\t" in line else line.split()
90 if len(parts) >= 2 and parts[1].isdigit():
91 sizes[canonical_contig(parts[0])] = int(parts[1])
92 return sizes
95def check_bounds(
96 find: Findings,
97 genome: dict[str, int],
98 lineno: int,
99 contig: str,
100 far: int,
101 source: str = "the genome file",
102) -> None:
103 """``far`` is the largest 1-based position the record touches."""
104 if not genome:
105 return
106 key = canonical_contig(contig)
107 if key not in genome:
108 find.add(
109 "fatal", "contig_not_in_genome", lineno,
110 f"{contig} is not in {source}, even ignoring the chr prefix",
111 )
112 elif far > genome[key]:
113 find.add(
114 "fatal", "past_contig_end", lineno,
115 f"{contig}:{far} is past the end of {contig} ({genome[key]} bp per "
116 f"{source}); wrong assembly, or an off-by-one",
117 )
120class SortState:
121 """Tracks the ordering rule that actually matters: no interleaved contigs,
122 ascending positions within each one."""
124 def __init__(self) -> None:
125 self.seen: set[str] = set()
126 self.current: str | None = None
127 self.position = -1
129 def check(self, find: Findings, lineno: int, contig: str, pos: int, why: str) -> None:
130 if contig != self.current:
131 if contig in self.seen:
132 find.add(
133 "warn", "interleaved_contigs", lineno,
134 f"{contig} reappears after {self.current}; records for one contig "
135 f"must be contiguous. {why}",
136 )
137 self.seen.add(contig)
138 self.current = contig
139 self.position = pos
140 return
141 if pos < self.position:
142 find.add(
143 "warn", "unsorted", lineno,
144 f"{contig}:{pos} follows {contig}:{self.position}. {why}",
145 )
146 self.position = pos
149# --------------------------------------------------------------------------
150# per-format checks
151# --------------------------------------------------------------------------
154def audit_bed(path: Path, find: Findings, genome: dict[str, int], flavour: str = "bed") -> None:
155 contigs: list[str] = []
156 widths: Counter = Counter()
157 total = zero_length = starts_at_zero = 0
158 order = SortState()
160 for lineno, raw in enumerate(path.read_text().splitlines(), start=1):
161 if not raw.strip() or raw.startswith(("#", "track ", "browser ")):
162 continue
163 fields = raw.split("\t")
164 widths[len(fields)] += 1
165 if len(fields) < 3:
166 if len(raw.split()) >= 3:
167 find.add(
168 "fatal", "space_separated", lineno,
169 "fields are separated by spaces, not tabs. BED is tab-delimited; "
170 "a contig name containing a space would be unparseable, and some "
171 "readers take the whole line as one field",
172 )
173 else:
174 find.add(
175 "fatal", "too_few_columns", lineno,
176 "BED needs at least 3 columns",
177 )
178 continue
179 contig = fields[0]
180 contigs.append(contig)
181 try:
182 start, end = int(fields[1]), int(fields[2])
183 except ValueError:
184 find.add(
185 "fatal", "non_integer_coordinate", lineno,
186 f"chromStart/chromEnd are {fields[1]!r}/{fields[2]!r}",
187 )
188 continue
189 total += 1
191 if start < 0:
192 find.add("fatal", "negative_start", lineno, f"chromStart is {start}")
193 if end < start:
194 find.add(
195 "fatal", "end_before_start", lineno,
196 f"chromEnd {end} < chromStart {start}",
197 )
198 if start == 0:
199 starts_at_zero += 1
200 if end == start:
201 zero_length += 1
203 if len(fields) >= 6 and fields[5] not in {"+", "-", "."}:
204 find.add(
205 "fatal", "bad_strand", lineno,
206 f"strand column is {fields[5]!r}; BED allows only +, - or .",
207 )
208 if flavour == "narrowpeak":
209 try:
210 summit = int(fields[9])
211 if summit != -1 and not 0 <= summit < end - start:
212 find.add("fatal", "bad_peak_offset", lineno, "peak must be -1 or an offset inside the interval")
213 except (IndexError, ValueError):
214 find.add("fatal", "bad_peak_offset", lineno, "narrowPeak needs an integer tenth column")
215 if flavour == "bed" and len(fields) >= 8:
216 try:
217 thick_start, thick_end = int(fields[6]), int(fields[7])
218 if not start <= thick_start <= thick_end <= end:
219 find.add(
220 "warn", "thick_outside_feature", lineno,
221 f"thickStart/thickEnd {thick_start}-{thick_end} fall outside "
222 f"the feature {start}-{end}",
223 )
224 except ValueError:
225 find.add("fatal", "non_integer_coordinate", lineno, "thickStart/thickEnd")
226 if flavour == "bed" and len(fields) >= 12:
227 audit_bed12(find, lineno, fields, start, end)
229 order.check(
230 find, lineno, contig, start,
231 "tabix and bedtools -sorted require coordinate-sorted input",
232 )
233 check_bounds(find, genome, lineno, contig, end)
235 if len(widths) > 1:
236 find.add(
237 "warn", "ragged_columns", "file",
238 "column count varies across the file ("
239 + ", ".join(f"{n} cols x{c}" for n, c in sorted(widths.items()))
240 + "); BED is positional, so parsers will misread the narrower rows",
241 )
242 if total and zero_length / total > 0.1:
243 find.add(
244 "warn", "many_zero_length", "file",
245 f"{zero_length}/{total} features have chromEnd == chromStart. In BED "
246 "that is a zero-length insertion point, not a single base. Single-base "
247 "features need chromEnd = chromStart + 1 -- this looks like 1-based "
248 "data written into a 0-based format",
249 )
250 if total >= 100 and starts_at_zero == 0:
251 find.add(
252 "info", "no_zero_start", "file",
253 f"no feature starts at 0 across {total} features. Weak on its own, but "
254 "consistent with 1-based coordinates; check one feature against the "
255 "reference sequence before trusting the file",
256 )
257 audit_naming(find, contigs)
260def audit_bed12(find: Findings, lineno: int, fields: list[str], start: int, end: int) -> None:
261 """The BED12 block rules are exact and routinely broken by hand-written files."""
262 try:
263 count = int(fields[9])
264 sizes = [int(v) for v in fields[10].rstrip(",").split(",") if v != ""]
265 offsets = [int(v) for v in fields[11].rstrip(",").split(",") if v != ""]
266 except (ValueError, IndexError):
267 find.add("fatal", "bad_block_fields", lineno, "blockCount/Sizes/Starts unparseable")
268 return
269 if not (count == len(sizes) == len(offsets)):
270 find.add(
271 "fatal", "block_count_mismatch", lineno,
272 f"blockCount {count} but {len(sizes)} sizes and {len(offsets)} starts",
273 )
274 return
275 if count < 1 or any(s <= 0 for s in sizes) or any(o < 0 for o in offsets):
276 find.add("fatal", "bad_block_geometry", lineno, "blocks need positive sizes/count and nonnegative offsets")
277 if any(o + s > end - start for o, s in zip(offsets, sizes)) or any(
278 a + s > b for a, s, b in zip(offsets, sizes, offsets[1:])
279 ):
280 find.add("fatal", "bad_block_geometry", lineno, "blocks must be ordered, nonoverlapping, and inside the feature")
281 if offsets and offsets[0] != 0:
282 find.add(
283 "fatal", "first_block_offset", lineno,
284 f"blockStarts[0] is {offsets[0]}; it must be 0, because block starts are "
285 "offsets from chromStart, not absolute coordinates",
286 )
287 if offsets and start + offsets[-1] + sizes[-1] != end:
288 find.add(
289 "fatal", "last_block_end", lineno,
290 f"chromStart + last blockStart + last blockSize = "
291 f"{start + offsets[-1] + sizes[-1]}, but chromEnd is {end}; the spec "
292 "requires them to be equal",
293 )
296def audit_gff(path: Path, find: Findings, genome: dict[str, int], flavour: str) -> None:
297 contigs: list[str] = []
298 gtf_attrs = gff3_attrs = 0
299 circular: set[str] = set()
301 for lineno, raw in enumerate(path.read_text().splitlines(), start=1):
302 if raw == "##FASTA":
303 break
304 if not raw.strip() or raw.startswith("#"):
305 continue
306 fields = raw.split("\t")
307 if len(fields) != 9:
308 find.add(
309 "fatal", "wrong_column_count", lineno,
310 f"{len(fields)} columns; GFF/GTF is exactly 9, tab separated",
311 )
312 continue
313 contig, _, feature, start_s, end_s, _, strand, phase, attrs = fields
314 if flavour == "gff3" and "Is_circular=true" in attrs.split(";"):
315 circular.add(canonical_contig(contig))
316 contigs.append(contig)
317 try:
318 start, end = int(start_s), int(end_s)
319 except ValueError:
320 find.add("fatal", "non_integer_coordinate", lineno, f"{start_s!r}/{end_s!r}")
321 continue
323 if start < 1:
324 find.add(
325 "fatal", "start_below_one", lineno,
326 f"start is {start}. GFF/GTF is 1-based; this may be a convention "
327 "error. Verify the source before shifting any records",
328 )
329 if end < start:
330 find.add(
331 "fatal", "end_before_start", lineno,
332 f"end {end} < start {start}. GFF/GTF always writes start <= end, "
333 "including on the minus strand; strand lives in column 7",
334 )
335 if strand not in {"+", "-", ".", "?"}:
336 find.add("fatal", "bad_strand", lineno, f"strand column is {strand!r}")
337 if phase not in {"0", "1", "2", "."}:
338 find.add(
339 "fatal", "bad_phase", lineno,
340 f"phase column is {phase!r}; allowed values are 0, 1, 2 and .",
341 )
342 if feature == "CDS" and phase == ".":
343 find.add(
344 "warn", "cds_without_phase", lineno,
345 "CDS features must declare phase; without it the reading frame is "
346 "undefined and translation will be wrong",
347 )
348 if "=" in attrs and '"' not in attrs:
349 gff3_attrs += 1
350 elif '"' in attrs:
351 gtf_attrs += 1
352 if canonical_contig(contig) in circular and genome.get(canonical_contig(contig), end) < end:
353 find.add("warn", "circular_extent_unchecked", lineno,
354 "circular GFF3 can extend beyond the landmark length; unwrap and validate feature semantics separately")
355 else:
356 check_bounds(find, genome, lineno, contig, end)
358 if flavour == "gtf" and gff3_attrs > gtf_attrs:
359 find.add(
360 "warn", "attribute_syntax", "file",
361 'attributes use GFF3 syntax (key=value;) but the file is named .gtf, '
362 'which expects key "value";. Parsers keyed on the extension will read '
363 "no attributes at all",
364 )
365 if flavour == "gff3" and gtf_attrs > gff3_attrs:
366 find.add(
367 "warn", "attribute_syntax", "file",
368 'attributes use GTF syntax (key "value";) but the file is named .gff3',
369 )
370 audit_naming(find, contigs)
373def audit_vcf(path: Path, find: Findings, genome: dict[str, int]) -> None:
374 contigs: list[str] = []
375 header_seen = False
376 order = SortState()
377 declared: dict[str, int] = {}
379 for lineno, raw in enumerate(path.read_text().splitlines(), start=1):
380 if raw.startswith("##contig="):
381 ident = re.search(r"ID=([^,>]+)", raw)
382 length = re.search(r"length=(\d+)", raw)
383 if ident and length:
384 declared[canonical_contig(ident.group(1))] = int(length.group(1))
385 continue
386 if raw.startswith("#CHROM"):
387 header_seen = True
388 continue
389 if raw.startswith("#") or not raw.strip():
390 continue
392 fields = raw.split("\t")
393 if len(fields) < 8:
394 find.add(
395 "fatal", "too_few_columns", lineno,
396 "VCF data lines need eight fixed fields: CHROM POS ID REF ALT QUAL FILTER INFO",
397 )
398 continue
399 contig, pos_s, _, ref, alt = fields[:5]
400 contigs.append(contig)
401 if not pos_s.isdigit():
402 find.add("fatal", "non_integer_position", lineno, f"POS is {pos_s!r}")
403 continue
404 pos = int(pos_s)
406 nonliteral = any(special_alt(one) for one in alt.split(","))
407 if pos < 1 and not nonliteral:
408 find.add(
409 "fatal", "pos_below_one", lineno,
410 f"POS is {pos}. Ordinary literal variants require POS >= 1; "
411 "telomeric breakends need separate structural validation",
412 )
413 if not ref or set(ref) - VALID_BASES:
414 find.add(
415 "fatal", "bad_ref_allele", lineno,
416 f"REF is {ref!r}; it must be a non-empty ACGTN string. A '-' or an "
417 "empty field is Ensembl/VEP notation, not VCF -- VCF represents "
418 "indels with an anchor base shared by REF and ALT",
419 )
420 for one in alt.split(","):
421 if not special_alt(one) and (not one or set(one) - VALID_BASES):
422 find.add(
423 "fatal", "bad_alt_allele", lineno,
424 f"ALT {one!r} is not a literal ACGTN allele or recognized special allele",
425 )
426 elif one == ref:
427 find.add("warn", "ref_equals_alt", lineno, f"REF and ALT are both {ref}")
428 elif (
429 not special_alt(one)
430 and len(ref) > 1
431 and len(one) > 1
432 and (ref[-1].upper() == one[-1].upper() or ref[0].upper() == one[0].upper())
433 ):
434 find.add(
435 "warn", "not_parsimonious", lineno,
436 f"{ref}>{one} share a removable flanking base, so this record is not "
437 "trimmed. Run normalize_variant.py before comparing or joining "
438 "on these alleles",
439 )
440 if "," in alt:
441 find.add(
442 "info", "multiallelic", lineno,
443 "multi-allelic record; split it before normalising or joining",
444 )
445 order.check(
446 find, lineno, contig, pos,
447 "tabix indexing requires coordinate-sorted input",
448 )
450 far = pos + max(len(ref) - 1, 0)
451 if nonliteral:
452 find.add("warn", "structural_extent_unchecked", lineno,
453 "nonliteral ALT: SV/gVCF extent, breakend syntax, and telomere bounds are not validated")
454 else:
455 check_bounds(
456 find,
457 genome or declared,
458 lineno,
459 contig,
460 far,
461 "the genome file" if genome else "the file's own ##contig headers",
462 )
464 if not header_seen:
465 find.add(
466 "warn", "missing_chrom_header", "file",
467 "no #CHROM line; the file is not a valid VCF and column meanings are "
468 "being guessed",
469 )
470 audit_naming(find, contigs)
473def audit_naming(find: Findings, contigs: list[str]) -> None:
474 style = naming_style(sorted(set(contigs)))
475 if style == "mixed":
476 prefixed = sorted({c for c in contigs if c.startswith("chr")})[:3]
477 plain = sorted({c for c in contigs if not c.startswith("chr")})[:3]
478 find.add(
479 "warn", "mixed_contig_naming", "file",
480 f"the file mixes chr-prefixed ({', '.join(prefixed)}) and plain "
481 f"({', '.join(plain)}) contig names; this can be valid for custom references. "
482 "Compare exact names with the intended sequence dictionary",
483 )
486def build_parser() -> argparse.ArgumentParser:
487 parser = argparse.ArgumentParser(
488 description="Audit a BED/GTF/GFF3/VCF file against its format's conventions."
489 )
490 parser.add_argument("file", type=Path)
491 parser.add_argument("--format", choices=("bed", "bedgraph", "narrowpeak", "broadpeak", "gtf", "gff3", "vcf"))
492 parser.add_argument(
493 "--genome", type=Path, help=".chrom.sizes or .fai to bounds-check against"
494 )
495 parser.add_argument("--max-examples", type=int, default=5)
496 parser.add_argument("--output-format", choices=("tsv", "json"), default="tsv")
497 parser.add_argument("-o", "--output")
498 return parser
501def main(argv: list[str] | None = None) -> int:
502 parser = build_parser()
503 args = parser.parse_args(argv)
504 if not args.file.exists():
505 print(f"error: no such file: {args.file}", file=sys.stderr)
506 return 2
507 if args.genome and not args.genome.exists():
508 print(f"error: no such file: {args.genome}", file=sys.stderr)
509 return 2
511 fmt = args.format or detect_format(args.file)
512 genome = load_genome(args.genome)
513 find = Findings(args.max_examples)
515 if fmt in {"bed", "bedgraph", "narrowpeak", "broadpeak"}:
516 audit_bed(args.file, find, genome, fmt)
517 elif fmt in {"gtf", "gff3"}:
518 audit_gff(args.file, find, genome, fmt)
519 else:
520 audit_vcf(args.file, find, genome)
522 if find.rows:
523 emit(find.rows, COLUMNS, args.output_format, args.output)
524 else:
525 print(f"{args.file.name}: no findings ({fmt} conventions)")
526 sys.stdout.flush()
528 if find.rows:
529 summary = ", ".join(
530 f"{count} {severity}"
531 for severity in ("fatal", "warn", "info")
532 if (count := find.by_severity[severity])
533 )
534 print(f"\n{args.file.name}: {summary}", file=sys.stderr)
535 return 1 if find.fatal else 0
538if __name__ == "__main__":
539 raise SystemExit(main())