Skills › Scientific Agent Skills › onekgpd › onekgpd_api.py
onekgpd_api.py
844 lines · skills/onekgpd/scripts/onekgpd_api.py at 5f2bc03
1# /// script
2# requires-python = ">=3.11"
3# dependencies = ["dnaerys>=0.2.1,<0.3.0"]
4# ///
5"""OneKGPd — individual-level queries over the 1000 Genomes Project.
7A single command-line wrapper exposing ten subcommands over the 1000 Genomes
8Project cohort (3,202 whole-genome-sequenced individuals, GRCh38): selecting and
9counting variants in a region (cohort-wide or within named individuals),
10selecting and counting the individuals who carry matching variants, homozygous-
11reference queries at a single position, pairwise relatedness, and dataset totals.
13Every command writes full JSON to a file (``--output``, or a temp file by
14default) and prints a concise human-readable summary to stdout. Coordinates are
15GRCh38, 1-based inclusive; resolve a gene/feature to coordinates against an
16authoritative source BEFORE querying.
18Examples
19--------
20 uv run scripts/onekgpd_api.py dataset-info
21 uv run scripts/onekgpd_api.py count-samples --chrom chr17 --start 43044295 --end 43170327 \
22 --consequence MISSENSE_VARIANT --alpha-missense-class AM_LIKELY_PATHOGENIC
23 uv run scripts/onekgpd_api.py kinship --sample1 NA19238 --sample2 NA19240
25MIT License. Author: Dnaerys.
26"""
28from __future__ import annotations
30import argparse
31import json
32import math
33import struct
34import os
35import sys
36import tempfile
37import time
38import warnings
39from typing import Any, NoReturn
41from dnaerys import (
42 AnnotationFilter,
43 DnaerysError,
44 DnaerysIncompleteResultWarning,
45 DnaerysClient,
46 Region,
47)
49# ---------------------------------------------------------------------------
50# Constants
51# ---------------------------------------------------------------------------
53DEFAULT_ENDPOINT = "db.dnaerys.org:443" # public 1000 Genomes instance (fixed)
54DEFAULT_TIMEOUT = 30.0 # seconds per RPC, including internal paging calls
55DEFAULT_VARIANT_LIMIT = 200 # hard cap when neither --limit nor --page-size given
56MAX_RETRIES = 3 # bounded retry attempts on retryable errors
57RETRY_BASE_DELAY = 1.0 # seconds; exponential backoff: 1, 2, ...
58PREVIEW_ROWS = 10 # rows shown in a stdout summary preview
60INCOMPLETE_NOTE = "[!] Result may be incomplete: some data was unreachable."
63# ---------------------------------------------------------------------------
64# Small generic helpers
65# ---------------------------------------------------------------------------
68def _fail(msg: str) -> NoReturn:
69 """Print a clean one-line message to stderr and exit non-zero."""
70 print(msg, file=sys.stderr)
71 sys.exit(1)
74def _split_csv(s: str) -> list[str]:
75 """Split a comma-separated value into a list of non-empty trimmed tokens."""
76 return [part.strip() for part in s.split(",") if part.strip()]
79def _positive_finite(value: str | float) -> float:
80 number = float(value)
81 if not math.isfinite(number) or number <= 0:
82 raise argparse.ArgumentTypeError("must be a positive finite number")
83 return number
86def _client(args) -> DnaerysClient:
87 """Keep the public TLS endpoint and GRCh38 explicit; bound each RPC."""
88 return DnaerysClient(
89 DEFAULT_ENDPOINT, tls=True, assembly="GRCh38",
90 default_timeout=getattr(args, "timeout", DEFAULT_TIMEOUT),
91 )
94def _sample_names(csv: str) -> list[str]:
95 names = _split_csv(csv)
96 if not names:
97 raise ValueError("--samples must name at least one individual; an empty list would query the whole cohort")
98 return names
101def _save_json(data: Any, prefix: str, output_path: str | None = None) -> str:
102 """Write *data* as indented JSON to *output_path* or a temp file; return path."""
103 if output_path:
104 path = output_path
105 with open(path, "w") as fh:
106 json.dump(data, fh, indent=2)
107 else:
108 fd, path = tempfile.mkstemp(prefix=f"onekgpd_{prefix}_", suffix=".json", text=True)
109 with os.fdopen(fd, "w") as fh:
110 json.dump(data, fh, indent=2)
111 return path
114def _emit(data: Any, prefix: str, summary_lines: list[str], output: str | None) -> None:
115 """Save full JSON to a file and print the summary + saved-path line to stdout."""
116 path = _save_json(data, prefix, output)
117 for line in summary_lines:
118 print(line)
119 print(f"[*] Full JSON saved to {path}")
122def _call_with_retry(fn):
123 """Run *fn*; on a retryable ``DnaerysError`` retry with bounded backoff.
125 The whole fetch-and-materialize closure is retried (a fresh client per
126 attempt), so streaming errors that surface during iteration are covered.
127 """
128 for attempt in range(1, MAX_RETRIES + 1):
129 try:
130 return fn()
131 except DnaerysError as e:
132 if e.is_retryable and attempt < MAX_RETRIES:
133 time.sleep(RETRY_BASE_DELAY * (2 ** (attempt - 1)))
134 continue
135 raise
138# ---------------------------------------------------------------------------
139# Input builders
140# ---------------------------------------------------------------------------
143def _parse_region_str(s: str) -> Region:
144 """Parse ``CHR:START-END`` into a ``Region`` (coordinate validation deferred)."""
145 try:
146 chrom, span = s.split(":", 1)
147 start_s, end_s = span.split("-", 1)
148 start, end = int(start_s), int(end_s)
149 except ValueError:
150 raise ValueError(
151 f"invalid --region {s!r}; expected CHR:START-END, "
152 "e.g. chr17:43044295-43170327"
153 )
154 return Region(chrom, start, end)
157def _build_regions(args) -> tuple[Region | None, list[Region] | None]:
158 """Return exactly one of (single Region, None) or (None, list[Region])."""
159 has_single = args.chrom is not None
160 has_multi = bool(args.region)
161 if has_single and has_multi:
162 raise ValueError("use either --chrom/--start/--end or --region, not both")
163 if not has_single and not has_multi:
164 raise ValueError(
165 "a region is required: pass --chrom/--start/--end or "
166 "one or more --region CHR:START-END"
167 )
168 if has_multi:
169 if args.ref or args.alt:
170 raise ValueError(
171 "--ref/--alt apply only to a single --chrom/--start/--end region"
172 )
173 return None, [_parse_region_str(r) for r in args.region]
174 if args.start is None or args.end is None:
175 raise ValueError("--chrom requires --start and --end")
176 return Region(args.chrom, args.start, args.end, ref=args.ref, alt=args.alt), None
179def _zygosity(args) -> tuple[bool, bool]:
180 """Return (hom, het); default both True, narrowed by --het-only/--hom-only."""
181 if getattr(args, "het_only", False):
182 return (False, True)
183 if getattr(args, "hom_only", False):
184 return (True, False)
185 return (True, True)
188_CSV_FIELDS = [
189 ("clin_significance", "clin_significance"),
190 ("consequence", "consequence"),
191 ("impact", "impact"),
192 ("variant_type", "variant_type"),
193 ("feature_type", "feature_type"),
194 ("bio_type", "bio_type"),
195 ("alpha_missense_class", "am_class"),
196]
197_FLOAT_FIELDS = [
198 ("af_lt", "af_lt"),
199 ("af_gt", "af_gt"),
200 ("gnomad_exomes_af_lt", "gnomad_exomes_af_lt"),
201 ("gnomad_exomes_af_gt", "gnomad_exomes_af_gt"),
202 ("gnomad_genomes_af_lt", "gnomad_genomes_af_lt"),
203 ("gnomad_genomes_af_gt", "gnomad_genomes_af_gt"),
204 ("alpha_missense_score_lt", "am_score_lt"),
205 ("alpha_missense_score_gt", "am_score_gt"),
206]
207_BOOL_FIELDS = [
208 ("biallelic_only", "biallelic_only"),
209 ("multiallelic_only", "multiallelic_only"),
210 ("exclude_males", "exclude_males"),
211 ("exclude_females", "exclude_females"),
212]
215def _build_annotation_filter(args) -> AnnotationFilter | None:
216 """Map the shared annotation flags to an ``AnnotationFilter`` (or None).
218 Note: ``--min-len-bp``/``--max-len-bp`` are NOT filter fields; they are passed
219 as the client method's ``variant_min_length``/``variant_max_length`` kwargs.
220 """
221 kwargs: dict[str, Any] = {}
222 for arg_name, field in _CSV_FIELDS:
223 raw = getattr(args, arg_name, None)
224 if raw is not None:
225 values = _split_csv(raw)
226 if not values:
227 raise ValueError(f"--{arg_name.replace('_', '-')} must contain at least one term")
228 kwargs[field] = values
229 for arg_name, field in _FLOAT_FIELDS:
230 val = getattr(args, arg_name, None)
231 if val is not None:
232 if not math.isfinite(val) or not 0 < val <= 1:
233 raise ValueError(
234 f"--{arg_name.replace('_', '-')} must be finite and in (0, 1]; "
235 "zero means an unset server filter. Use an explicit positive bound "
236 "or post-filter complete variant results for exact > 0."
237 )
238 if struct.unpack("!f", struct.pack("!f", val))[0] == 0:
239 raise ValueError(f"--{arg_name.replace('_', '-')} underflows the API float32 bound to zero")
240 kwargs[field] = val
241 for arg_name, field in _BOOL_FIELDS:
242 if getattr(args, arg_name, False):
243 kwargs[field] = True
245 # am-class vs am-score cannot be one argparse group (am-score is a range pair).
246 if kwargs.get("am_class") and ("am_score_lt" in kwargs or "am_score_gt" in kwargs):
247 _fail(
248 "Error: --alpha-missense-class cannot be combined with "
249 "--alpha-missense-score-lt/--alpha-missense-score-gt"
250 )
252 if not kwargs:
253 return None
254 return AnnotationFilter(**kwargs)
257# ---------------------------------------------------------------------------
258# Display / serialization helpers
259# ---------------------------------------------------------------------------
262def _chr_to_str(chrom) -> str:
263 """Render a ``Chromosome`` enum as ``chr17`` / ``chrX`` / ``chrMT``."""
264 return "chr" + chrom.name[len("CHR"):]
267def _region_one_label(r: Region) -> str:
268 return f"{_chr_to_str(r.chr)}:{r.start}-{r.end}"
271def _region_label(region: Region | None, regions: list[Region] | None) -> str:
272 if region is not None:
273 return _region_one_label(region)
274 return ", ".join(_region_one_label(r) for r in regions or [])
277def _zyg_label(hom: bool, het: bool) -> str:
278 if hom and het:
279 return "hom+het"
280 if het:
281 return "het only"
282 return "hom only"
285def _variant_to_dict(v) -> dict:
286 """Serialize a ``Variant`` to its in-scope output keys (enum rendered as text)."""
287 return {
288 "chr": _chr_to_str(v.chr),
289 "start": v.start,
290 "end": v.end,
291 "ref": v.ref,
292 "alt": v.alt,
293 "af": v.af,
294 "ac": v.ac,
295 "an": v.an,
296 "hom_samples": v.hom_samples,
297 "het_samples": v.het_samples,
298 "mis_samples": v.mis_samples,
299 "hom_samples_fx": v.hom_samples_fx,
300 "het_samples_fx": v.het_samples_fx,
301 "mis_samples_fx": v.mis_samples_fx,
302 "hom_samples_mxy": v.hom_samples_mxy,
303 "het_samples_mxy": v.het_samples_mxy,
304 "mis_samples_mxy": v.mis_samples_mxy,
305 "gnomad_exomes_af": v.gnomad_exomes_af,
306 "gnomad_genomes_af": v.gnomad_genomes_af,
307 "am_score": v.am_score,
308 "amino_acids": v.amino_acids,
309 "biallelic": v.biallelic,
310 }
313def _filters_echo(args) -> dict:
314 """Echo the annotation/length flags that were actually set (for provenance)."""
315 out: dict[str, Any] = {}
316 for arg_name, _ in _FLOAT_FIELDS:
317 val = getattr(args, arg_name, None)
318 if val is not None:
319 out[arg_name] = val
320 for arg_name, _ in _CSV_FIELDS:
321 raw = getattr(args, arg_name, None)
322 if raw:
323 out[arg_name] = _split_csv(raw)
324 for arg_name, _ in _BOOL_FIELDS:
325 if getattr(args, arg_name, False):
326 out[arg_name] = True
327 if getattr(args, "min_len_bp", None) is not None:
328 out["variant_min_length"] = args.min_len_bp
329 if getattr(args, "max_len_bp", None) is not None:
330 out["variant_max_length"] = args.max_len_bp
331 return out
334def _request_echo(args, *, samples, region, regions, hom, het) -> dict:
335 """Build the request-provenance block included in every region-command JSON."""
336 req: dict[str, Any] = {}
337 if region is not None:
338 req["region"] = _region_one_label(region)
339 if region.ref:
340 req["ref"] = region.ref
341 if region.alt:
342 req["alt"] = region.alt
343 if regions is not None:
344 req["regions"] = [_region_one_label(r) for r in regions]
345 if samples is not None:
346 req["samples"] = samples
347 req["zygosity"] = _zyg_label(hom, het)
348 filt = _filters_echo(args)
349 if filt:
350 req["filters"] = filt
351 return req
354# ---------------------------------------------------------------------------
355# Command workers / handlers
356# ---------------------------------------------------------------------------
359def _run_count_variants(args, samples: list[str] | None) -> None:
360 region, regions = _build_regions(args)
361 hom, het = _zygosity(args)
362 ann = _build_annotation_filter(args)
364 def fetch():
365 with _client(args) as client:
366 return client.count_variants(
367 region=region,
368 regions=regions,
369 samples=samples,
370 hom=hom,
371 het=het,
372 annotations=ann,
373 variant_min_length=args.min_len_bp,
374 variant_max_length=args.max_len_bp,
375 )
377 result = _call_with_retry(fetch)
378 incomplete = result.metadata.affected
379 data = {
380 "command": args.command,
381 "count": result.count,
382 "request": _request_echo(
383 args, samples=samples, region=region, regions=regions, hom=hom, het=het
384 ),
385 "result_incomplete": incomplete,
386 }
387 summary = [
388 f"{result.count:,} variants match in {_region_label(region, regions)} "
389 f"({_zyg_label(hom, het)})"
390 ]
391 if incomplete:
392 summary.append(INCOMPLETE_NOTE)
393 _emit(data, args.command.replace("-", "_"), summary, args.output)
396def _run_select_variants(args, samples: list[str] | None) -> None:
397 region, regions = _build_regions(args)
398 hom, het = _zygosity(args)
399 ann = _build_annotation_filter(args)
400 page_size = args.page_size
401 limit = args.limit
403 def fetch():
404 with _client(args) as client:
405 if page_size is not None:
406 pq = client.paginate_variants(
407 page_size=page_size,
408 region=region,
409 regions=regions,
410 samples=samples,
411 hom=hom,
412 het=het,
413 annotations=ann,
414 variant_min_length=args.min_len_bp,
415 variant_max_length=args.max_len_bp,
416 )
417 collected = []
418 for page in pq:
419 collected.extend(page.variants)
420 return collected, pq.metadata, False
421 stream = client.select_variants(
422 region=region,
423 regions=regions,
424 samples=samples,
425 hom=hom,
426 het=het,
427 annotations=ann,
428 variant_min_length=args.min_len_bp,
429 variant_max_length=args.max_len_bp,
430 limit=limit,
431 )
432 collected = stream.to_list()
433 truncated = limit is not None and len(collected) >= limit
434 return collected, stream.metadata, truncated
436 variants, meta, truncated = _call_with_retry(fetch)
437 incomplete = meta.affected
438 request = _request_echo(
439 args, samples=samples, region=region, regions=regions, hom=hom, het=het
440 )
441 if page_size is not None:
442 request["page_size"] = page_size
443 else:
444 request["limit"] = limit
445 data = {
446 "command": args.command,
447 "count_returned": len(variants),
448 "truncated": truncated,
449 "request": request,
450 "result_incomplete": incomplete,
451 "variants": [_variant_to_dict(v) for v in variants],
452 }
453 summary = [
454 f"Returned {len(variants)} variants in {_region_label(region, regions)} "
455 f"({_zyg_label(hom, het)})"
456 ]
457 for v in variants[:PREVIEW_ROWS]:
458 summary.append(
459 f" {_chr_to_str(v.chr)}:{v.start} {v.ref}>{v.alt} "
460 f"af={v.af:.6g} am_score={v.am_score:.6g} aa={v.amino_acids or '-'}"
461 )
462 if not variants:
463 summary.append(" (no variants matched)")
464 if truncated:
465 summary.append(
466 f"[!] Truncated at --limit {limit}; raise --limit or use "
467 "--page-size to retrieve the full set."
468 )
469 if incomplete:
470 summary.append(INCOMPLETE_NOTE)
471 _emit(data, args.command.replace("-", "_"), summary, args.output)
474def cmd_count_variants(args) -> None:
475 _run_count_variants(args, samples=None)
478def cmd_count_variants_in_samples(args) -> None:
479 _run_count_variants(args, samples=_sample_names(args.samples))
482def cmd_select_variants(args) -> None:
483 _run_select_variants(args, samples=None)
486def cmd_select_variants_in_samples(args) -> None:
487 _run_select_variants(args, samples=_sample_names(args.samples))
490def cmd_count_samples(args) -> None:
491 region, regions = _build_regions(args)
492 hom, het = _zygosity(args)
493 ann = _build_annotation_filter(args)
495 def fetch():
496 with _client(args) as client:
497 return client.count_samples(
498 region=region,
499 regions=regions,
500 hom=hom,
501 het=het,
502 annotations=ann,
503 variant_min_length=args.min_len_bp,
504 variant_max_length=args.max_len_bp,
505 )
507 result = _call_with_retry(fetch)
508 incomplete = result.metadata.affected
509 data = {
510 "command": args.command,
511 "count": result.count,
512 "request": _request_echo(
513 args, samples=None, region=region, regions=regions, hom=hom, het=het
514 ),
515 "result_incomplete": incomplete,
516 }
517 summary = [
518 f"{result.count:,} individuals carry a matching variant in "
519 f"{_region_label(region, regions)} ({_zyg_label(hom, het)})"
520 ]
521 if incomplete:
522 summary.append(INCOMPLETE_NOTE)
523 _emit(data, "count_samples", summary, args.output)
526def cmd_select_samples(args) -> None:
527 region, regions = _build_regions(args)
528 hom, het = _zygosity(args)
529 ann = _build_annotation_filter(args)
531 def fetch():
532 with _client(args) as client:
533 return client.select_samples(
534 region=region,
535 regions=regions,
536 hom=hom,
537 het=het,
538 annotations=ann,
539 variant_min_length=args.min_len_bp,
540 variant_max_length=args.max_len_bp,
541 skip=args.skip,
542 limit=args.limit,
543 )
545 result = _call_with_retry(fetch)
546 names = list(result.samples)
547 incomplete = result.metadata.affected
548 request = _request_echo(
549 args, samples=None, region=region, regions=regions, hom=hom, het=het
550 )
551 if args.skip is not None:
552 request["skip"] = args.skip
553 if args.limit is not None:
554 request["limit"] = args.limit
555 data = {
556 "command": args.command,
557 "count": len(names),
558 "samples": names,
559 "request": request,
560 "result_incomplete": incomplete,
561 }
562 summary = [
563 f"{len(names)} individuals carry a matching variant in "
564 f"{_region_label(region, regions)} ({_zyg_label(hom, het)})"
565 ]
566 for name in names[:PREVIEW_ROWS]:
567 summary.append(f" {name}")
568 if not names:
569 summary.append(" (none)")
570 if incomplete:
571 summary.append(INCOMPLETE_NOTE)
572 _emit(data, "select_samples", summary, args.output)
575def cmd_count_samples_hom_ref(args) -> None:
576 Region(args.chrom, args.position, args.position) # validate 1-based position
577 def fetch():
578 with _client(args) as client:
579 return client.count_samples_hom_ref(chr=args.chrom, position=args.position)
581 result = _call_with_retry(fetch)
582 count = result.count
583 incomplete = result.metadata.affected
584 present = None if incomplete else count != -1
585 pos = f"{args.chrom}:{args.position}"
586 data = {
587 "command": args.command,
588 "count": count,
589 "variant_present": present,
590 "result_incomplete": incomplete,
591 "request": {"chrom": args.chrom, "position": args.position},
592 }
593 if incomplete:
594 summary = [f"Homozygous-reference response at {pos}: {count}; presence and count are not definitive.", INCOMPLETE_NOTE]
595 elif count == -1:
596 summary = [
597 f"No variant exists at {pos} in the dataset; "
598 "homozygous-reference count is undefined here."
599 ]
600 elif count == 0:
601 summary = [
602 f"A variant exists at {pos}, but no individual is homozygous reference."
603 ]
604 else:
605 summary = [f"{count:,} individuals are homozygous reference at {pos}."]
606 _emit(data, "count_samples_hom_ref", summary, args.output)
609def cmd_select_samples_hom_ref(args) -> None:
610 Region(args.chrom, args.position, args.position) # validate 1-based position
611 def fetch():
612 with _client(args) as client:
613 return client.select_samples_hom_ref(chr=args.chrom, position=args.position)
615 result = _call_with_retry(fetch)
616 names = list(result.samples)
617 incomplete = result.metadata.affected
618 pos = f"{args.chrom}:{args.position}"
619 data = {
620 "command": args.command,
621 "result_incomplete": incomplete,
622 "count": len(names),
623 "samples": names,
624 "request": {"chrom": args.chrom, "position": args.position},
625 }
626 summary = [f"{len(names)} individuals are homozygous reference at {pos}"]
627 for name in names[:PREVIEW_ROWS]:
628 summary.append(f" {name}")
629 if not names:
630 summary.append(" (none)")
631 if incomplete:
632 summary.append(INCOMPLETE_NOTE)
633 _emit(data, "select_samples_hom_ref", summary, args.output)
636def cmd_kinship(args) -> None:
637 def fetch():
638 with _client(args) as client:
639 return client.kinship_duo(sample1=args.sample1, sample2=args.sample2)
641 result = _call_with_retry(fetch)
642 if not result.pairs:
643 _fail(f"Error: no relatedness result returned for {args.sample1}, {args.sample2}")
644 pair = result.pairs[0]
645 degree = pair.degree.name
646 phi = pair.phi_bwf
647 incomplete = result.metadata.affected
648 data = {
649 "command": "kinship",
650 "sample1": pair.sample1,
651 "sample2": pair.sample2,
652 "degree": degree,
653 "phi_bwf": phi,
654 "result_incomplete": incomplete,
655 }
656 summary = [
657 f"{pair.sample1} <-> {pair.sample2}: {degree} "
658 f"(KING kinship coefficient phi = {phi:.4f})"
659 ]
660 if incomplete:
661 summary.append(INCOMPLETE_NOTE)
662 _emit(data, "kinship", summary, args.output)
665def cmd_dataset_info(args) -> None:
666 def fetch():
667 with _client(args) as client:
668 return client.dataset_info()
670 info = _call_with_retry(fetch)
671 data = {
672 "command": "dataset-info",
673 "samples_total": info.samples_total,
674 "females_total": info.females_total,
675 "males_total": info.males_total,
676 "variants_total": info.variants_total,
677 "assembly": info.assembly.name,
678 "cohorts": [
679 {
680 "cohort_name": c.cohort_name,
681 "samples_count": c.samples_count,
682 "female_count": c.female_count,
683 "male_count": c.male_count,
684 "synthetic": c.synthetic,
685 }
686 for c in info.cohorts
687 ],
688 }
689 summary = [
690 f"1000 Genomes Project - {info.samples_total:,} individuals "
691 f"(F: {info.females_total:,}, M: {info.males_total:,}); "
692 f"{info.variants_total:,} variants; {info.assembly.name}"
693 ]
694 for c in info.cohorts:
695 summary.append(
696 f" {c.cohort_name}: {c.samples_count:,} "
697 f"(F {c.female_count:,} / M {c.male_count:,})"
698 )
699 _emit(data, "dataset_info", summary, args.output)
702# ---------------------------------------------------------------------------
703# Argument parsing
704# ---------------------------------------------------------------------------
707def _add_annotation_flags(p: argparse.ArgumentParser) -> None:
708 g = p.add_argument_group("annotation filters (different fields AND; CSV within a field OR)")
709 g.add_argument("--af-lt", type=float, help="Keep variants with 1000 Genomes AF < this value.")
710 g.add_argument("--af-gt", type=float, help="Keep variants with 1000 Genomes AF > this value.")
711 g.add_argument("--gnomad-exomes-af-lt", type=float, help="gnomAD v4.1 exomes AF < this value.")
712 g.add_argument("--gnomad-exomes-af-gt", type=float, help="gnomAD v4.1 exomes AF > this value.")
713 g.add_argument("--gnomad-genomes-af-lt", type=float, help="gnomAD v4.1 genomes AF < this value.")
714 g.add_argument("--gnomad-genomes-af-gt", type=float, help="gnomAD v4.1 genomes AF > this value.")
715 g.add_argument("--clin-significance", metavar="CSV", help="ClinVar significance terms (CSV). See references/annotation_vocabularies.md.")
716 g.add_argument("--consequence", metavar="CSV", help="Sequence Ontology consequence terms (CSV).")
717 g.add_argument("--impact", metavar="CSV", help="VEP impact (CSV: HIGH,MODERATE,LOW,MODIFIER).")
718 g.add_argument("--variant-type", metavar="CSV", help="SO variant class terms (CSV).")
719 g.add_argument("--feature-type", metavar="CSV", help="VEP feature types (CSV).")
720 g.add_argument("--bio-type", metavar="CSV", help="VEP biotypes (CSV).")
721 g.add_argument("--alpha-missense-class", metavar="CSV", help="AM_LIKELY_BENIGN,AM_LIKELY_PATHOGENIC,AM_AMBIGUOUS (CSV).")
722 g.add_argument("--alpha-missense-score-lt", type=float, help="AlphaMissense score < this value.")
723 g.add_argument("--alpha-missense-score-gt", type=float, help="AlphaMissense score > this value.")
725 bm = p.add_mutually_exclusive_group()
726 bm.add_argument("--biallelic-only", action="store_true", help="Keep only biallelic sites.")
727 bm.add_argument("--multiallelic-only", action="store_true", help="Keep only multiallelic sites.")
729 ef = p.add_mutually_exclusive_group()
730 ef.add_argument("--exclude-males", action="store_true", help="Exclude male samples.")
731 ef.add_argument("--exclude-females", action="store_true", help="Exclude female samples.")
733 z = p.add_mutually_exclusive_group()
734 z.add_argument("--het-only", action="store_true", help="Heterozygous carriage only (default: both).")
735 z.add_argument("--hom-only", action="store_true", help="Homozygous carriage only (default: both).")
738def build_parser() -> argparse.ArgumentParser:
739 conn_parser = argparse.ArgumentParser(add_help=False)
740 conn_parser.add_argument(
741 "--timeout", type=_positive_finite, default=DEFAULT_TIMEOUT,
742 help="Positive finite seconds per RPC (default 30); paging and retries may take longer overall.",
743 )
744 conn_parser.add_argument(
745 "--output",
746 help="Write full JSON to this path (default: a temp file in the system temp dir).",
747 )
749 region_parser = argparse.ArgumentParser(add_help=False)
750 region_parser.add_argument("--chrom", help="Chromosome, e.g. chr17, 17, X, MT (single-region mode).")
751 region_parser.add_argument("--start", type=int, help="1-based inclusive start (with --chrom).")
752 region_parser.add_argument("--end", type=int, help="1-based inclusive end (with --chrom).")
753 region_parser.add_argument("--ref", help="Reference allele to narrow to one allele (single-region only).")
754 region_parser.add_argument("--alt", help="Alternate allele to narrow to one allele (single-region only).")
755 region_parser.add_argument(
756 "--region", action="append", metavar="CHR:START-END",
757 help="A region as CHR:START-END; repeat for multiple regions (multi-region mode).",
758 )
759 region_parser.add_argument("--min-len-bp", type=int, help="Minimum alternate-allele length (bp).")
760 region_parser.add_argument("--max-len-bp", type=int, help="Maximum alternate-allele length (bp).")
762 annot_parser = argparse.ArgumentParser(add_help=False)
763 _add_annotation_flags(annot_parser)
765 parser = argparse.ArgumentParser(
766 prog="onekgpd_api.py",
767 description="Individual-level queries over the 1000 Genomes Project "
768 "(3,202 WGS individuals, GRCh38).",
769 )
770 sub = parser.add_subparsers(dest="command", required=True)
772 region_parents = [conn_parser, region_parser, annot_parser]
774 # dataset-info
775 p = sub.add_parser("dataset-info", parents=[conn_parser], help="Dataset totals (samples, sex split, variants, assembly).")
776 p.set_defaults(func=cmd_dataset_info)
778 # count-variants
779 p = sub.add_parser("count-variants", parents=region_parents, help="Count variants in a region, cohort-wide.")
780 p.set_defaults(func=cmd_count_variants)
782 # select-variants
783 p = sub.add_parser("select-variants", parents=region_parents, help="Select variants in a region, cohort-wide.")
784 lp = p.add_mutually_exclusive_group()
785 lp.add_argument("--limit", type=int, default=DEFAULT_VARIANT_LIMIT, help=f"Hard cap on returned variants (default {DEFAULT_VARIANT_LIMIT}).")
786 lp.add_argument("--page-size", type=int, help="Retrieve ALL matching variants in pages of this size (full walk).")
787 p.set_defaults(func=cmd_select_variants)
789 # count-variants-in-samples
790 p = sub.add_parser("count-variants-in-samples", parents=region_parents, help="Count variants in a region within named individuals.")
791 p.add_argument("--samples", required=True, metavar="CSV", help="Comma-separated sample names.")
792 p.set_defaults(func=cmd_count_variants_in_samples)
794 # select-variants-in-samples
795 p = sub.add_parser("select-variants-in-samples", parents=region_parents, help="Select variants in a region within named individuals.")
796 p.add_argument("--samples", required=True, metavar="CSV", help="Comma-separated sample names.")
797 lp = p.add_mutually_exclusive_group()
798 lp.add_argument("--limit", type=int, default=DEFAULT_VARIANT_LIMIT, help=f"Hard cap on returned variants (default {DEFAULT_VARIANT_LIMIT}).")
799 lp.add_argument("--page-size", type=int, help="Retrieve ALL matching variants in pages of this size (full walk).")
800 p.set_defaults(func=cmd_select_variants_in_samples)
802 # count-samples
803 p = sub.add_parser("count-samples", parents=region_parents, help="Count individuals carrying a matching variant in a region.")
804 p.set_defaults(func=cmd_count_samples)
806 # select-samples
807 p = sub.add_parser("select-samples", parents=region_parents, help="List individuals carrying a matching variant in a region.")
808 p.add_argument("--skip", type=int, help="Skip the first N individuals.")
809 p.add_argument("--limit", type=int, help="Return at most N individuals.")
810 p.set_defaults(func=cmd_select_samples)
812 # count-samples-hom-ref
813 p = sub.add_parser("count-samples-hom-ref", parents=[conn_parser], help="Count individuals homozygous reference at a single position.")
814 p.add_argument("--chrom", required=True, help="Chromosome, e.g. chr17, 17, X, MT.")
815 p.add_argument("--position", type=int, required=True, help="1-based position.")
816 p.set_defaults(func=cmd_count_samples_hom_ref)
818 # select-samples-hom-ref
819 p = sub.add_parser("select-samples-hom-ref", parents=[conn_parser], help="List individuals homozygous reference at a single position.")
820 p.add_argument("--chrom", required=True, help="Chromosome, e.g. chr17, 17, X, MT.")
821 p.add_argument("--position", type=int, required=True, help="1-based position.")
822 p.set_defaults(func=cmd_select_samples_hom_ref)
824 # kinship
825 p = sub.add_parser("kinship", parents=[conn_parser], help="Relatedness (degree + KING coefficient) between two individuals.")
826 p.add_argument("--sample1", required=True, help="First sample name.")
827 p.add_argument("--sample2", required=True, help="Second sample name.")
828 p.set_defaults(func=cmd_kinship)
830 return parser
833def main(argv: list[str] | None = None) -> None:
834 # Incompleteness is surfaced via result metadata and the summary, not warnings.
835 warnings.simplefilter("ignore", DnaerysIncompleteResultWarning)
836 args = build_parser().parse_args(argv)
837 try:
838 args.func(args)
839 except (DnaerysError, ValueError, OSError) as e:
840 _fail(f"Error: {e}")
843if __name__ == "__main__":
844 main()