Skills › Scientific Agent Skills › genomic-coordinates › normalize_variant.py
normalize_variant.py
315 lines · skills/genomic-coordinates/scripts/normalize_variant.py at 5f2bc03
1#!/usr/bin/env python3
2"""Normalise VCF-style variants: check REF, trim, and left-align.
4Two variant records can describe exactly the same change to the genome and share
5no field values at all. Comparing, joining, or deduplicating variants without
6normalising first silently loses real matches. This implements the parsimony +
7left-alignment procedure of Tan, Abecasis & Kang (2015), which is what
8``bcftools norm`` and ``vt normalize`` implement.
10 python3 normalize_variant.py --fasta ref.fa chr1 7 CAC C
11 python3 normalize_variant.py --fasta ref.fa --input variants.vcf
12 python3 normalize_variant.py --fasta ref.fa --compare chr1:7:CAC:C chr1:3:CAC:C
14Exit codes: 0 no literal-allele failure (may include unvalidated nonliteral
15pass-through rows), 1 mismatch, invalid/incomplete normalization, or an
16unverifiable comparison, 2 usage or reference error. Inspect ref_check per row.
17"""
19from __future__ import annotations
21import argparse
22import re
23import sys
24from pathlib import Path
26sys.path.insert(0, str(Path(__file__).resolve().parent))
28from _common import Reference, ReferenceError, emit, iter_data_lines # noqa: E402
30COLUMNS = [
31 "input",
32 "normalized",
33 "type",
34 "pos_shift",
35 "ref_check",
36 "changed",
37 "detail",
38]
40DEFAULT_WINDOW = 1000
41BASES = set("ACGTN")
44def special_alt(alt: str) -> bool:
45 """Recognize nonliteral alleles; pass-through is not structural validation."""
46 return (alt in {"*", "."}
47 or re.fullmatch(r"<[^<>\s,]+>", alt) is not None
48 or "[" in alt or "]" in alt
49 or (len(alt) > 1 and (alt.startswith(".") or alt.endswith("."))))
52class VariantError(ValueError):
53 """The record cannot be normalised as written."""
56def classify(ref: str, alt: str) -> str:
57 if len(ref) == len(alt) == 1:
58 return "snv"
59 if len(ref) == len(alt):
60 return "mnv"
61 if len(ref) == 1 and len(alt) > 1 and (alt.startswith(ref) or alt.endswith(ref)):
62 return "insertion"
63 if len(alt) == 1 and len(ref) > 1 and (ref.startswith(alt) or ref.endswith(alt)):
64 return "deletion"
65 return "complex"
68def normalize(
69 reference: Reference,
70 contig: str,
71 pos: int,
72 ref: str,
73 alt: str,
74 window: int = DEFAULT_WINDOW,
75) -> dict:
76 """Normalise a single bi-allelic record. ``pos`` is 1-based, as in VCF.
78 Returns the normalised record plus how far it moved and whether the stated
79 REF actually matched the reference sequence. ``shifted`` is positive when
80 left-alignment walked the anchor left through a repeat, and negative when
81 trimming redundant flanking bases moved it right onto a parsimonious
82 representation (``3 CA>CT`` is really the SNV ``4 A>T``).
83 """
84 if pos < 1:
85 raise VariantError(f"POS {pos} is not 1-based; VCF positions start at 1")
86 if window < 1:
87 raise VariantError("window must be a positive integer")
88 if not ref or set(ref.upper()) - BASES:
89 raise VariantError("REF must be a non-empty ACGTN sequence")
90 if "," in alt:
91 raise VariantError("multiple ALT alleles: use --split before allele-wise normalization")
92 if special_alt(alt):
93 return {
94 "pos": pos,
95 "ref": ref,
96 "alt": alt,
97 "shifted": 0,
98 "ref_check": "skipped",
99 "type": "symbolic",
100 "detail": "nonliteral ALT: passed through without REF, structure, or equivalence validation",
101 }
102 if not alt or set(alt.upper()) - BASES:
103 raise VariantError("ALT must be a non-empty ACGTN sequence or a supported nonliteral allele")
105 ref, alt = ref.upper(), alt.upper()
106 observed = reference.fetch(contig, pos - 1, pos - 1 + len(ref))
107 if observed != ref:
108 return {
109 "pos": pos,
110 "ref": ref,
111 "alt": alt,
112 "shifted": 0,
113 "ref_check": "MISMATCH",
114 "type": classify(ref, alt),
115 "detail": (
116 f"REF says {ref} but the reference has {observed or '(past contig end)'} "
117 f"at {contig}:{pos}. Do not normalise this record -- the variants and "
118 "the FASTA may differ in assembly, sequence, orientation, or coordinates"
119 ),
120 }
121 if ref == alt:
122 raise VariantError(f"REF and ALT are both {ref}; this record asserts no change")
124 start_pos = pos
125 limit = max(0, pos - 1 - window)
126 incomplete = False
128 # Right-trim and left-extend until the alleles no longer share a final base.
129 while ref[-1] == alt[-1]:
130 if len(ref) == 1 or len(alt) == 1:
131 if pos - 1 <= limit:
132 incomplete = pos > 1
133 break
134 base = reference.fetch(contig, pos - 2, pos - 1)
135 if not base:
136 break
137 if base not in BASES:
138 raise VariantError("left alignment reached an ambiguous reference base outside ACGTN")
139 pos -= 1
140 ref, alt = base + ref, base + alt
141 ref, alt = ref[:-1], alt[:-1]
143 # Left-trim any shared leading bases, keeping one anchor base for indels.
144 while len(ref) > 1 and len(alt) > 1 and ref[0] == alt[0]:
145 ref, alt = ref[1:], alt[1:]
146 pos += 1
148 shifted = start_pos - pos
149 detail = ""
150 if incomplete:
151 detail = (
152 f"left-alignment stopped at the {window} bp window; the repeat may extend "
153 "further. Re-run with a larger --window to confirm"
154 )
155 return {
156 "pos": pos,
157 "ref": ref,
158 "alt": alt,
159 "shifted": shifted,
160 "ref_check": "incomplete" if incomplete else "ok",
161 "type": classify(ref, alt),
162 "detail": detail,
163 }
166def parse_spec(text: str) -> tuple[str, int, str, str]:
167 """Parse ``contig:pos:ref:alt``."""
168 parts = text.rsplit(":", 3)
169 if len(parts) != 4:
170 raise VariantError(f"expected contig:pos:ref:alt, got {text!r}")
171 contig, pos, ref, alt = parts
172 if not pos.isdigit():
173 raise VariantError(f"POS {pos!r} in {text!r} is not a number")
174 return contig, int(pos), ref, alt
177def read_records(path: str, split: bool) -> list[tuple[str, int, str, str]]:
178 """Read CHROM/POS/REF/ALT from a VCF, or from a bare 4-column TSV.
180 Five or more columns are read as VCF (CHROM POS ID REF ALT); exactly four as
181 CHROM POS REF ALT.
182 """
183 records: list[tuple[str, int, str, str]] = []
184 for lineno, line in iter_data_lines(path):
185 fields = line.split("\t")
186 if len(fields) >= 5:
187 contig, pos, ref, alt = fields[0], fields[1], fields[3], fields[4]
188 elif len(fields) == 4:
189 contig, pos, ref, alt = fields[0], fields[1], fields[2], fields[3]
190 else:
191 raise SystemExit(f"{path}:{lineno}: need CHROM, POS, REF, ALT columns")
192 if not pos.isdigit():
193 raise SystemExit(f"{path}:{lineno}: POS {pos!r} is not a number")
194 for one in alt.split(",") if split else [alt]:
195 records.append((contig, int(pos), ref, one))
196 return records
199def build_parser() -> argparse.ArgumentParser:
200 parser = argparse.ArgumentParser(
201 description="Check, trim, and left-align VCF-style variants."
202 )
203 parser.add_argument("contig", nargs="?")
204 parser.add_argument("pos", nargs="?", type=int)
205 parser.add_argument("ref", nargs="?")
206 parser.add_argument("alt", nargs="?")
207 parser.add_argument("--fasta", required=True, help="reference FASTA (.fai used if present)")
208 parser.add_argument("--input", help="VCF, or a TSV of contig/pos/ref/alt")
209 parser.add_argument(
210 "--compare",
211 nargs="+",
212 metavar="CONTIG:POS:REF:ALT",
213 help="normalise two or more records and report whether they are the same variant",
214 )
215 parser.add_argument(
216 "--split",
217 action="store_true",
218 help="split comma-separated ALTs into one record each before normalising",
219 )
220 parser.add_argument("--window", type=int, default=DEFAULT_WINDOW,
221 help="maximum left shift (default 1000); hitting it is incomplete and exits 1")
222 parser.add_argument("--format", choices=("tsv", "json"), default="tsv")
223 parser.add_argument("-o", "--output")
224 return parser
227def run(records, reference, window) -> list[dict]:
228 rows = []
229 for contig, pos, ref, alt in records:
230 label = f"{contig}:{pos}:{ref}:{alt}"
231 try:
232 result = normalize(reference, contig, pos, ref, alt, window)
233 except (VariantError, ReferenceError) as exc:
234 rows.append(
235 {
236 "input": label,
237 "normalized": "",
238 "type": "",
239 "pos_shift": "",
240 "ref_check": "error",
241 "changed": "",
242 "detail": str(exc),
243 }
244 )
245 continue
246 norm = f"{contig}:{result['pos']}:{result['ref']}:{result['alt']}"
247 rows.append(
248 {
249 "input": label,
250 "normalized": norm,
251 "type": result["type"],
252 "pos_shift": result["shifted"],
253 "ref_check": result["ref_check"],
254 "changed": "yes" if norm != label else "no",
255 "detail": result["detail"],
256 }
257 )
258 return rows
261def main(argv: list[str] | None = None) -> int:
262 parser = build_parser()
263 args = parser.parse_args(argv)
264 if args.window < 1:
265 parser.error("--window must be positive")
266 if args.compare and len(args.compare) < 2:
267 parser.error("--compare needs at least two records")
269 try:
270 reference = Reference(args.fasta)
271 except ReferenceError as exc:
272 print(f"error: {exc}", file=sys.stderr)
273 return 2
275 records: list[tuple[str, int, str, str]] = []
276 if args.compare:
277 try:
278 records = [parse_spec(spec) for spec in args.compare]
279 except VariantError as exc:
280 parser.error(str(exc))
281 elif args.input:
282 records = read_records(args.input, args.split)
283 elif args.contig and args.pos is not None and args.ref and args.alt:
284 alts = args.alt.split(",") if args.split else [args.alt]
285 records = [(args.contig, args.pos, args.ref, a) for a in alts]
286 else:
287 parser.error("give contig pos ref alt, or --input, or --compare")
289 rows = run(records, reference, args.window)
290 emit(rows, COLUMNS, args.format, args.output)
291 sys.stdout.flush()
293 if args.compare:
294 if any(row["ref_check"] != "ok" for row in rows):
295 print("\nverdict: cannot compare -- at least one record failed", file=sys.stderr)
296 return 1
297 keys = {row["normalized"] for row in rows}
298 if len(keys) == 1:
299 print(
300 f"\nverdict: identical -- all {len(rows)} records normalise to "
301 f"{keys.pop()}",
302 file=sys.stderr,
303 )
304 else:
305 print(
306 f"\nverdict: distinct -- {len(keys)} different variants after "
307 "normalisation",
308 file=sys.stderr,
309 )
311 return 1 if any(row["ref_check"] in {"MISMATCH", "error", "incomplete"} for row in rows) else 0
314if __name__ == "__main__":
315 raise SystemExit(main())