Skills › Scientific Agent Skills › genomic-coordinates › convert_coords.py
convert_coords.py
189 lines · skills/genomic-coordinates/scripts/convert_coords.py at 5f2bc03
1#!/usr/bin/env python3
2"""Convert intervals between genomic coordinate conventions.
4Every conversion goes through one canonical form (0-based half-open), so the
5answer never depends on remembering which pair of formats is involved.
7 python3 convert_coords.py --from bed --to gff chr1 999 1000
8 python3 convert_coords.py --from ucsc --to bed "chr7:5,530,601-5,530,625"
9 python3 convert_coords.py --from bed --to ensembl --input peaks.bed
10 python3 convert_coords.py --list
12Exit codes: 0 fine, 1 at least one interval is degenerate or invalid,
132 usage error.
14"""
16from __future__ import annotations
18import argparse
19import sys
20from pathlib import Path
22sys.path.insert(0, str(Path(__file__).resolve().parent))
24from _common import ( # noqa: E402
25 CONVENTIONS,
26 REGION_STRING_FORMATS,
27 Convention,
28 emit,
29 format_region,
30 from_canonical,
31 get_convention,
32 iter_data_lines,
33 parse_region,
34 to_canonical,
35)
37COLUMNS = ["contig", "input", "output", "length", "status", "detail"]
40def convert_one(
41 contig: str,
42 start: int,
43 end: int,
44 src: Convention,
45 dst: Convention,
46 strand: str | None = None,
47) -> dict:
48 """Convert one interval and describe anything the caller should look at."""
49 start0, end0 = to_canonical(src, start, end)
50 length = end0 - start0
52 status, detail = "ok", ""
53 if start0 < 0:
54 status, detail = "invalid", f"start {start} is below the origin of {src.label}"
55 elif length < 0 or (not src.half_open and end < start):
56 status = "invalid"
57 detail = (
58 f"end precedes start in {src.label} coordinates; "
59 "an interval file with start > end is corrupt, not merely mis-converted"
60 )
61 elif length == 0:
62 if dst.half_open:
63 status = "zero_length"
64 detail = "zero-length interval (an insertion point between two bases)"
65 else:
66 status = "unrepresentable"
67 detail = (
68 f"zero-length interval needs format-specific feature semantics in {dst.label}; "
69 f"the arithmetic gives end = start - 1. Keep it in a half-open "
70 "format, or record the insertion against its anchor base as VCF does"
71 )
73 out_start, out_end = from_canonical(dst, start0, end0)
75 def render(conv: Convention, s: int, e: int) -> str:
76 if conv.name in REGION_STRING_FORMATS:
77 return format_region(contig, s, e, strand, conv.name)
78 return f"{s}-{e}"
80 return {
81 "contig": contig,
82 "input": render(src, start, end),
83 "output": render(dst, out_start, out_end) if status not in {"invalid", "unrepresentable"} else "",
84 "length": length,
85 "status": status,
86 "detail": detail,
87 "start": out_start,
88 "end": out_end,
89 "from": src.name,
90 "to": dst.name,
91 }
94def read_intervals(path: str, conv: Convention) -> list[tuple[str, int, int, str | None]]:
95 """Read intervals from a file: columnar for file formats, one region per line
96 for the region-string formats."""
97 out: list[tuple[str, int, int, str | None]] = []
98 supported = REGION_STRING_FORMATS | {"bed", "bedgraph", "gff", "gff3", "gtf", "vcf", "granges", "pyranges", "python"}
99 if conv.name not in supported:
100 raise SystemExit(f"--input does not parse {conv.name} files; extract contig/start/end and pass an explicit triple")
101 for lineno, line in iter_data_lines(path):
102 try:
103 if conv.name in REGION_STRING_FORMATS:
104 contig, start, end, strand = parse_region(line.split("\t")[0], conv)
105 else:
106 fields = line.split("\t")
107 if len(fields) < 3:
108 raise ValueError("need at least contig, start, end columns")
109 if conv.name == "gff" or conv.name == "gff3" or conv.name == "gtf":
110 contig, start, end = fields[0], int(fields[3]), int(fields[4])
111 strand = fields[6] if len(fields) > 6 else None
112 elif conv.name == "vcf":
113 contig, pos, ref = fields[0], int(fields[1]), fields[3]
114 if len(fields) < 8 or not ref or set(ref.upper()) - set("ACGTN"):
115 raise ValueError("expected a VCF record with eight fixed fields and literal REF")
116 if any(not a or set(a.upper()) - set("ACGTN") for a in fields[4].split(",")):
117 raise ValueError("only literal VCF alleles have a simple REF-span conversion; handle SV/gVCF/breakends separately")
118 start, end, strand = pos, pos + len(ref) - 1, None
119 else:
120 contig, start, end = fields[0], int(fields[1]), int(fields[2])
121 strand = fields[5] if len(fields) > 5 and fields[5] in "+-." else None
122 except (ValueError, IndexError) as exc:
123 raise SystemExit(f"{path}:{lineno}: {exc}") from exc
124 out.append((contig, start, end, strand))
125 return out
128def build_parser() -> argparse.ArgumentParser:
129 parser = argparse.ArgumentParser(
130 description="Convert intervals between genomic coordinate conventions.",
131 formatter_class=argparse.RawDescriptionHelpFormatter,
132 )
133 parser.add_argument("contig", nargs="?", help="contig name, or a region string")
134 parser.add_argument("start", nargs="?", type=int)
135 parser.add_argument("end", nargs="?", type=int)
136 parser.add_argument("--from", dest="src", help="source convention")
137 parser.add_argument("--to", dest="dst", help="target convention")
138 parser.add_argument("--input", help="BED/bedGraph, GFF/GTF, literal-allele VCF, three-column TSV for granges/pyranges/python, or explicit region strings")
139 parser.add_argument("--format", choices=("tsv", "json"), default="tsv")
140 parser.add_argument("-o", "--output", help="write here instead of stdout")
141 parser.add_argument(
142 "--list", action="store_true", help="print the convention table and exit"
143 )
144 return parser
147def main(argv: list[str] | None = None) -> int:
148 parser = build_parser()
149 args = parser.parse_args(argv)
151 if args.list:
152 width = max(len(name) for name in CONVENTIONS)
153 print(f"{'format'.ljust(width)} convention notes")
154 for name, conv in CONVENTIONS.items():
155 print(f"{name.ljust(width)} {conv.label.ljust(20)} {conv.note}")
156 return 0
158 if not args.src or not args.dst:
159 parser.error("--from and --to are both required")
160 try:
161 src = get_convention(args.src)
162 dst = get_convention(args.dst)
163 except ValueError as exc:
164 parser.error(str(exc))
166 intervals: list[tuple[str, int, int, str | None]] = []
167 if args.input:
168 intervals = read_intervals(args.input, src)
169 elif args.contig is not None:
170 if args.start is None:
171 try:
172 intervals = [parse_region(args.contig, src)]
173 except ValueError as exc:
174 parser.error(str(exc))
175 else:
176 end = args.end
177 if end is None:
178 end = args.start if not src.half_open else args.start + 1
179 intervals = [(args.contig, args.start, end, None)]
180 else:
181 parser.error("give a region, a contig/start/end triple, or --input")
183 rows = [convert_one(c, s, e, src, dst, strand) for c, s, e, strand in intervals]
184 emit(rows, COLUMNS, args.format, args.output)
185 return 1 if any(r["status"] in {"invalid", "unrepresentable"} for r in rows) else 0
188if __name__ == "__main__":
189 raise SystemExit(main())