Skills › Scientific Agent Skills › pathogen-variant-surveillance › mutation_profile.py
mutation_profile.py
233 lines · skills/pathogen-variant-surveillance/scripts/mutation_profile.py at 5f2bc03
1#!/usr/bin/env python3
2# /// script
3# requires-python = ">=3.11"
4# dependencies = []
5# ///
6"""Mutations carried by a lineage, or the difference between two lineages.
8Use this to answer "what distinguishes this lineage" and "does my assay target
9still match" from current sequences, instead of from a lineage's founding
10description -- which drifts as sublineages accumulate substitutions.
12 python3 mutation_profile.py XFG.1.1 --gene S --min-proportion 0.9
13 python3 mutation_profile.py XFG.23.1.3 --versus XFG.1.1 --gene S
14 python3 mutation_profile.py 2.3.4.4b --instance h5n1 --gene HA --min-proportion 0.8
15 python3 mutation_profile.py XFG --sublineages --nucleotide --since 2026-01-01
17``proportion`` is over ``coverage`` -- the sequences that actually resolved that
18site -- not over every matching sequence. A site with poor coverage can show a
19high proportion on very few sequences, so read the coverage column before quoting a
20proportion.
21"""
22from __future__ import annotations
24import argparse
25import sys
26from datetime import date
28from lapis_client import (
29 LapisError,
30 count,
31 data_version,
32 describe_instance,
33 emit,
34 lineage_filter,
35 mutations,
36 pick_date_field,
37 pick_lineage_field,
38 range_keys,
39 request,
40 resolve_base_url,
41 surveillance_filters,
42)
44PROFILE_COLUMNS = ("mutation", "gene", "position", "from", "to", "proportion", "count", "coverage")
45DIFF_COLUMNS = ("mutation", "gene", "position", "verdict", "prop_a", "prop_b", "n_a", "n_b", "coverage_a", "coverage_b")
48def parse_where(pairs: list[str], schema: dict) -> dict[str, str]:
49 """Turn ``KEY=VALUE`` arguments into filters, checking keys against the schema."""
50 filters: dict[str, str] = {}
51 types = schema.get("types", {})
52 for pair in pairs:
53 key, sep, value = pair.partition("=")
54 if not sep:
55 raise LapisError(f"--where expects KEY=VALUE, got {pair!r}")
56 if key.split(".")[0] not in types:
57 raise LapisError(f"{key!r} is not a field on this instance")
58 filters[key] = value
59 return filters
62def profile(base_url: str, filters: dict, *, amino_acid: bool, min_proportion: float,
63 gene: str | None) -> dict[str, dict]:
64 """Mutation rows keyed by mutation string, optionally restricted to one gene."""
65 rows = mutations(base_url, filters, amino_acid=amino_acid, min_proportion=min_proportion)
66 keyed: dict[str, dict] = {}
67 for row in rows:
68 if gene and str(row.get("sequenceName") or "").upper() != gene.upper():
69 continue
70 keyed[row["mutation"]] = row
71 return keyed
74def build_parser() -> argparse.ArgumentParser:
75 parser = argparse.ArgumentParser(description="Mutation profile of a lineage.")
76 parser.add_argument("lineage", help="lineage or clade name; 'NAME*' includes descendants")
77 parser.add_argument("--versus", help="second lineage to diff against")
78 parser.add_argument("--instance", default="sars-cov-2", help="registry name (default: sars-cov-2)")
79 parser.add_argument("--base-url", help="any other LAPIS deployment")
80 parser.add_argument("--lineage-field", help="override the auto-detected lineage column")
81 parser.add_argument("--date-field", help="override the auto-detected collection-date column")
82 parser.add_argument("--gene", help="restrict to one gene or segment (e.g. S, HA, seg4)")
83 parser.add_argument("--nucleotide", action="store_true",
84 help="nucleotide instead of amino-acid mutations")
85 parser.add_argument("--min-proportion", type=float, default=0.8,
86 help="report mutations at or above this proportion (default: 0.8)")
87 parser.add_argument("--sublineages", action="store_true", help="include descendant lineages")
88 parser.add_argument("--since", help="restrict to sequences collected on or after YYYY-MM-DD")
89 parser.add_argument("--where", action="append", default=[], metavar="KEY=VALUE",
90 help="extra filter, repeatable")
91 parser.add_argument("--format", choices=("table", "tsv", "json"), default="table")
92 return parser
95def main(argv: list[str] | None = None) -> int:
96 args = build_parser().parse_args(argv)
98 try:
99 if not 0 < args.min_proportion <= 1:
100 raise LapisError("--min-proportion must be in (0, 1]")
101 base_url = resolve_base_url(args.instance, args.base_url)
102 schema = describe_instance(base_url)
103 lineage_field, has_index = pick_lineage_field(schema, args.lineage_field)
104 where = surveillance_filters(schema, parse_where(args.where, schema))
105 if lineage_field in where:
106 raise LapisError("set the lineage with the positional argument, not --where")
107 if args.since:
108 date.fromisoformat(args.since)
109 date_field = pick_date_field(schema, "collection", args.date_field)
110 where[range_keys(date_field)[0]] = args.since
111 if args.gene:
112 reference = request(base_url, "sample/referenceGenome")
113 names = [r["name"] for r in reference["nucleotideSequences" if args.nucleotide else "genes"]]
114 matches = [n for n in names if n.upper() == args.gene.upper()]
115 if len(matches) != 1:
116 raise LapisError(f"unknown {'segment' if args.nucleotide else 'gene'} {args.gene!r}; choose from {names}")
117 args.gene = matches[0]
119 # Fetch unpruned mutation rows for comparisons. A missing row still
120 # has unknown coverage and cannot establish a zero proportion.
121 diff_floor = 0.0
123 primary = lineage_filter(args.lineage, has_index, args.sublineages)
124 base_filters = {**where, lineage_field: primary}
125 n_primary = count(base_url, base_filters)
126 if n_primary == 0:
127 print(
128 f"error: no sequences match {lineage_field}={primary}"
129 + (f" with {where}" if where else "")
130 + ". Check the name with resolve_lineage.py.",
131 file=sys.stderr,
132 )
133 return 1
135 first = profile(
136 base_url,
137 base_filters,
138 amino_acid=not args.nucleotide,
139 min_proportion=args.min_proportion if not args.versus else diff_floor,
140 gene=args.gene,
141 )
143 if not args.versus:
144 rows = [
145 {
146 "mutation": m["mutation"],
147 "gene": m.get("sequenceName", ""),
148 "position": m.get("position", ""),
149 "from": m.get("mutationFrom", ""),
150 "to": m.get("mutationTo", ""),
151 "proportion": f"{m['proportion']:.3f}",
152 "count": m.get("count", ""),
153 "coverage": m.get("coverage", ""),
154 }
155 for m in sorted(
156 first.values(),
157 key=lambda r: (str(r.get("sequenceName")), int(r.get("position") or 0)),
158 )
159 ]
160 summary = (
161 f"\n# {schema['name']} via {base_url} | data version {data_version(base_url)}"
162 f"\n# {lineage_field}={primary} | {n_primary} sequences"
163 f" | {'nucleotide' if args.nucleotide else 'amino acid'} mutations"
164 f" at proportion >= {args.min_proportion}"
165 f"\n# proportion is over per-site coverage, not over all {n_primary} sequences"
166 )
167 else:
168 secondary = lineage_filter(args.versus, has_index, args.sublineages)
169 n_secondary = count(base_url, {**where, lineage_field: secondary})
170 if n_secondary == 0:
171 print(f"error: no sequences match {lineage_field}={secondary}", file=sys.stderr)
172 return 1
173 second = profile(
174 base_url,
175 {**where, lineage_field: secondary},
176 amino_acid=not args.nucleotide,
177 min_proportion=diff_floor,
178 gene=args.gene,
179 )
181 rows = []
182 for mutation in sorted(set(first) | set(second)):
183 a = first.get(mutation, {})
184 b = second.get(mutation, {})
185 pa = float(a["proportion"]) if a and a.get("coverage", 0) > 0 else None
186 pb = float(b["proportion"]) if b and b.get("coverage", 0) > 0 else None
187 if max(pa or 0.0, pb or 0.0) < args.min_proportion:
188 continue
189 if pa is None or pb is None:
190 verdict = "not_comparable"
191 elif pa >= args.min_proportion and pb < args.min_proportion:
192 verdict = "above_a_only"
193 elif pb >= args.min_proportion and pa < args.min_proportion:
194 verdict = "above_b_only"
195 elif pa >= args.min_proportion and pb >= args.min_proportion:
196 verdict = "above_both"
197 else:
198 continue
199 source = a or b
200 rows.append(
201 {
202 "mutation": mutation,
203 "gene": source.get("sequenceName", ""),
204 "position": source.get("position", ""),
205 "verdict": verdict,
206 "prop_a": f"{pa:.3f}" if pa is not None else None,
207 "prop_b": f"{pb:.3f}" if pb is not None else None,
208 "n_a": a.get("count"),
209 "n_b": b.get("count"),
210 "coverage_a": a.get("coverage"),
211 "coverage_b": b.get("coverage"),
212 }
213 )
214 rows.sort(key=lambda r: (str(r["gene"]), int(r["position"] or 0)))
215 summary = (
216 f"\n# {schema['name']} via {base_url} | data version {data_version(base_url)}"
217 f"\n# A = {primary} ({n_primary} sequences), B = {secondary} ({n_secondary})"
218 f"\n# verdict is relative to a {args.min_proportion} proportion threshold; not evolutionary gain/loss"
219 f"\n# absent rows have unknown proportion/coverage, never an inferred zero"
220 )
221 summary += f"\n# filters {where} | gene/segment {args.gene or 'all'} | lineage field {lineage_field}"
222 except (LapisError, ValueError) as exc:
223 print(f"error: {exc}", file=sys.stderr)
224 return 2
226 print(emit(rows, DIFF_COLUMNS if args.versus else PROFILE_COLUMNS, args.format))
227 sys.stdout.flush()
228 print(summary, file=sys.stderr)
229 return 0
232if __name__ == "__main__":
233 sys.exit(main())