# Copyright 2026 Google LLC # # Licensed under the Apache License, Version 2.0 (the "License"); # you may not use this file except in compliance with the License. # You may obtain a copy of the License at # # http://www.apache.org/licenses/LICENSE-2.0 # # Unless required by applicable law or agreed to in writing, software # distributed under the License is distributed on an "AS IS" BASIS, # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. # See the License for the specific language governing permissions and # limitations under the License. """A command-line tool to query the Ensembl REST API. This script provides subcommands for gene/transcript/protein lookup, ID resolution and cross-referencing, sequence retrieval, and variant effect prediction (VEP). All rich data is saved to a temporary JSON file; a concise human-readable summary is printed to stdout. """ # /// script # requires-python = ">=3.10" # dependencies = [ # "scienceskillscommon", # ] # [tool.uv.sources] # scienceskillscommon = { path = "../../scienceskillscommon" } # /// import argparse import json import os import sys import tempfile from science_skills.skills.scienceskillscommon import http_client BASE_URL = "https://rest.ensembl.org" GRCH37_URL = "https://grch37.rest.ensembl.org" VEP_PLUGINS = ( "?AlphaMissense=1&Conservation=1&DosageSensitivity=1&IntAct=1" "&MaveDB=1&OpenTargets=1&LoF=loftee&NMD=1&UTRAnnotator=1" "&mutfunc=1&LOEUF=1" ) _CLIENT_REGULAR = http_client.HttpClient(BASE_URL, qps=15) _CLIENT_GRCH37 = http_client.HttpClient(GRCH37_URL, qps=15) def _get_client(assembly=None): """Return the correct client for the requested assembly.""" if assembly and assembly.upper() == "GRCH37": print("[*] Using GRCh37 assembly.") return _CLIENT_GRCH37 return _CLIENT_REGULAR def _get_species(args): """Return the species, defaulting to 'human' if not specified.""" species = args.species if not species: species = "human" print("[*] No species specified. Defaulting to 'human'.") return species def _save_json(data, prefix, output_path=None): """Write *data* as pretty-printed JSON to a temp file and return the path.""" if output_path: path = output_path with open(path, "w") as fh: json.dump(data, fh, indent=2) else: fd, path = tempfile.mkstemp( prefix=f"ensembl_{prefix}_", suffix=".json", text=True ) with os.fdopen(fd, "w") as fh: json.dump(data, fh, indent=2) return path def _try_fallback(url: str, query: str, client=None) -> list[dict[str, str]]: """Tries to resolve gene symbol via cross-references if initial lookup fails. Args: url: The URL for the cross-reference lookup. query: The original gene symbol query. client: The HttpClient to use. Defaults to _CLIENT_REGULAR. Returns: A list of dictionaries, where each dictionary contains gene information resolved from the cross-reference. """ if client is None: client = _CLIENT_REGULAR fallback_data = client.fetch_json(url) if not fallback_data: return [] matches = [] for item in fallback_data: if item.get("type") != "gene": continue ens_id = item.get("id") url_lookup = f"/lookup/id/{ens_id}" try: gene_data = client.fetch_json(url_lookup) matches.append(gene_data) except http_client.HttpError: print(f"[!] Failed to fetch details for resolved ID {ens_id}") matches.append({ "id": ens_id, "biotype": "N/A", "seq_region_name": "?", "start": "?", "end": "?", "strand": "?", }) if matches: print(f"[*] Resolved via synonym '{query}' to {len(matches)} gene(s):") return matches def cmd_resolve_gene(args): """Resolve a symbol / alias / RefSeq ID to one or more ENSG IDs.""" query = args.query species = _get_species(args) client = _get_client(args.assembly) url = f"/lookup/symbol/{species}/{query}?expand=0" try: data = client.fetch_json(url) matches = data if isinstance(data, list) else [data] print(f"[*] Symbol '{query}' resolved ({len(matches)} match(es)):") except http_client.HttpError as e: if e.status_code == 404 or ( e.status_code == 400 and e.body and b"No valid lookup found for symbol" in e.body ): print(f"[*] Symbol '{query}' not found. Trying synonym resolution...") url_fallback = f"/xrefs/symbol/{species}/{query}" try: data = _try_fallback(url_fallback, query, client=client) except http_client.HttpError as exc: raise e from exc if not data: raise e matches = data else: raise e for m in matches: eid = m.get("id", "N/A") biotype = m.get("biotype", "N/A") chrom = m.get("seq_region_name", "?") start = m.get("start", "?") end = m.get("end", "?") strand = m.get("strand", "?") print( f" ENSG: {eid} | biotype: {biotype} " f"| location: chr{chrom}:{start}-{end} (strand {strand})" ) path = _save_json(data, f"resolve_{query}", output_path=args.output) print(f"[*] Full JSON saved to {path}") def cmd_map_id(args): """Cross-reference an Ensembl ID to an external database.""" eid = args.id ext_db = args.external_db client = _get_client(args.assembly) params = f"?external_db={ext_db}" if ext_db else "" url = f"/xrefs/id/{eid}{params}" data = client.fetch_json(url) if not data: print(f"[*] No cross-references found for {eid}.") else: db_label = f" in {ext_db}" if ext_db else "" print(f"[*] {len(data)} cross-reference(s) for {eid}{db_label}:") for entry in data[:10]: primary = entry.get("primary_id", "N/A") display = entry.get("display_id", "N/A") dbname = entry.get("db_display_name", entry.get("dbname", "N/A")) print(f" {dbname}: {primary} ({display})") if len(data) > 10: print(f" … and {len(data) - 10} more (see JSON).") path = _save_json(data, f"mapid_{eid}", output_path=args.output) print(f"[*] Full JSON saved to {path}") def cmd_get_sequence(args): """Fetch raw genomic DNA for a coordinate window.""" coords = args.coords species = _get_species(args) assembly = args.assembly # Normalise the region string: accept chr17:100-200, 17:100-200, # 17:100..200 region = coords.replace(",", "").lower().removeprefix("chr") region = region.replace("-", "..") client = _get_client(assembly) url = f"/sequence/region/{species}/{region}?" if assembly: url += f"coord_system_version={assembly}&" headers = {"Accept": "text/plain"} seq = client.fetch_text(url, headers=headers) if args.output: path = args.output with open(path, "w") as fh: fh.write(str(seq)) else: fd, path = tempfile.mkstemp(prefix="ensembl_seq_", suffix=".txt", text=True) with os.fdopen(fd, "w") as fh: fh.write(str(seq)) length = len(str(seq)) print(f"[*] Fetched genomic sequence for {coords} (length: {length} bp).") print(f"[*] Sequence saved to {path}") def cmd_gene_summary(args): """Return high-level metadata for a gene by its ENSG ID.""" ensg = args.ensg_id client = _get_client(args.assembly) url = f"/lookup/id/{ensg}" data = client.fetch_json(url) symbol = data.get("display_name", "N/A") biotype = data.get("biotype", "N/A") desc = data.get("description", "N/A") chrom = data.get("seq_region_name", "?") start = data.get("start", "?") end = data.get("end", "?") strand = "+" if data.get("strand", 1) == 1 else "-" assembly = data.get("assembly_name", "N/A") print(f"[*] Gene summary for {ensg}:") print(f" Symbol: {symbol}") print(f" Biotype: {biotype}") print(f" Description: {desc}") print(f" Location: chr{chrom}:{start}-{end} ({strand})") print(f" Assembly: {assembly}") path = _save_json(data, f"gene_{ensg}", output_path=args.output) print(f"[*] Full JSON saved to {path}") def cmd_transcripts(args): """List transcripts for a gene, with optional MANE / canonical filtering.""" ensg = args.ensg_id client = _get_client(args.assembly) url = f"/lookup/id/{ensg}?expand=1;mane=1" data = client.fetch_json(url) transcripts = data.get("Transcript", []) if not transcripts: print(f"[*] No transcripts found for {ensg}.") path = _save_json(data, f"transcripts_{ensg}", output_path=args.output) print(f"[*] Full JSON saved to {path}") return # Apply filters filtered = transcripts if args.only_mane: filtered = [t for t in transcripts if t.get("MANE")] elif args.only_canonical: filtered = [t for t in transcripts if t.get("is_canonical") == 1] if not filtered and (args.only_mane or args.only_canonical): label = "MANE Select" if args.only_mane else "Canonical" print( f"[*] No {label} transcript found for {ensg}. " f"Total transcripts: {len(transcripts)}." ) return filter_label = "" if args.only_mane: filter_label = " (MANE Select only)" elif args.only_canonical: filter_label = " (Canonical only)" print( f"[*] {len(filtered)} transcript(s) for {ensg}{filter_label} " f"(from {len(transcripts)} total):" ) print() print("| Transcript ID | Biotype | TSL | Length (bp) | Flags |") print("| --- | --- | --- | --- | --- |") for t in filtered: tid = t.get("id", "N/A") biotype = t.get("biotype", "N/A") tsl = t.get("TSL", {}) tsl_val = tsl.get("value") if isinstance(tsl, dict) else tsl if tsl_val is None: tsl_val = "N/A" length = t.get("length", "N/A") flags = [] if t.get("is_canonical") == 1: flags.append("Canonical") mane = t.get("MANE") if mane: for m in mane: flags.append(m.get("type", "MANE")) flag_str = ", ".join(flags) if flags else "-" print(f"| {tid} | {biotype} | {tsl_val} | {length} | {flag_str} |") path = _save_json(data, f"transcripts_{ensg}", output_path=args.output) print(f"\n[*] Full JSON saved to {path}") def cmd_canonical_tss(args): """Retrieve the TSS for the canonical transcript of a gene.""" query = args.gene species = _get_species(args) client = _get_client(args.assembly) # 1. Resolve to ID if symbol ensg = query if not query.lower().startswith("ens"): url = f"/lookup/symbol/{species}/{query}?expand=0" data = client.fetch_json(url) if isinstance(data, list) and len(data) > 1: print( f"[!] Warning: '{query}' resolved to {len(data)} genes. " f"Using first match: {data[0].get('id')} " f"({data[0].get('display_name', 'N/A')}). " f"Other matches: {', '.join(d.get('id', '?') for d in data[1:])}" ) ensg = data.get("id") if not isinstance(data, list) else data[0].get("id") if not ensg: print( f"[!] Could not resolve symbol {query} to Ensembl ID.", file=sys.stderr, ) sys.exit(1) # 2. Fetch transcripts url = f"/lookup/id/{ensg}?expand=1;mane=1" data = client.fetch_json(url) transcripts = data.get("Transcript", []) if not transcripts: print(f"[!] No transcripts found for {ensg}.", file=sys.stderr) sys.exit(1) canonical = None for t in transcripts: if t.get("is_canonical") == 1: canonical = t break if not canonical: print(f"[!] No Canonical transcript found for {ensg}.", file=sys.stderr) sys.exit(1) strand = canonical.get("strand", 1) start = canonical.get("start") end = canonical.get("end") tid = canonical.get("id") if strand == 1: tss = start else: tss = end chrom = data.get("seq_region_name", "?") print(f"[*] Gene {ensg} (chr{chrom})") print(f"[*] Canonical Transcript: {tid}") print(f"[*] Strand: {'+' if strand == 1 else '-'}") print(f"[*] TSS Coordinate: {tss}") path = _save_json(canonical, f"canonical_tss_{ensg}", output_path=args.output) print(f"[*] Full canonical transcript JSON saved to {path}") def cmd_transcript_structure(args): """Return exon, CDS, and UTR layout for a transcript.""" enst = args.transcript_id client = _get_client(args.assembly) url = f"/lookup/id/{enst}?expand=1;mane=1" data = client.fetch_json(url) exons = data.get("Exon", []) strand = data.get("strand", 1) trans_start = data.get("start") trans_end = data.get("end") print(f"[*] Transcript structure for {enst}:") print(f" Biotype: {data.get('biotype', 'N/A')}") print( f" Genomic span: chr{data.get('seq_region_name', '?')}:" f"{trans_start}-{trans_end} (strand {'+' if strand == 1 else '-'})" ) print(f" Exons: {len(exons)}") translation = data.get("Translation") utr5 = None utr3 = None if translation: cds_start = translation.get("start") cds_end = translation.get("end") cds_length_aa = translation.get("length", "N/A") ensp = translation.get("id", "N/A") print(f" CDS: {cds_start}-{cds_end} ({cds_length_aa} aa, {ensp})") # Compute UTRs utr5 = ( (cds_end + 1, trans_end) if cds_end and trans_end and cds_end < trans_end else None ) utr3 = ( (trans_start, cds_start - 1) if cds_start and trans_start and cds_start > trans_start else None ) if strand == 1: utr3, utr5 = utr5, utr3 if utr5: print(f" 5' UTR: {utr5[0]}-{utr5[1]}") if utr3: print(f" 3' UTR: {utr3[0]}-{utr3[1]}") else: print(" (non-coding – no CDS/UTR)") if exons: print() print("| Exon # | ID | Start | End | Length (bp) |") print("| --- | --- | --- | --- | --- |") sorted_exons = sorted(exons, key=lambda e: e.get("start", 0)) for i, ex in enumerate(sorted_exons, 1): eid = ex.get("id", "N/A") estart = ex.get("start", "?") eend = ex.get("end", "?") elen = ( eend - estart + 1 if isinstance(estart, int) and isinstance(eend, int) else "?" ) print(f"| {i} | {eid} | {estart} | {eend} | {elen} |") # Enrich the saved data with computed UTR info if translation: data["_computed_utrs"] = {} if utr5: data["_computed_utrs"]["5_prime"] = {"start": utr5[0], "end": utr5[1]} if utr3: data["_computed_utrs"]["3_prime"] = {"start": utr3[0], "end": utr3[1]} path = _save_json(data, f"structure_{enst}", output_path=args.output) print(f"\n[*] Full JSON saved to {path}") def cmd_protein_info(args): """Fetch ENSP ID and sequence length for a transcript.""" enst = args.transcript_id client = _get_client(args.assembly) url = f"/lookup/id/{enst}?expand=1" data = client.fetch_json(url) translation = data.get("Translation") if not translation: print(f"[*] {enst} has no translation (likely non-coding).") path = _save_json(data, f"protein_{enst}", output_path=args.output) print(f"[*] Full JSON saved to {path}") return ensp = translation.get("id", "N/A") length = translation.get("length", "N/A") print(f"[*] Protein for {enst}:") print(f" ENSP: {ensp}") print(f" Length: {length} aa") path = _save_json(data, f"protein_{enst}", output_path=args.output) print(f"[*] Full JSON saved to {path}") def cmd_protein_sequence(args): """Fetch the amino acid sequence (FASTA) for a transcript or protein ID.""" target = args.id client = _get_client(args.assembly) url = f"/sequence/id/{target}?type=protein" headers = {"Accept": "text/x-fasta"} fasta = client.fetch_text(url, headers=headers) if args.output: path = args.output with open(path, "w") as fh: fh.write(str(fasta)) else: fd, path = tempfile.mkstemp( prefix=f"ensembl_protseq_{target}_", suffix=".fasta", text=True ) with os.fdopen(fd, "w") as fh: fh.write(str(fasta)) # Count sequence length (exclude header lines) lines = str(fasta).strip().split("\n") seq = "".join(l for l in lines if not l.startswith(">")) print(f"[*] Protein sequence for {target}: {len(seq)} aa") print(f"[*] FASTA saved to {path}") def cmd_vep(args): """Predict variant consequences using the Ensembl VEP.""" variant = args.variant_str species = _get_species(args) client = _get_client(args.assembly) if variant.startswith("rs"): url = f"/vep/{species}/id/{variant}{VEP_PLUGINS}" else: parts = variant.split(":") if len(parts) == 4: chrom, pos, ref, alt = parts end = int(pos) + len(ref) - 1 region = f"{chrom}:{pos}-{end}:1" url = f"/vep/{species}/region/{region}/{alt}{VEP_PLUGINS}" else: # Fallback: treat as HGVS url = f"/vep/{species}/hgvs/{variant}{VEP_PLUGINS}" data = client.fetch_json(url) if not isinstance(data, list) or not data: print("[!] No VEP results returned.") path = _save_json( data if data else {}, f"vep_{variant.replace(':', '_')}", output_path=args.output, ) print(f"[*] JSON saved to {path}") return top = data[0] mcv = top.get("most_severe_consequence", "Unknown") input_var = top.get("input", variant) t_conseq = top.get("transcript_consequences", []) print(f"[*] Variant: {input_var}") print(f"[*] Most severe consequence: {mcv}") print(f"[*] Found {len(t_conseq)} transcript consequences.") # Build the predictions table open_keys = { "am_class": "AlphaMissense Class", "am_pathogenicity": "AlphaMissense Pathogenicity", "conservation": "Conservation", "phaplo": "Dosage Sensitivity (Haplo)", "ptriplo": "Dosage Sensitivity (Triplo)", "lof": "Loss of Function (Loftee)", "nmd": "Nonsense-mediated Decay", "utr_annotator": "UTR Annotator", "mutfunc": "Mutfunc", "loeuf": "Loss of Function (LOEUF)", } rows = [] for tc in t_conseq: tid = tc.get("transcript_id", "Unknown") gene = tc.get("gene_symbol", "Unknown") terms = tc.get("consequence_terms", []) if terms: rows.append((tid, gene, "Consequence", ", ".join(terms))) aa = tc.get("amino_acids") if aa: rows.append((tid, gene, "Amino Acids", aa)) sift = tc.get("sift_prediction") if sift: score = tc.get("sift_score") val = f"{sift} ({score})" if score is not None else sift rows.append((tid, gene, "SIFT", val)) poly = tc.get("polyphen_prediction") if poly: score = tc.get("polyphen_score") val = f"{poly} ({score})" if score is not None else poly rows.append((tid, gene, "PolyPhen", val)) for key, label in open_keys.items(): val = tc.get(key) if val is None and "alphamissense" in tc: val = tc["alphamissense"].get(key) if val is not None: if isinstance(val, (list, dict)): val = json.dumps(val) rows.append((tid, gene, label, str(val))) if rows: print("\n[*] VEP Predictions Table:") print("| Transcript ID | Gene | Method/Metric | Value |") print("| --- | --- | --- | --- |") for r in rows: print(f"| {r[0]} | {r[1]} | {r[2]} | {r[3]} |") else: print("\n[*] No detailed predictions in transcript consequences.") safe = variant.replace(":", "_").replace(">", "_") path = _save_json(data, f"vep_{safe}", output_path=args.output) print(f"\n[*] Full JSON saved to {path}") def main(): """Parse CLI arguments and dispatch to the appropriate subcommand.""" parent_parser = argparse.ArgumentParser(add_help=False) parent_parser.add_argument( "--output", help="Output file path (optional)", ) parent_parser.add_argument( "--assembly", default=None, help="Assembly (e.g. GRCh38, GRCh37). Default: GRCh38.", ) parser = argparse.ArgumentParser( description="Query the Ensembl REST API.", formatter_class=argparse.RawDescriptionHelpFormatter, ) sub = parser.add_subparsers(dest="command", required=True) # ---- resolve-gene ---- p = sub.add_parser( "resolve-gene", parents=[parent_parser], help="Resolve a gene symbol / alias / RefSeq ID to ENSG ID(s).", ) p.add_argument("query", help="Gene symbol, alias, or RefSeq ID") p.add_argument( "--species", default=None, help="Species (defaults to 'human' if not specified)", ) p.set_defaults(func=cmd_resolve_gene) # ---- map-id ---- p = sub.add_parser( "map-id", parents=[parent_parser], help="Cross-reference an Ensembl ID to external databases.", ) p.add_argument("id", help="Ensembl ID (ENSG, ENST, ENSP)") p.add_argument( "--external-db", dest="external_db", default=None, help="Filter by external DB (e.g., UniProt, HGNC, RefSeq_mRNA)", ) p.set_defaults(func=cmd_map_id) # ---- get-sequence ---- p = sub.add_parser( "get-sequence", parents=[parent_parser], help="Fetch raw genomic DNA for a coordinate window.", ) p.add_argument( "coords", help="Genomic region, e.g. 17:7661779-7687550 or chr17:7661779-7687550", ) p.add_argument( "--species", default=None, help="Species (defaults to 'human' if not specified)", ) p.set_defaults(func=cmd_get_sequence) # ---- gene-summary ---- p = sub.add_parser( "gene-summary", parents=[parent_parser], help="Get high-level metadata for a gene (symbol, biotype, location).", ) p.add_argument("ensg_id", help="Ensembl gene ID (e.g. ENSG00000141510)") p.set_defaults(func=cmd_gene_summary) # ---- transcripts ---- p = sub.add_parser( "transcripts", parents=[parent_parser], help="List transcripts for a gene, with optional MANE/canonical filter.", ) p.add_argument("ensg_id", help="Ensembl gene ID") grp = p.add_mutually_exclusive_group() grp.add_argument( "--only-mane", action="store_true", help="Return only the MANE Select transcript (human only).", ) grp.add_argument( "--only-canonical", action="store_true", help="Return only the Ensembl Canonical transcript.", ) p.set_defaults(func=cmd_transcripts) # ---- canonical-tss ---- p = sub.add_parser( "canonical-tss", parents=[parent_parser], help="Get TSS coordinate for the canonical transcript of a gene.", ) p.add_argument("gene", help="Gene symbol or Ensembl ID") p.add_argument( "--species", default=None, help="Species (defaults to 'human' if not specified)", ) p.set_defaults(func=cmd_canonical_tss) # ---- transcript-structure ---- p = sub.add_parser( "transcript-structure", parents=[parent_parser], help="Get exon, CDS, and UTR layout for a transcript.", ) p.add_argument( "transcript_id", help="Ensembl transcript ID (e.g. ENST00000269305)" ) p.set_defaults(func=cmd_transcript_structure) # ---- protein-info ---- p = sub.add_parser( "protein-info", parents=[parent_parser], help="Get ENSP ID and sequence length for a transcript.", ) p.add_argument( "transcript_id", help="Ensembl transcript ID (e.g. ENST00000269305)" ) p.set_defaults(func=cmd_protein_info) # ---- protein-sequence ---- p = sub.add_parser( "protein-sequence", parents=[parent_parser], help="Fetch the amino acid sequence (FASTA) for an ENST or ENSP.", ) p.add_argument("id", help="Ensembl transcript or protein ID") p.set_defaults(func=cmd_protein_sequence) # ---- vep ---- p = sub.add_parser( "vep", parents=[parent_parser], help="Predict variant consequences (VEP) with open-license plugins.", ) p.add_argument( "variant_str", help="Variant as chr:pos:ref:alt (e.g. 9:21971147:T:C) or rsID", ) p.add_argument( "--species", default=None, help="Species (defaults to 'human' if not specified)", ) p.set_defaults(func=cmd_vep) args = parser.parse_args() args.func(args) if __name__ == "__main__": main()