Skills › Scientific Agent Skills › genomic-coordinates › check_contigs.py
check_contigs.py
400 lines · skills/genomic-coordinates/scripts/check_contigs.py at 5f2bc03
1#!/usr/bin/env python3
2"""Identify the assembly behind a file, and check that two files can be joined.
4Contig-name and declared-length conflicts can reveal incompatible inputs.
5Matching lengths cannot establish sequence identity or rule out masking.
7 python3 check_contigs.py --identify ref.fa.fai
8 python3 check_contigs.py variants.vcf annotation.gtf ref.fa.fai
9 python3 check_contigs.py peaks.bed --genome hg38.chrom.sizes
11Exit codes: 0 no detected conflict (possibly unknown), 1 conflict, 2 usage error.
12"""
14from __future__ import annotations
16import argparse
17import re
18import sys
19from pathlib import Path
21sys.path.insert(0, str(Path(__file__).resolve().parent))
23from _common import canonical_contig, emit, naming_style # noqa: E402
25# Primary-chromosome lengths, read from the UCSC bigZips chrom.sizes for each
26# assembly (hg19, hg38, hs1) and cross-checked against the NCBI assembly report
27# for GRCh37.p13. UCSC tables rechecked 2026-10-01.
28#
29# GRCh37 and hg19 are the same assembly for every contig here except chrM: UCSC
30# kept the older NC_001807 mitochondrion (16,571 bp) while GRCh37 adopted the
31# rCRS (16,569 bp). That two-base difference is the only signal in the primary
32# chromosomes that tells the two apart, and it is why an hg19 BAM and a GRCh37
33# VCF can disagree about mitochondrial variants.
34BUILDS: dict[str, dict[str, int]] = {
35 "GRCh37/hg19": {
36 "1": 249250621, "2": 243199373, "3": 198022430, "4": 191154276,
37 "5": 180915260, "6": 171115067, "7": 159138663, "8": 146364022,
38 "9": 141213431, "10": 135534747, "11": 135006516, "12": 133851895,
39 "13": 115169878, "14": 107349540, "15": 102531392, "16": 90354753,
40 "17": 81195210, "18": 78077248, "19": 59128983, "20": 63025520,
41 "21": 48129895, "22": 51304566, "X": 155270560, "Y": 59373566,
42 },
43 "GRCh38/hg38": {
44 "1": 248956422, "2": 242193529, "3": 198295559, "4": 190214555,
45 "5": 181538259, "6": 170805979, "7": 159345973, "8": 145138636,
46 "9": 138394717, "10": 133797422, "11": 135086622, "12": 133275309,
47 "13": 114364328, "14": 107043718, "15": 101991189, "16": 90338345,
48 "17": 83257441, "18": 80373285, "19": 58617616, "20": 64444167,
49 "21": 46709983, "22": 50818468, "X": 156040895, "Y": 57227415,
50 },
51 "T2T-CHM13v2.0/hs1": {
52 "1": 248387328, "2": 242696752, "3": 201105948, "4": 193574945,
53 "5": 182045439, "6": 172126628, "7": 160567428, "8": 146259331,
54 "9": 150617247, "10": 134758134, "11": 135127769, "12": 133324548,
55 "13": 113566686, "14": 101161492, "15": 99753195, "16": 96330374,
56 "17": 84276897, "18": 80542538, "19": 61707364, "20": 66210255,
57 "21": 45090682, "22": 51324926, "X": 154259566, "Y": 62460029,
58 },
59}
61MITO = {16571: "hg19 (UCSC NC_001807 chrM)", 16569: "GRCh37/38 (rCRS MT)"}
63COLUMNS = ["file", "kind", "contigs", "naming", "assembly", "detail"]
66# --------------------------------------------------------------------------
67# readers
68# --------------------------------------------------------------------------
71def read_sizes(path: Path) -> tuple[dict[str, int], str]:
72 """.fai, .chrom.sizes, or .genome -- name in column 1, length in column 2."""
73 sizes = {}
74 for line in path.read_text().splitlines():
75 if not line.strip() or line.startswith("#"):
76 continue
77 parts = line.split("\t") if "\t" in line else line.split()
78 if len(parts) >= 2 and parts[1].isdigit():
79 sizes[parts[0]] = int(parts[1])
80 return sizes, "sizes"
83def read_fasta_names(path: Path) -> tuple[dict[str, int], str]:
84 sizes: dict[str, int] = {}
85 name, count = None, 0
86 with path.open() as handle:
87 for line in handle:
88 if line.startswith(">"):
89 if name is not None:
90 sizes[name] = count
91 name, count = line[1:].split()[0], 0
92 elif name is not None:
93 count += len(line.strip())
94 if name is not None:
95 sizes[name] = count
96 return sizes, "fasta"
99def read_vcf(path: Path) -> tuple[dict[str, int], str]:
100 """Prefer the ##contig header; fall back to the CHROM column."""
101 header: dict[str, int] = {}
102 seen: dict[str, int] = {}
103 with path.open() as handle:
104 for line in handle:
105 if line.startswith("##contig="):
106 ident = re.search(r"ID=([^,>]+)", line)
107 length = re.search(r"length=(\d+)", line)
108 if ident:
109 header[ident.group(1)] = int(length.group(1)) if length else 0
110 elif not line.startswith("#") and line.strip():
111 fields = line.split("\t")
112 if len(fields) >= 2 and fields[1].isdigit():
113 pos = int(fields[1]) + (len(fields[3]) - 1 if len(fields) >= 4 else 0)
114 seen[fields[0]] = max(seen.get(fields[0], 0), pos)
115 if header:
116 return header, "vcf header"
117 return seen, "vcf records (max REF end, not contig length)"
120def read_vcf_record_extents(path: Path) -> dict[str, int]:
121 """Check records even when a header exists; header lengths are not observations."""
122 seen: dict[str, int] = {}
123 with path.open() as handle:
124 for line in handle:
125 if line.startswith("#") or not line.strip():
126 continue
127 fields = line.rstrip().split("\t")
128 if len(fields) >= 4 and fields[1].isdigit():
129 far = int(fields[1]) + len(fields[3]) - 1
130 seen[fields[0]] = max(seen.get(fields[0], 0), far)
131 return seen
134def read_sam_header(path: Path) -> tuple[dict[str, int], str]:
135 sizes: dict[str, int] = {}
136 with path.open() as handle:
137 for line in handle:
138 if not line.startswith("@"):
139 break
140 if line.startswith("@SQ"):
141 ident = re.search(r"SN:(\S+)", line)
142 length = re.search(r"LN:(\d+)", line)
143 if ident and length:
144 sizes[ident.group(1)] = int(length.group(1))
145 return sizes, "sam header"
148def read_intervals(path: Path, start_col: int, end_col: int) -> tuple[dict[str, int], str]:
149 """BED/GTF/GFF: names from column 1, plus the largest coordinate observed."""
150 sizes: dict[str, int] = {}
151 with path.open() as handle:
152 for line in handle:
153 if line.rstrip() == "##FASTA":
154 break
155 if not line.strip() or line.startswith(("#", "track", "browser")):
156 continue
157 fields = line.rstrip("\n").split("\t")
158 if len(fields) <= end_col:
159 continue
160 try:
161 far = max(int(fields[start_col]), int(fields[end_col]))
162 except ValueError:
163 continue
164 sizes[fields[0]] = max(sizes.get(fields[0], 0), far)
165 return sizes, "intervals (max coordinate, not contig length)"
168def load(path: Path) -> tuple[dict[str, int], str]:
169 suffixes = [s.lower() for s in path.suffixes]
170 name = path.name.lower()
171 if name.endswith((".fai", ".chrom.sizes", ".sizes", ".genome")):
172 return read_sizes(path)
173 if ".vcf" in suffixes:
174 return read_vcf(path)
175 if ".sam" in suffixes or name.endswith(".header") or name.endswith(".header.txt"):
176 return read_sam_header(path)
177 if ".gtf" in suffixes or ".gff" in suffixes or ".gff3" in suffixes:
178 return read_intervals(path, 3, 4)
179 if ".bed" in suffixes or ".narrowpeak" in suffixes or ".broadpeak" in suffixes:
180 return read_intervals(path, 1, 2)
181 if ".fa" in suffixes or ".fasta" in suffixes or ".fna" in suffixes:
182 return read_fasta_names(path)
183 raise SystemExit(
184 f"error: cannot tell what kind of file {path.name} is. Supported: "
185 ".fai, .chrom.sizes, .vcf, .sam (header), .bed, .gtf, .gff, .fasta"
186 )
189# --------------------------------------------------------------------------
190# identification
191# --------------------------------------------------------------------------
194def identify(sizes: dict[str, int], exact: bool) -> tuple[str, str]:
195 """Match observed contig lengths against the known assemblies.
197 ``exact`` is False for interval files, where the numbers are the largest
198 coordinate seen rather than the contig length -- there the useful question is
199 only whether anything overflows a candidate assembly.
200 """
201 folded = {canonical_contig(n): length for n, length in sizes.items()}
202 if len(folded) != len(sizes):
203 return "CONFLICT", "multiple contig names collapse to the same comparison alias"
204 if exact:
205 folded = {name: length for name, length in folded.items() if length > 0}
206 mito = folded.get("M")
208 scores = []
209 for build, lengths in BUILDS.items():
210 shared = [c for c in lengths if c in folded]
211 if not shared:
212 continue
213 if exact:
214 agree = sum(1 for c in shared if folded[c] == lengths[c])
215 else:
216 agree = sum(1 for c in shared if 0 < folded[c] <= lengths[c])
217 scores.append((agree / len(shared), agree, len(shared), build))
218 if not scores:
219 return "unknown", "no primary chromosomes to match against"
220 scores.sort(reverse=True)
221 fraction, agree, total, build = scores[0]
223 if not exact:
224 fits = [b for f, _, _, b in scores if f == 1.0]
225 if len(fits) == 1:
226 return f"consistent with {fits[0]}", "no coordinate overflows this assembly"
227 if not fits:
228 over = [c for c in folded if c in BUILDS[build] and folded[c] > BUILDS[build][c]]
229 return "CONFLICT", (
230 f"coordinates run past the end of every known assembly (closest is "
231 f"{build}, overflowing on {', '.join(sorted(over)[:4])}); wrong "
232 "assembly, or 1-based data written into a 0-based file"
233 )
234 return "ambiguous", f"coordinates fit any of: {', '.join(fits)}"
236 if fraction < 0.9:
237 return "unknown", (
238 f"closest is {build} at {agree}/{total} chromosomes -- not a match"
239 )
240 detail = f"{agree}/{total} primary chromosome lengths match (length evidence only, not sequence identity)"
241 if build == "GRCh37/hg19" and mito is not None:
242 detail += f"; chrM is {mito} bp, i.e. {MITO.get(mito, 'a non-standard mitochondrion')}"
243 build = "hg19" if mito == 16571 else "GRCh37" if mito == 16569 else build
244 elif mito is not None and mito not in MITO:
245 detail += f"; chrM length {mito} is not a standard human mitochondrion"
246 if agree != total:
247 wrong = [c for c in BUILDS[scores[0][3]] if c in folded and folded[c] != BUILDS[scores[0][3]][c]]
248 detail += f"; disagrees on {', '.join(sorted(wrong))}"
249 return "CONFLICT", f"closest is {build}; {detail}"
250 return build, detail
253# --------------------------------------------------------------------------
254# comparison
255# --------------------------------------------------------------------------
258def compare(loaded: list[tuple[Path, dict[str, int], str, bool]]) -> list[str]:
259 """Report detectable name, length, and observed-bound conflicts.
261 Each entry carries ``exact``: True when its numbers are declared contig
262 lengths, False when they are only the largest coordinate observed. Two exact
263 files must agree exactly; an inexact one can only ever be caught overflowing.
264 """
265 problems: list[str] = []
266 for i in range(len(loaded)):
267 for j in range(i + 1, len(loaded)):
268 path_a, sizes_a, _, exact_a = loaded[i]
269 path_b, sizes_b, _, exact_b = loaded[j]
270 fold_a = {canonical_contig(n): v for n, v in sizes_a.items()}
271 fold_b = {canonical_contig(n): v for n, v in sizes_b.items()}
272 shared = set(fold_a) & set(fold_b)
273 aliases_a = {canonical_contig(n): n for n in sizes_a}
274 aliases_b = {canonical_contig(n): n for n in sizes_b}
275 renamed = sorted(c for c in shared if aliases_a[c] != aliases_b[c])
276 if renamed:
277 problems.append(f"contig naming differs: {path_a.name} and {path_b.name} use different exact names for {', '.join(renamed[:5])}; verify sequence identity before renaming")
278 if not shared:
279 problems.append(
280 f"{path_a.name} and {path_b.name} share no contigs at all, "
281 "even ignoring the chr prefix"
282 )
283 continue
285 if exact_a and exact_b:
286 clash = sorted(
287 c for c in shared
288 if fold_a[c] and fold_b[c] and fold_a[c] != fold_b[c]
289 )
290 if clash:
291 example = clash[0]
292 problems.append(
293 f"{path_a.name} and {path_b.name} disagree on contig length "
294 f"for {', '.join(clash[:5])}"
295 f"{'...' if len(clash) > 5 else ''} "
296 f"(e.g. {example}: {fold_a[example]} vs {fold_b[example]}) -- "
297 "these are different assemblies and their coordinates are "
298 "not comparable"
299 )
300 else:
301 # One side holds only the largest coordinate observed, so equality
302 # proves nothing. Overflowing a declared contig length still does.
303 pairs = []
304 if exact_b and not exact_a:
305 pairs.append((path_a.name, fold_a, path_b.name, fold_b))
306 if exact_a and not exact_b:
307 pairs.append((path_b.name, fold_b, path_a.name, fold_a))
308 for name_hi, hi, name_lo, lo in pairs:
309 over = sorted(c for c in shared if lo[c] and hi[c] > lo[c])
310 if over:
311 example = over[0]
312 problems.append(
313 f"{name_hi} has coordinates past the end of {name_lo} on "
314 f"{', '.join(over[:5])}{'...' if len(over) > 5 else ''} "
315 f"(e.g. {example}: {hi[example]} > {lo[example]}) -- "
316 "wrong assembly, or an off-by-one from a 1-based source"
317 )
319 for source, observed, source_exact, target, declared, target_exact in (
320 (path_a, fold_a, exact_a, path_b, fold_b, exact_b),
321 (path_b, fold_b, exact_b, path_a, fold_a, exact_a),
322 ):
323 missing = sorted(set(observed) - set(declared))
324 # Reference supersets are normal. Observed record contigs must
325 # exist in a supplied reference/header, in either file order.
326 if target_exact and not source_exact and missing:
327 problems.append(f"{source.name} has contigs absent from {target.name}: {', '.join(missing[:5])}")
328 return problems
331def build_parser() -> argparse.ArgumentParser:
332 parser = argparse.ArgumentParser(
333 description="Identify assemblies and check contig compatibility between files."
334 )
335 parser.add_argument("files", nargs="+", type=Path)
336 parser.add_argument(
337 "--identify",
338 action="store_true",
339 help="only report the assembly of each file, skipping the cross-file checks",
340 )
341 parser.add_argument(
342 "--genome",
343 type=Path,
344 help="a .chrom.sizes/.fai to check every other file against",
345 )
346 parser.add_argument("--format", choices=("tsv", "json"), default="tsv")
347 parser.add_argument("-o", "--output")
348 return parser
351def main(argv: list[str] | None = None) -> int:
352 parser = build_parser()
353 args = parser.parse_args(argv)
355 paths = list(args.files)
356 if args.genome:
357 paths.append(args.genome)
358 for path in paths:
359 if not path.exists():
360 print(f"error: no such file: {path}", file=sys.stderr)
361 return 2
363 loaded = []
364 rows = []
365 for path in paths:
366 sizes, kind = load(path)
367 exact = kind in {"sizes", "fasta", "vcf header", "sam header"}
368 loaded.append((path, sizes, kind, exact))
369 if kind == "vcf header":
370 loaded.append((path, read_vcf_record_extents(path), "vcf record extents", False))
371 assembly, detail = identify(sizes, exact)
372 rows.append(
373 {
374 "file": path.name,
375 "kind": kind,
376 "contigs": len(sizes),
377 "naming": naming_style(list(sizes)),
378 "assembly": assembly,
379 "detail": detail,
380 }
381 )
383 emit(rows, COLUMNS, args.format, args.output)
384 sys.stdout.flush()
386 failed = any(r["assembly"] == "CONFLICT" for r in rows)
387 if not args.identify and len(loaded) > 1:
388 problems = compare(loaded)
389 if problems:
390 failed = True
391 print("\nincompatibilities:", file=sys.stderr)
392 for problem in problems:
393 print(f" - {problem}", file=sys.stderr)
394 else:
395 print("\nno contig incompatibilities found", file=sys.stderr)
396 return 1 if failed else 0
399if __name__ == "__main__":
400 raise SystemExit(main())