ENSEMBL BioMart批量数据导出：REST API与biomaRt

1
#!/usr/bin/env python3
2
"""
3
fetch_gene_info.py - 通过 BioMart REST API 获取基因基础注释
4
Debian 12 实测 2025-12-05
5
"""
6

7
import requests
8
import json
9
import time
10

11
BASE_URL = "https://rest.ensembl.org"
12

13
def query_biomart(dataset, filters, attributes):
14
    """
15
    通用 BioMart 查询函数
16

17
    参数:
18
        dataset: 数据集名，如 "hsapiens_gene_ensembl"
19
        filters: 字典，过滤条件 {"filter_name": "value"}
20
        attributes: 列表，输出属性 ["attr1", "attr2"]
21

22
    返回:
23
        list of dict，查询结果
24
    """
25
    url = f"{BASE_URL}/biomart/martservice/results"
26

27
    # 构造 XML payload（BioMart 的 REST API 用 XML 描述查询）
28
    filter_xml = "".join(
29
        f'<Filter name="{k}" value="{v}"/>'
30
        for k, v in filters.items()
31
    )
32
    attr_xml = "".join(
33
        f'<Attribute name="{a}"/>'
34
        for a in attributes
35
    )
36

37
    payload = f"""<?xml version="1.0" encoding="UTF-8"?>
38
<!DOCTYPE Query>
39
<Query virtualSchemaName="default" formatter="TSV"
40
       header="1" uniqueRows="1"
41
       datasetConfigVersion="0.6">
42
    <Dataset name="{dataset}" interface="default">
43
        {filter_xml}
44
        {attr_xml}
45
    </Dataset>
46
</Query>"""
47

48
    headers = {"Content-Type": "application/x-www-form-urlencoded"}
49

50
    resp = requests.post(
51
        url,
52
        data={"query": payload},
53
        headers=headers,
54
        timeout=120
55
    )
56

57
    if resp.status_code != 200:
58
        raise Exception(f"API error {resp.status_code}: {resp.text}")
59

60
    # TSV → 字典列表
61
    lines = resp.text.strip().split("\n")
62
    if len(lines) < 2:
63
        return []
64

65
    header = lines[0].split("\t")
66
    results = []
67
    for line in lines[1:]:
68
        values = line.split("\t")
69
        results.append(dict(zip(header, values)))
70

71
    return results
72

73

74
# ===== 示例1：获取TP53及其靶基因的基础信息 =====
75
gene_list = ["TP53", "MDM2", "CDKN1A", "BAX", "BCL2", "GADD45A"]
76

77
filters = {"hgnc_symbol": ",".join(gene_list)}
78
attributes = [
79
    "ensembl_gene_id",
80
    "external_gene_name",
81
    "chromosome_name",
82
    "start_position",
83
    "end_position",
84
    "strand",
85
    "gene_biotype",
86
    "description"
87
]
88

89
print("=== 查询基因基础信息 ===")
90
results = query_biomart("hsapiens_gene_ensembl", filters, attributes)
91

92
for row in results:
93
    gene = row.get("Gene name", "N/A")
94
    chrom = row.get("Chromosome/scaffold name", "N/A")
95
    start = row.get("Gene start (bp)", "N/A")
96
    biotype = row.get("Gene type", "N/A")
97
    desc = row.get("Gene description", "N/A")
98

99
    # 计算基因长度
100
    if start != "N/A" and row.get("Gene end (bp)", "N/A") != "N/A":
101
        length = int(row["Gene end (bp)"]) - int(start) + 1
102
        length_str = f"{length:,}bp"
103
    else:
104
        length_str = "N/A"
105

106
    print(f"  {gene}: {chrom}:{start}, {biotype}, {length_str}")
107
    print(f"    {desc[:80]}...")
108

109
time.sleep(0.5)  # 礼貌性延迟，别打爆 API

1
# ===== 示例2：获取基因的GO注释 =====
2

3
filters = {"hgnc_symbol": ",".join(gene_list)}
4
attributes = [
5
    "ensembl_gene_id",
6
    "external_gene_name",
7
    "go_id",
8
    "name_1006",       # GO term 名称
9
    "namespace_1003",  # GO 类别（biological_process/等）
10
    "go_linkage_type"  # IEA/IDA/IMP...（证据类型）
11
]
12

13
print("\n=== 查询 GO 注释 ===")
14
results = query_biomart("hsapiens_gene_ensembl", filters, attributes)
15

16
# 按基因分组整理
17
from collections import defaultdict
18
go_by_gene = defaultdict(list)
19

20
for row in results:
21
    gene = row.get("Gene name", "N/A")
22
    go_id = row.get("GO term accession", "")
23
    go_name = row.get("GO term name", "")
24
    go_ns = row.get("GO domain", "")
25
    evidence = row.get("GO term evidence code", "")
