chore: import zh skill ensembl-database
This commit is contained in:
@@ -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 清单整理
|
||||
- 原作者、版权和许可证信息以上游仓库为准
|
||||
@@ -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` 错误(例如使用指数退避)。
|
||||
@@ -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 一行命令提取特定字段。
|
||||
@@ -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()
|
||||
Reference in New Issue
Block a user