commit 95742963fcf22e6176b42085013b3913afb5bb0d Author: wehub-skill-sync Date: Mon Jul 13 21:35:31 2026 +0800 chore: import zh skill ensembl-database diff --git a/README.wehub.md b/README.wehub.md new file mode 100644 index 0000000..f2dd5f6 --- /dev/null +++ b/README.wehub.md @@ -0,0 +1,9 @@ +# WeHub 来源说明 + +- Skill 名称:`ensembl-database` +- 中文类目:序列拉取与跨库 ID 解析 +- 上游仓库:`google-deepmind__science-skills` +- 上游路径:`skills/ensembl_database/SKILL.md` +- 上游链接:https://github.com/google-deepmind/science-skills/blob/HEAD/skills/ensembl_database/SKILL.md +- 本仓库为 WeHub 中文 Skill 汉化包,基于 skill 市场筛选 Top200 清单整理 +- 原作者、版权和许可证信息以上游仓库为准 diff --git a/SKILL.md b/SKILL.md new file mode 100644 index 0000000..cea28dc --- /dev/null +++ b/SKILL.md @@ -0,0 +1,142 @@ +--- +name: ensembl-database +description: > + 查询 Ensembl 数据库以解析基因、转录本和蛋白质 ID,获取基因组或蛋白质序列, + 检索基因结构(外显子),以及获取变异后果与效应预测(VEP)。使用本技能作为 + 主要的 ID 翻译器、基因组序列数据库和变异效应预测工具。 +--- + +# Ensembl 数据库:ID 映射与基因组特征 + +## 前置条件 + +1. **`uv`**:阅读 `uv` 技能并遵循其设置说明,以确保 `uv` 已安装并位于 PATH 中。 +2. **用户通知**:如果本技能目录中尚未存在 `LICENSE_NOTIFICATION.txt`,则 (1) 显著提示用户查看 https://useast.ensembl.org/index.html 和 https://github.com/Ensembl/ensembl-rest/wiki 上的条款,然后 (2) 创建记录通知文本和时间戳的文件。 + +## 概述 + +Ensembl 数据库是基因组注释的资源。本技能允许你与 Ensembl REST API 交互,以解析模糊的符号、交叉引用 ID(RefSeq、HGNC、UniProt、ENSG)、获取原始序列,以及检索详细的转录本结构。 + +**关键概念:** + +- **ENSG(基因):** 人类基因的稳定标识符。其他物种将具有不同的三位字母物种代码。 +- **ENST(转录本):** 转录本(剪接亚型)的稳定标识符。 +- **ENSP(蛋白质):** 翻译后蛋白质的稳定标识符。 +- **MANE Select:** Ensembl 和 NCBI 共同认定的共识主要转录本。 +- **Canonical:** Ensembl 的代表性转录本(当 MANE 不可用或非人类物种时使用)。 + +## 核心规则 + +- **使用封装脚本**:始终执行提供的辅助脚本来查询数据库,而非直接访问数据库。这些脚本会自动优雅地强制执行所需的速率限制。 +- **默认物种:** 如果提示中物种缺失或模糊,默认使用 `"human"`。你必须明确向用户标记此默认行为,以确保其知晓。 +- **主要转录本:** 列出基因的转录本时,仅返回 MANE Select 转录本(人类)或 Canonical 转录本(其他物种),除非用户明确要求所有替代亚型。当存在多个转录本且你默认使用主要转录本时,必须向用户标记说明。 +- **基因组组装处理:** 默认组装为 GRCh38。对于 GRCh37 请求,必须使用 `--assembly GRCh37` 标志。当使用非默认组装时,必须明确向用户标记说明。 +- **输出位置:** 脚本默认将完整的 JSON/FASTA 输出写入 `/tmp` 下的临时文件,或使用 `--output` 标志写入用户指定的文件。同时会在标准输出打印简明摘要。 +- **通知**:如果使用了本技能,请确保在输出中提及。 + +### 可用命令 + +**1. 解析基因 ID** — 将符号、别名或 RefSeq ID 解析为一个或多个 ENSG ID。如果未找到主要符号,则自动回退到解析同义词。 + +```bash +uv run scripts/ensembl_api.py resolve-gene TP53 --species human --output tp53.json +uv run scripts/ensembl_api.py resolve-gene PCL2 --output pcl2.json # 回退到同义词解析 +``` + +**2. 映射 ID 到外部数据库** — 将 Ensembl ID 交叉引用到 UniProt、HGNC、RefSeq 等。 + +```bash +uv run scripts/ensembl_api.py map-id ENSG00000141510 --external-db UniProt --output uniprot_map.json +uv run scripts/ensembl_api.py map-id ENST00000269305 --external-db RefSeq_mRNA --output refseq_map.json +``` + +**3. 获取基因组序列** — 获取坐标窗口的原始 DNA。支持通过 `--assembly GRCh37` 使用 GRCh37。 + +```bash +uv run scripts/ensembl_api.py get-sequence 17:7661779-7687550 --species human --output seq.txt +uv run scripts/ensembl_api.py get-sequence chr9:21971100-21971200 --assembly GRCh37 --output seq_grch37.txt +``` + +**4. 基因摘要** — 高层级元数据:符号、生物类型、描述、染色体位置。 + +```bash +uv run scripts/ensembl_api.py gene-summary ENSG00000141510 --output gene_summary.json +``` + +**5. 列出转录本** — 基因的所有转录本,可选 `--only-mane` 或 `--only-canonical` 过滤。输出包含转录本支持等级(TSL)。 + +```bash +uv run scripts/ensembl_api.py transcripts ENSG00000141510 --only-mane --output transcripts_mane.json +uv run scripts/ensembl_api.py transcripts ENSG00000141510 --only-canonical --output transcripts_canonical.json +uv run scripts/ensembl_api.py transcripts ENSG00000141510 --output transcripts_all.json +``` + +**5b. 规范转录起始位点(Canonical TSS)** — 获取基因规范转录本的转录起始位点(TSS)的单一坐标。 + +> [!NOTE] 与标准的 `transcripts` 命令不同,`canonical-tss` 接受符号(例如 `TP53`)和 Ensembl ID,并自动解析它们。它还会计算链方向(TSS 对于 `+` 链为 `Start`,对于 `-` 链为 `End`),直接输出单一的整数坐标。 + +```bash +uv run scripts/ensembl_api.py canonical-tss TP53 --output tp53_tss.json +uv run scripts/ensembl_api.py canonical-tss ENSG00000141510 --output tss.json +``` + +**6. 转录本结构** — 转录本的外显子坐标、CDS 边界以及计算出的 5'/3' UTR 区域。 + +```bash +uv run scripts/ensembl_api.py transcript-structure ENST00000269305 --output structure.json +``` + +**7. 蛋白质信息** — 转录本的 ENSP ID 和序列长度。 + +```bash +uv run scripts/ensembl_api.py protein-info ENST00000269305 --output protein_info.json +``` + +**8. 蛋白质序列** — 转录本(ENST)或蛋白质(ENSP)ID 的氨基酸 FASTA 序列。 + +```bash +uv run scripts/ensembl_api.py protein-sequence ENST00000269305 --output protein.fasta +uv run scripts/ensembl_api.py protein-sequence ENSP00000269305 --output protein_ensp.fasta +``` + +**9. 变异后果(VEP)** — 预测基因组变异的分子后果。包含开放许可插件:AlphaMissense、Conservation、DosageSensitivity、IntAct、MaveDB、OpenTargets、LoF(Loftee)、NMD、UTRAnnotator、mutfunc、LOEUF。 + +```bash +uv run scripts/ensembl_api.py vep 9:21971147:T:C --species human --output vep.json +uv run scripts/ensembl_api.py vep rs699 --species human --output vep_rs699.json +``` + +VEP 标准输出示例: + +``` +[*] 变异:9:21971147:T>C +[*] 最严重的后果:missense_variant +[*] 发现 15 个转录本后果。 + +[*] VEP 预测: + + - ENST00000304494(CDKN2A):后果 = missense_variant + - ENST00000304494(CDKN2A):氨基酸 = N/S + - ENST00000304494(CDKN2A):SIFT = 有害(0.01) + - ENST00000304494(CDKN2A):AlphaMissense 分类 = likely_benign + - ENST00000304494(CDKN2A):AlphaMissense 致病性 = 0.2129 + - ENST00000304494(CDKN2A):保守性 = 2.05 + - ENST00000304494(CDKN2A):剂量敏感性(单倍体不足) = 0.889228328567991 + - ENST00000304494(CDKN2A):剂量敏感性(三倍体) = 0.135514349094646 + - ENST00000304494(CDKN2A):功能缺失(LOEUF) = 0.791 +``` + +**展示 VEP 结果:** 运行 VEP 命令后,你必须将标准输出中的完整 VEP 预测列表呈现给用户。该列表包含标准的 VEP 预测(后果、氨基酸、SIFT、PolyPhen)和开放许可插件结果(AlphaMissense、保守性、剂量敏感性、LOEUF、Loftee LoF、NMD、UTRAnnotator、Mutfunc)。不要只做总结——展示完整列表,以便用户看到所有预测。如果列表很长(许多转录本),则完整显示 MANE Select / canonical 转录本行,并注明完整数据在 JSON 输出中。 + +## 解析输出 + +如果用户需要未在标准输出中总结的详细嵌套结构数据(例如转录本外显子 2 的精确整数坐标): + +1. 找到 JSON 文件(通过 `--output` 指定的文件或脚本打印的临时文件路径)。 +2. 使用 `jq` 等终端工具或编写一个简短的、一次性 Python 脚本来提取所请求的特定数据点。如果 JSON 文件非常大,**不要**尝试将整个文件读入你的上下文。 + +## 自定义查询 + +如果你需要进行脚本不支持 API 调用(例如获取蛋白质结构域注释、组装之间的坐标映射、同源搜索、连锁不平衡或表型查询),请阅读 `references/ensembl_rest_api_reference.md` 以获取可用端点、参数和响应字段的完整参考。 + +**关键:** 编写自定义脚本或使用提供脚本的替代方案时,你**必须**遵守 Ensembl REST API 速率限制(最大每秒 15 个请求),并优雅地处理 `429 Too Many Requests` 错误(例如使用指数退避)。 diff --git a/references/ensembl_rest_api_reference.md b/references/ensembl_rest_api_reference.md new file mode 100644 index 0000000..4b4834e --- /dev/null +++ b/references/ensembl_rest_api_reference.md @@ -0,0 +1,277 @@ +# Ensembl REST API 参考文档 + +本文档提供 Ensembl REST API(`https://rest.ensembl.org`)的简明参考。当 `ensembl_api.py` 脚本无法覆盖特定用例时,可用此文档构建自定义查询。 + +## 通用约定 + +- **基础 URL:** `https://rest.ensembl.org`(GRCh38)。GRCh37 请使用:`https://grch37.rest.ensembl.org` +- **内容协商:** 设置 `Content-Type` 请求头来控制响应格式: + - `application/json` — 结构化 JSON(大部分端点的默认格式) + - `text/plain` — 原始序列字符串 + - `text/x-fasta` — FASTA 格式的序列 + - `text/x-gff3` — GFF3 注释输出 +- **速率限制:** 最大 15 次请求/秒。收到 HTTP 429 时,请遵守 `Retry-After` 响应头。 +- **区域格式:** `CHR:START..END:STRAND`,其中 STRAND 为 `1`(正向)或 `-1`(反向)。大多数端点也支持连字符格式(`START-END`)。 +- **查询参数分隔符:** Ensembl 使用 `;`(分号)分隔查询参数,例如 `?expand=1;mane=1`。标准的 `&` 同样适用。 + +--- + +## 查询端点 + +### `GET /lookup/id/{id}` +查询任意 Ensembl 稳定 ID(基因、转录本、蛋白质)并获取元数据。 + +- **`expand`**(0/1):包含子对象(基因的 Transcript 数组,转录本的 Exon 数组,编码转录本的 Translation) +- **`mane`**(0/1):在转录本上包含 MANE Select/Plus Clinical 注释 +- **`db_type`**(字符串):数据库(默认:`core`)。可选值:`core`、`otherfeatures` +- **`format`**(字符串):`full`(默认)或 `condensed` +- **`species`**(字符串):当 ID 不明确时覆盖物种信息 + +**关键响应字段(基因):** +`id`、`display_name`(符号)、`biotype`、`description`、`seq_region_name`(染色体)、`start`、`end`、`strand`、`assembly_name`、`Transcript[]`(展开后)。 + +**关键响应字段(转录本,展开后):** +`id`、`biotype`、`display_name`、`is_canonical`(0 或 1)、`length`、`MANE[]`(包含 `type` 的数组,值为 `MANE_Select` 或 `MANE_Plus_Clinical`)、`TSL`(Transcript Support Level 对象,含 `value`)、`Exon[]`、`Translation`(含 `id`、`start`、`end`、`length`)。 + +### `GET /lookup/symbol/{species}/{symbol}` +将基因符号解析为其 Ensembl 稳定 ID。 + +- **`expand`**(0/1):包含子对象 + +返回与 `/lookup/id/` 相同的结构。 + +### `POST /lookup/id` +批量查询:以 JSON 体形式发送 `{"ids": ["ENSG...", "ENST..."]}`。返回以 ID 为键的字典。 + +--- + +## 交叉引用端点 + +### `GET /xrefs/id/{id}` +获取 Ensembl ID 的外部数据库引用。 + +- **`external_db`**(字符串):按数据库名称过滤(例如 `UniProt`、`HGNC`、`RefSeq_mRNA`、`UCSC`、`EntrezGene`) +- **`all_levels`**(0/1):包含来自父/子特征的交叉引用 + +**响应:** 对象数组,包含 `primary_id`、`display_id`、`db_display_name`、`dbname`、`description`、`info_type`。 + +### `GET /xrefs/symbol/{species}/{symbol}` +查找与外部符号匹配的 Ensembl ID。 + +### `GET /xrefs/name/{species}/{name}` +更广泛的搜索——跨所有外部数据库查询任意名称。 + +--- + +## 序列端点 + +### `GET /sequence/id/{id}` +通过稳定 ID 获取 Ensembl 特征的序列。 + +- **`type`**(字符串):`genomic`(默认)、`cdna`、`cds`、`protein` +- **`expand_5prime`**(整数):向上游延伸 N 个碱基 +- **`expand_3prime`**(整数):向下游延伸 N 个碱基 +- **`mask`**(字符串):屏蔽方式:`hard` 或 `soft` + +设置 `Accept: text/x-fasta` 获取 FASTA 格式输出,`text/plain` 获取原始字符串。 + +### `GET /sequence/region/{species}/{region}` +获取坐标窗口内的基因组 DNA。 + +- **`coord_system_version`**(字符串):组装版本(例如 `GRCh38`) +- **`expand_5prime`**(整数):向上游延伸 N 个碱基 +- **`expand_3prime`**(整数):向下游延伸 N 个碱基 +- **`mask`**(字符串):`hard` 或 `soft` 重复序列屏蔽 +- **`mask_feature`**(0/1):应用特征级屏蔽 + +区域格式:`CHR:START..END:STRAND`(例如 `X:1000000..1000100:1`)。 + +### `POST /sequence/region/{species}` +批量:发送 `{"regions": ["X:1000..2000", "7:100..200"]}`。 + +--- + +## 重叠端点 + +### `GET /overlap/region/{species}/{region}` +查找与基因组区域重叠的特征。适用于查找某基因座的基因、窗口内的变异或调控特征。 + +- **`feature`**(字符串):要返回的特征类型。可重复使用以指定多个类型。可选值:`gene`、`transcript`、`cds`、`exon`、`repeat`、`simple`、`misc`、`variation`、`somatic_variation`、`structural_variation`、`somatic_structural_variation`、`constrained`、`regulatory`、`motif`、`chipseq`、`array_probe` +- **`biotype`**(字符串):按生物类型过滤(例如 `protein_coding`) +- **`variant_set`**(字符串):用于变异过滤的短集合名称 + +**示例:** 查找区域内所有基因和转录本: +``` +/overlap/region/human/7:140424943-140624564?feature=gene;feature=transcript +``` + +### `GET /overlap/id/{id}` +与 Ensembl 特征(基因、转录本等)重叠的特征。`feature` 参数选项同上。 + +### `GET /overlap/translation/{id}` +与翻译产物重叠的蛋白质特征。用于结构域注释。 + +- **`feature`**(字符串):`protein_feature`、`residue_overlap`、`translation_exon` +- **`type`**(字符串):按源数据库名称过滤(例如 `Pfam`、`Gene3D`、`CDD`、`Smart`、`SuperFamily`、`PANTHER`、`Prosite_patterns`、`PRINTS`、`MobiDBLite`) + +**示例:** 获取某个蛋白质的 Pfam 结构域注释: +``` +/overlap/translation/ENSP00000269305?feature=protein_feature;type=Pfam +``` + +**响应字段(protein_feature):** `type`(源数据库)、`id`(结构域登录号)、`description`、`start`(氨基酸起始位置)、`end`。 + +--- + +## 比较基因组学端点 + +### `GET /homology/id/{species}/{id}` +获取某个基因的同源物(直系同源/旁系同源)。 + +- **`type`**(字符串):`orthologues`、`paralogues`、`projections`、`all` +- **`target_species`**(字符串):限定到特定目标物种 +- **`target_taxon`**(整数):按 NCBI 分类单元 ID 限定 +- **`sequence`**(字符串):`none`、`cdna`、`protein` — 包含比对后的序列 + +### `GET /homology/symbol/{species}/{symbol}` +同上,但按基因符号而非 ID 查询。 + +### `GET /genetree/id/{id}` +通过 Ensembl Compara 树 ID 获取完整的基因树。 + +### `GET /genetree/member/id/{species}/{id}` +获取包含特定基因的基因树。 + +--- + +## 变异端点 + +### `GET /variation/{species}/{id}` +获取已知变异(通过 rsID 或 Ensembl 变异 ID)的详细信息。 + +**响应字段:** `name`(rsID)、`mappings[]`(包含 `location`、`allele_string`、`start`、`end`、`strand`)、`ancestral_allele`、`minor_allele`、`MAF`、`clinical_significance[]`、`source`。 + +### `GET /variant_recoder/{species}/{id}` +在不同格式(HGVS、VCF、SPDI、rsID)之间转换变异表示。 + +**响应字段:** `spdi[]`、`hgvsg[]`、`hgvsc[]`、`hgvsp[]`、`vcf_string[]`、`id[]`(rsID)。 + +--- + +## VEP(变异效应预测器)端点 + +### `GET /vep/{species}/region/{region}/{allele}` +预测基因组变异的影响。 + +区域格式:`CHR:START-END:STRAND`(例如 `9:21971147-21971147:1`)。Allele 为替代等位基因字符串。 + +### `GET /vep/{species}/id/{id}` +通过 rsID 预测变异影响。 + +### `GET /vep/{species}/hgvs/{hgvs_notation}` +通过 HGVS 命名法预测变异影响。 + +**插件参数(以查询参数形式附加):** + +- `AlphaMissense=1` — AlphaMissense 致病性预测 +- `Conservation=1` — PhyloP 保守性评分 +- `DosageSensitivity=1` — 单倍剂量不足/三倍剂量敏感 +- `LoF=loftee` — LOFTEE 功能缺失评估 +- `LOEUF=1` — 功能缺失观察值与期望值上限比值 +- `NMD=1` — 无义介导的 mRNA 降解预测 +- `UTRAnnotator=1` — 5'/3' UTR 变异注释 +- `mutfunc=1` — 功能影响预测 +- `IntAct=1` — 蛋白质相互作用影响 +- `MaveDB=1` — 多重分析评分 +- `OpenTargets=1` — Open Targets 遗传学数据 + +**关键响应字段:** `most_severe_consequence`、`transcript_consequences[]`(包含 `gene_symbol`、`transcript_id`、`consequence_terms[]`、`amino_acids`、`sift_prediction`、`sift_score`、`polyphen_prediction`、`polyphen_score`、`am_class`、`am_pathogenicity`、`conservation`、`lof`、`loeuf`)。 + +### `POST /vep/{species}/region` +批量 VEP:以类似 VCF 的格式发送 `{"variants": ["1 100 . A T . . ."]}`。 + +--- + +## 映射端点 + +### `GET /map/{species}/{asm_one}/{region}/{asm_two}` +在不同组装版本之间转换坐标(例如 GRCh37 → GRCh38)。 + +**示例:** +``` +/map/human/GRCh37/17:43044295-43125370/GRCh38 +``` + +**响应:** `mappings[]`,包含 `original` 和 `mapped` 坐标块。 + +### `GET /map/cdna/{id}/{region}` +将 cDNA 坐标映射到某个转录本的基因组坐标。 + +### `GET /map/cds/{id}/{region}` +将 CDS 坐标映射到基因组坐标。 + +### `GET /map/translation/{id}/{region}` +将蛋白质(氨基酸)位置映射到基因组坐标。 + +--- + +## 表型端点 + +### `GET /phenotype/gene/{species}/{gene}` +某个基因(按符号或 Ensembl ID)的表型注释。 + +### `GET /phenotype/region/{species}/{region}` +某个基因组区域中与表型相关的变异。 + +### `GET /phenotype/term/{species}/{term}` +查找与某个表型术语(本体论 ID 或描述字符串,例如 `coffee consumption`)相关的变异/基因。 + +--- + +## 连锁不平衡端点 + +### `GET /ld/{species}/{id}/{population_name}` +计算某个变异周围窗口内变异的 LD。 + +- **`window_size`**(整数):窗口大小,单位 kb(默认:500) +- **`r2`**(浮点数):最小 r² 阈值 +- **`d_prime`**(浮点数):最小 D' 阈值 + +### `GET /ld/{species}/pairwise/{id1}/{id2}` +两个特定变异之间的成对 LD。 + +### `GET /ld/{species}/region/{region}/{population_name}` +某个区域内所有变异对的 LD。 + +群体名称:例如 `1000GENOMES:phase_3:CEU`、`1000GENOMES:phase_3:YRI`。 + +--- + +## 信息端点 + +### `GET /info/assembly/{species}` +组装元数据(核型、顶层区域、坐标系统)。 + +### `GET /info/assembly/{species}/{region_name}` +特定染色体/区域的详细信息(长度、条带)。 + +### `GET /info/species` +列出 Ensembl REST API 中所有可用的物种。 + +### `GET /info/external_dbs/{species}` +列出所有可用于交叉引用的外部数据库名称。有助于查找正确的 `external_db` 参数值。 + +### `GET /info/biotypes/{species}` +列出某个物种的所有生物类型分类。 + +--- + +## 自定义查询的最佳实践 + +1. **始终设置 `Content-Type`**:JSON 响应设为 `application/json`,序列端点设为 `text/plain` 或 `text/x-fasta`。 +2. **在查询端点使用 `expand=1`**:通过一次调用获取子特征,避免为每个转录本/外显子分别发起请求。 +3. **优先使用批量端点**(`POST /lookup/id`、`POST /sequence/region`):当查询多个 ID 或区域时,减少 HTTP 往返次数。 +4. **使用交叉引用过滤前检查 `info/external_dbs`**:数据库名称字符串区分大小写,且并非总是显而易见(例如 `UniProt_gn` 而非 `UniProt`)。 +5. **区域大小限制**:`/overlap/region/` 端点对大多数特征类型有最大 5 Mb 的区域大小限制。较大区域请拆分。 +6. **当坐标为 GRCh37(hg19)组装时**,使用 `grch37.rest.ensembl.org` 作为基础 URL。大多数端点在该旧版服务器上支持相同的路径。 +7. **将响应保存到临时文件**——不要试图将大型 JSON 响应读入上下文。使用 `jq` 或 Python 一行命令提取特定字段。 diff --git a/scripts/ensembl_api.py b/scripts/ensembl_api.py new file mode 100644 index 0000000..1ab7225 --- /dev/null +++ b/scripts/ensembl_api.py @@ -0,0 +1,779 @@ +# 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()