26

27
    if go_id:
28
        go_by_gene[gene].append({
29
            "id": go_id,
30
            "name": go_name,
31
            "namespace": go_ns,
32
            "evidence": evidence
33
        })
34

35
for gene, go_terms in go_by_gene.items():
36
    print(f"\n  {gene} ({len(go_terms)} GO terms):")
37
    # 按 namespace 分组
38
    by_ns = defaultdict(list)
39
    for t in go_terms:
40
        by_ns[t["namespace"]].append(t)
41

42
    for ns, terms in sorted(by_ns.items()):
43
        print(f"    [{ns}]:")
44
        for t in terms[:3]:  # 每类只显示3个
45
            print(f"      {t['id']} {t['name']} ({t['evidence']})")

1
# ===== 示例3：人类→小鼠同源基因 =====
2

3
filters = {"hgnc_symbol": ",".join(gene_list)}
4
attributes = [
5
    "ensembl_gene_id",
6
    "external_gene_name",
7
    "mmusculus_homolog_ensembl_gene",
8
    "mmusculus_homolog_associated_gene_name",
9
    "mmusculus_homolog_orthology_type",
10
    "mmusculus_homolog_perc_id_r1"
11
]
12

13
print("\n=== 查询人类→小鼠同源基因 ===")
14
results = query_biomart("hsapiens_gene_ensembl", filters, attributes)
15

16
for row in results:
17
    hgnc = row.get("Gene name", "N/A")
18
    mouse_id = row.get("Mouse gene stable ID", "N/A")
19
    mouse_name = row.get("Mouse gene name", "N/A")
20
    ortho_type = row.get("Mouse orthology type", "N/A")
21
    identity = row.get("Mouse %id. query gene identical to target gene", "N/A")
22

23
    if mouse_id:
24
        print(f"  {hgnc} → {mouse_name} ({ortho_type}, {identity}% identity)")

1
#!/usr/bin/env Rscript
2
# fetch_genes.R - 用 biomaRt 批量获取基因注释
3
# Debian 12 实测 2025-12-05
4

5
library(biomaRt)
6

7
# ===== 1. 连接 BioMart =====
8
# 列出可用的 mart
9
listEnsembl()
10

11
# 选择 ENSEMBL Genes mart
12
ensembl <- useEnsembl(
13
    biomart = "genes",
14
    dataset = "hsapiens_gene_ensembl"
15
)
16

17
# 确认连接成功
18
ensembl
19
# Object of class 'Mart':
20
#   Using the ENSEMBL Genes dataset
21
#   Using the hsapiens_gene_ensembl dataset

1
# ===== 2. 按基因名查询基础信息 =====
2
gene_list <- c("TP53", "MDM2", "CDKN1A", "BAX", "BCL2", "GADD45A")
3

4
# 查看可用属性
5
attrs <- listAttributes(ensembl)
6
head(attrs[, c("name", "description")], 10)
7

8
# 执行查询
9
gene_info <- getBM(
10
    attributes = c(
11
        "ensembl_gene_id",
12
        "external_gene_name",
13
        "chromosome_name",
14
        "start_position",
15
        "end_position",
16
        "strand",
17
        "gene_biotype",
18
        "description"
19
    ),
20
    filters = "hgnc_symbol",
21
    values = gene_list,
22
    mart = ensembl
23
)
24

25
# 计算基因长度
26
gene_info$gene_length_bp <- gene_info$end_position - gene_info$start_position + 1
27

28
print(gene_info[, c("external_gene_name", "chromosome_name",
29
                     "gene_biotype", "gene_length_bp")])
30

31
#   external_gene_name chromosome_name   gene_biotype gene_length_bp
32
# 1             CDKN1A               6 protein_coding          50481
33
# 2                BAX              19 protein_coding           7096
34
# 3               BCL2              18 protein_coding         198810
35
# 4             GADD45A               1 protein_coding           3017
36
# 5               MDM2              12 protein_coding          36675
37
# 6               TP53              17 protein_coding          26383

1
# ===== 3. 批量获取 GO 注释 =====
2
go_terms <- getBM(
3
    attributes = c(
4
        "ensembl_gene_id",
5
        "external_gene_name",
6
        "go_id",
7
        "name_1006",       # GO term 名称
8
        "namespace_1003",  # GO 类别
9
        "go_linkage_type"  # 证据代码
10
    ),
11
    filters = "hgnc_symbol",
12
    values = gene_list,
13
    mart = ensembl
14
)
15

16
head(go_terms, 10)
17

18
# 统计每个基因的 GO 数量和证据类型分布
19
library(dplyr)
20

21
go_summary <- go_terms %>%
22
    group_by(external_gene_name, namespace_1003) %>%
