Skills › Scientific Agent Skills › pathogen-variant-surveillance › resolve_lineage.py
resolve_lineage.py
207 lines · skills/pathogen-variant-surveillance/scripts/resolve_lineage.py at 5f2bc03
1#!/usr/bin/env python3
2# /// script
3# requires-python = ">=3.11"
4# dependencies = []
5# ///
6"""Resolve a lineage name against the live nomenclature before trusting it.
8Pango names are not stable identifiers. They are minted continuously, aliased
9through a key that must be fetched to be read, and **withdrawn or redesignated
10after the fact** -- so a remembered lineage fact is not merely stale, it can be
11actively wrong. This script answers, for each name: does it still exist, what
12does it expand to, what is it descended from, and how many sequences carry it.
14 python3 resolve_lineage.py XFG.23.1.3 PQ.17 PC.2
15 python3 resolve_lineage.py XFG --descendants
16 python3 resolve_lineage.py 2.3.4.4b --instance h5n1
18Exit code is 1 when any name is withdrawn, unknown or unverified, so it gates
19a manuscript's lineage list.
20"""
21from __future__ import annotations
23import argparse
24import re
25import sys
27from lapis_client import (
28 LapisError,
29 SnapshotChanged,
30 children_map,
31 count,
32 data_version,
33 descendants,
34 describe_instance,
35 emit,
36 fetch_lineage_notes,
37 fetch_pango_aliases,
38 lineage_definition,
39 lineage_field_candidates,
40 parent_chain,
41 pango_provenance,
42 pick_lineage_field,
43 recombinant_parents,
44 resolve_base_url,
45 surveillance_filters,
46 unalias_full,
47)
49COLUMNS = (
50 "query",
51 "status",
52 "unaliased",
53 "parent",
54 "recombinant_of",
55 "descendants",
56 "sequences",
57 "detail",
58)
59PANGO_FIELDS = {"pangoLineage", "nextcladePangoLineage"}
60REDESIGNATED = re.compile(r"[Rr]edesignated as ([A-Za-z][A-Za-z0-9.]*)")
63def normalise(name: str) -> str:
64 """Uppercase the alphabetic head so ``xfg.1`` matches ``XFG.1``."""
65 head, sep, tail = name.strip().partition(".")
66 return f"{head.upper()}{sep}{tail}"
69def build_parser() -> argparse.ArgumentParser:
70 parser = argparse.ArgumentParser(
71 description="Check lineage names against the live nomenclature.",
72 )
73 parser.add_argument("names", nargs="+", help="lineage or clade names to resolve")
74 parser.add_argument("--instance", default="sars-cov-2", help="registry name (default: sars-cov-2)")
75 parser.add_argument("--base-url", help="any other LAPIS deployment")
76 parser.add_argument("--lineage-field", help="override the auto-detected lineage column")
77 parser.add_argument("--descendants", action="store_true",
78 help="list every descendant lineage instead of counting them")
79 parser.add_argument("--no-counts", action="store_true",
80 help="skip the sequence counts (one fewer request per name)")
81 parser.add_argument("--format", choices=("table", "tsv", "json"), default="table")
82 return parser
85def main(argv: list[str] | None = None) -> int:
86 args = build_parser().parse_args(argv)
88 try:
89 base_url = resolve_base_url(args.instance, args.base_url)
90 schema = describe_instance(base_url)
91 lineage_field, has_index = pick_lineage_field(schema, args.lineage_field)
92 # Fetched here, not inside the provenance f-string: a network failure at
93 # print time would raise after the table had already been written.
94 version = data_version(base_url)
95 except LapisError as exc:
96 print(f"error: {exc}", file=sys.stderr)
97 return 2
99 definition: dict = {}
100 if has_index:
101 try:
102 definition = lineage_definition(base_url, lineage_field)
103 except LapisError as exc:
104 print(f"error: no lineage definition for {lineage_field}: {exc}", file=sys.stderr)
105 return 2
107 aliases: dict = {}
108 notes: dict = {}
109 if lineage_field in PANGO_FIELDS:
110 try:
111 aliases = fetch_pango_aliases()
112 notes = fetch_lineage_notes()
113 except LapisError as exc:
114 print(f"error: pango-designation unreachable: {exc}", file=sys.stderr)
115 return 2
117 # Inverted once: the SARS-CoV-2 definition holds thousands of entries, so
118 # rebuilding it per name makes a list of lineages quadratic.
119 children = children_map(definition)
121 rows: list[dict] = []
122 failures = 0
123 for raw in args.names:
124 name = normalise(raw) if lineage_field in PANGO_FIELDS else raw.strip()
125 note = notes.get(name, {})
126 in_definition = name in definition
128 if note.get("status") == "withdrawn":
129 status = "withdrawn"
130 elif note.get("status") == "designated" or (in_definition and lineage_field not in PANGO_FIELDS):
131 status = "current"
132 elif notes or definition:
133 status = "unknown"
134 else:
135 status = "unverified"
136 if status in {"withdrawn", "unknown", "unverified"}:
137 failures += 1
139 detail = note.get("note", "")
140 successor = REDESIGNATED.search(detail)
141 if successor:
142 detail = f"now {successor.group(1)}; {detail}"
143 elif status == "unknown":
144 detail = "no such name in the live nomenclature for this instance"
146 chain = parent_chain(name, definition) if definition else []
147 kids = descendants(name, definition, children) if definition else []
149 sequences: object = ""
150 if not args.no_counts:
151 try:
152 query = f"{name}*" if has_index else name
153 sequences = count(base_url, surveillance_filters(schema, {lineage_field: query}))
154 except SnapshotChanged as exc:
155 print(f"error: {exc}", file=sys.stderr)
156 return 2
157 except LapisError as exc:
158 # An indexed column rejects an unknown lineage outright rather
159 # than answering 0 -- which is the one place a typo is caught
160 # for you. An unindexed column would have returned 0 instead.
161 sequences = "n/a" if "not a valid lineage" in str(exc) else "error"
162 if sequences == "error":
163 failures += 1
165 rows.append(
166 {
167 "query": raw,
168 "status": status,
169 "unaliased": unalias_full(name, aliases) if aliases else "",
170 "parent": chain[0] if chain else "",
171 "recombinant_of": "+".join(recombinant_parents(unalias_full(name, aliases), aliases)) if aliases else "",
172 "descendants": ", ".join(kids) if args.descendants else len(kids),
173 "sequences": sequences,
174 "detail": detail,
175 }
176 )
178 print(emit(rows, COLUMNS, args.format))
179 sys.stdout.flush()
181 print(
182 f"\n# {schema['name']} via {base_url} | data version {version}"
183 f"\n# lineage column {lineage_field}"
184 f"{' with a lineage index' if has_index else ' with no lineage index'}"
185 + (
186 f"\n# nomenclature from pango-designation ({len(notes)} names, "
187 f"{sum(1 for v in notes.values() if v['status'] == 'withdrawn')} withdrawn)"
188 f"\n# source content {pango_provenance()}"
189 if notes
190 else ""
191 )
192 + ("\n# 'sequences' counts the lineage and its descendants" if has_index else "")
193 + f"\n# record filters {surveillance_filters(schema, {})}",
194 file=sys.stderr,
195 )
196 others = [n for n, _ in lineage_field_candidates(schema) if n != lineage_field]
197 if others:
198 print(
199 f"# other lineage-like columns here: {', '.join(others)} "
200 f"(select one with --lineage-field)",
201 file=sys.stderr,
202 )
203 return 1 if failures else 0
206if __name__ == "__main__":
207 sys.exit(main())