23
    summarise(
24
        n_terms = n(),
25
        experimental = sum(go_linkage_type %in% c("IDA", "IMP", "IGI", "IPI")),
26
        computational = sum(go_linkage_type %in% c("IEA", "ISS")),
27
        .groups = "drop"
28
    )
29

30
print(go_summary)

1
# ===== 4. 获取基因序列 =====
2
# 查看可用的序列类型
3
listFilters(ensembl)[grep("seq", listFilters(ensembl)$name), ]
4

5
# 获取 cDNA 序列
6
cdna_seqs <- getSequence(
7
    id = gene_info$ensembl_gene_id,
8
    type = "ensembl_gene_id",
9
    seqType = "cdna",
10
    mart = ensembl
11
)
12

13
# 写入 FASTA
14
write.fasta <- function(ids, seqs, file) {
15
    con <- file(file, "w")
16
    for (i in seq_along(ids)) {
17
        writeLines(paste0(">", ids[i]), con)
18
        writeLines(seqs[i], con)
19
    }
20
    close(con)
21
}
22

23
# 只在有序列数据时才写入
24
if (nrow(cdna_seqs) > 0) {
25
    write.fasta(cdna_seqs$ensembl_gene_id, cdna_seqs$cdna,
26
                "gene_cdna.fasta")
27
    cat(sprintf("Wrote %d cDNA sequences to gene_cdna.fasta\n",
28
                nrow(cdna_seqs)))
29
}
30

31
# 获取蛋白序列
32
protein_seqs <- getSequence(
33
    id = gene_info$ensembl_gene_id,
34
    type = "ensembl_gene_id",
35
    seqType = "peptide",
36
    mart = ensembl
37
)

1
# ===== 5. 人类→小鼠同源基因 =====
2
mouse_homologs <- getBM(
3
    attributes = c(
4
        "ensembl_gene_id",
5
        "external_gene_name",
6

7
        # 小鼠同源基因属性
8
        "mmusculus_homolog_ensembl_gene",
9
        "mmusculus_homolog_associated_gene_name",
10
        "mmusculus_homolog_orthology_type",
11
        "mmusculus_homolog_perc_id_r1",
12
        "mmusculus_homolog_perc_id",
13

14
        # 同源关系类型
15
        "mmusculus_homolog_orthology_confidence"
16
    ),
17
    filters = "hgnc_symbol",
18
    values = gene_list,
19
    mart = ensembl
20
)
21

22
print(mouse_homologs)

1
# ===== 6. ID 批量转换 =====
2
# 模拟一批 Ensembl ID
3
ensembl_ids <- c(
4
    "ENSG00000141510",  # TP53
5
    "ENSG00000135679",  # MDM2
6
    "ENSG00000124762",  # CDKN1A
7
    "ENSG00000087088",  # BAX
8
    "ENSG00000171791",  # BCL2
9
    "ENSG00000116717"   # GADD45A
10
)
11

12
id_map <- getBM(
13
    attributes = c(
14
        "ensembl_gene_id",
15
        "external_gene_name",
16
        "entrezgene_id",
17
        "uniprotswissprot",
18
        "hgnc_id"
19
    ),
20
    filters = "ensembl_gene_id",
21
    values = ensembl_ids,
22
    mart = ensembl
23
)
24

25
print(id_map)
26

27
#    ensembl_gene_id external_gene_name entrezgene_id uniprotswissprot hgnc_id
28
# 1 ENSG00000087088                BAX           581           Q07812   990
29
# 2 ENSG00000116717           GADD45A          1647           P24522  4095
30
# 3 ENSG00000124762            CDKN1A          1026           P38936  1784
31
# 4 ENSG00000135679              MDM2          4193               NA  6973
32
# 5 ENSG00000141510              TP53          7157           P04637 11998
33
# 6 ENSG00000171791              BCL2           596           P10415   990

1
genes_in_region <- getBM(
2
    attributes = c(
3
        "ensembl_gene_id",
4
        "external_gene_name",
5
        "chromosome_name",
6
        "start_position",
7
        "end_position",
8
        "gene_biotype"
9
    ),
10
    filters = c("chromosome_name", "start", "end"),
11
    values = list(
12
        chromosome_name = "1",
13
        start = "1000000",
14
        end = "2000000"
15
    ),
16
    mart = ensembl
17
)
18

19
cat(sprintf("Found %d genes in chr1:1000000-2000000\n", nrow(genes_in_region)))

1
import time
2

3
# 每次请求后至少等0.5秒
4
time.sleep(0.5)
5

6
# 或者用重试逻辑
7
import time
8
max_retries = 3
9
for attempt in range(max_retries):
10
    resp = requests.post(url, ...)
11
    if resp.status_code == 429:
12
        wait = 2 ** attempt  # 指数退避
13
        print(f"Rate limited, waiting {wait}s...")
14
        time.sleep(wait)
15
    elif resp.status_code == 200:
16
        break