生信数据库批量查询¶

白话：API 就是数据库留给程序员的"后门"。
网页版数据库是给人看的（点按钮、翻页面），
API 是给程序用的（发个网址请求，直接返回数据）。

类比：
- 网页查询 = 去银行柜台排队办业务（得一个个等）
- API 查询 = 用手机银行（自动化、批量、快速）

REST API 的核心：
- 用 URL 表示你要查什么（https://rest.uniprot.org/uniprotkb/P12345）
- 用 HTTP 方法表示你要做什么（GET=查询，POST=提交）
- 返回格式通常是 JSON 或 XML（结构化数据，Python 能直接解析）

白话：NCBI 旗下有几十个数据库（PubMed、SRA、GEO、GenBank...），
Entrez 是它们的"总管家"，让你用同一套接口查所有库。

类比：Entrez 就像购物 App 的"搜索框"，
虽然背后有服装、食品、电子等多个店铺（数据库），
但你只需要在一个搜索框里打字就能跨店搜索。

关键限制：
- 每秒最多 3 次请求（没有 API Key 时）
- 注册 API Key 后可以每秒 10 次
- 必须设置 email（NCBI 要求，不设会被封 IP）

批量查询：问数据库"给我这些东西的信息"（获取元数据/注释）
批量下载：把实际的大文件（FASTQ、基因组序列）拉到本地

举例：
- 批量查询：查 100 个 SRA 样本的物种、测序平台、文库类型
- 批量下载：把这 100 个样本的 FASTQ 文件全下载下来

工具对应：
- 批量查询 → Python + requests/Biopython
- 批量下载 → prefetch + fasterq-dump（SRA）、wget/curl（通用）

白话：元数据就是"关于数据的数据"，即描述信息。

举例（SRA 样本的元数据）：
- 样本来源：人肠道
- 疾病状态：T2D 患者 / 健康对照
- 测序平台：Illumina NovaSeq 6000
- 文库类型：WGS（全基因组鸟枪法）
- Read 长度：150bp paired-end
- 数据量：5 Gbp

为什么重要：
你做宏基因组分析，第一步永远是筛选合适的样本。
批量查元数据就是帮你快速"选品"——哪些样本符合你的分析条件。

# ============================================
# SRA 批量查询：查找 T2D 相关宏基因组数据
# ============================================
from Bio import Entrez  # Biopython 的 Entrez 模块
import time  # 用于控制请求频率

# 必须设置邮箱，NCBI 要求（不设会被封 IP）
Entrez.email = "yourname@example.com"
# 可选：注册 API Key 后每秒可发 10 次请求（默认只有 3 次）
# Entrez.api_key = "你的API_KEY"

# 第一步：搜索（esearch）—— 找到符合条件的记录 ID
handle = Entrez.esearch(
    db="sra",  # 指定搜索 SRA 数据库
    term="type 2 diabetes AND metagenome AND gut",  # 搜索词
    retmax=50,  # 最多返回 50 条结果
    usehistory="y"  # 使用历史服务器（处理大量结果时必须开）
)
results = Entrez.read(handle)  # 解析返回的 XML
handle.close()

print(f"找到 {results['Count']} 条相关记录")  # 打印总数
id_list = results['IdList']  # 获取记录 ID 列表

# 第二步：获取详情（efetch）—— 根据 ID 拿到具体信息
for sra_id in id_list[:5]:  # 先拿前 5 条看看
    handle = Entrez.efetch(
        db="sra",  # 从 SRA 库获取
        id=sra_id,  # 指定记录 ID
        rettype="full",  # 返回完整信息
        retmode="xml"  # XML 格式
    )
    record = handle.read()  # 读取原始 XML 文本
    handle.close()
    print(record[:500])  # 打印前 500 字符看看结构
    time.sleep(0.4)  # 每次请求间隔 0.4 秒，避免被封

# ============================================
# GEO 批量查询：查找 T2D 相关转录组数据集
# ============================================
from Bio import Entrez  # Entrez 统一检索
import time

Entrez.email = "yourname@example.com"

# 搜索 GEO DataSets 中与 T2D 相关的表达谱
handle = Entrez.esearch(
    db="gds",  # GEO DataSets 数据库
    term="type 2 diabetes[Title] AND Homo sapiens[Organism]",  # 标题含 T2D + 人类
    retmax=20  # 返回 20 条
)
results = Entrez.read(handle)
handle.close()

print(f"找到 {results['Count']} 个相关数据集")

# 获取数据集详细信息
for gds_id in results['IdList'][:3]:  # 取前 3 条
    handle = Entrez.esummary(
        db="gds",  # 从 GDS 获取摘要
        id=gds_id  # 指定 ID
    )
    summary = Entrez.read(handle)
    handle.close()
    # 打印数据集标题和描述
    for item in summary:
        print(f"标题: {item.get('title', 'N/A')}")
        print(f"摘要: {item.get('summary', 'N/A')[:200]}")
        print("---")
    time.sleep(0.4)  # 限速

# ============================================
# PubMed 批量查询：获取 T2D 肠道菌群最新文献
# ============================================
from Bio import Entrez
import time

Entrez.email = "yourname@example.com"

# 搜索 PubMed
handle = Entrez.esearch(
    db="pubmed",  # PubMed 数据库
    term="type 2 diabetes gut microbiome metagenomics",  # 搜索词
    retmax=10,  # 返回 10 条
    sort="relevance"  # 按相关度排序（也可以用 "pub_date" 按时间）
)
results = Entrez.read(handle)
handle.close()

print(f"找到 {results['Count']} 篇相关文献")

# 批量获取文献摘要
ids = ",".join(results['IdList'])  # 把 ID 列表拼成逗号分隔的字符串
handle = Entrez.efetch(
    db="pubmed",
    id=ids,  # 一次传入多个 ID，批量获取
    rettype="abstract",  # 获取摘要
    retmode="xml"  # XML 格式（方便解析）
)
records = Entrez.read(handle)
handle.close()

# 解析并打印每篇文献
for article in records['PubmedArticle']:
    # 层层解析 XML 结构拿到标题和摘要
    medline = article['MedlineCitation']
    title = medline['Article']['ArticleTitle']
    # 摘要可能不存在（综述类文章）
    abstract_list = medline['Article'].get('Abstract', {}).get('AbstractText', ['无摘要'])
    abstract = " ".join(str(a) for a in abstract_list)
    pmid = str(medline['PMID'])
    print(f"[PMID:{pmid}] {title}")
    print(f"  摘要前100字: {abstract[:100]}...")
    print()

# ============================================
# Ensembl REST API：批量获取基因信息
# ============================================
import requests  # HTTP 请求库
import time

# Ensembl REST API 基础地址
BASE_URL = "https://rest.ensembl.org"

def get_gene_info(gene_symbol, species="homo_sapiens"):
    """根据基因名获取 Ensembl 基因信息"""
    # lookup/symbol 端点：用基因符号查找
    url = f"{BASE_URL}/lookup/symbol/{species}/{gene_symbol}"
    headers = {"Content-Type": "application/json"}  # 要求返回 JSON
    response = requests.get(url, headers=headers)

    if response.status_code == 200:  # 200 表示成功
        return response.json()  # 解析 JSON 返回字典
    else:
        print(f"查询 {gene_symbol} 失败，状态码: {response.status_code}")
        return None

# T2D 相关基因列表（从你的差异分析结果中来）
t2d_genes = ["TCF7L2", "PPARG", "KCNJ11", "SLC30A8", "INS", "IRS1"]

# 批量查询
results = []
for gene in t2d_genes:
    info = get_gene_info(gene)
    if info:
        results.append({
            "symbol": gene,
            "ensembl_id": info.get("id"),  # Ensembl 基因 ID
            "chromosome": info.get("seq_region_name"),  # 染色体
            "start": info.get("start"),  # 起始位置
            "end": info.get("end"),  # 终止位置
            "biotype": info.get("biotype"),  # 生物类型（protein_coding 等）
            "description": info.get("description")  # 基因描述
        })
        print(f"{gene}: {info.get('id')} | Chr{info.get('seq_region_name')}:{info.get('start')}-{info.get('end')}")
    time.sleep(0.2)  # Ensembl 限制每秒 15 次请求

# 批量查询（POST 方式，更高效）
def batch_lookup_genes(gene_ids):
    """用 POST 方式批量查询多个 Ensembl ID"""
    url = f"{BASE_URL}/lookup/id"
    headers = {"Content-Type": "application/json", "Accept": "application/json"}
    # POST 请求体：传入 ID 列表
    data = {"ids": gene_ids}
    response = requests.post(url, headers=headers, json=data)

    if response.status_code == 200:
        return response.json()
    else:
        print(f"批量查询失败: {response.status_code}")
        return None

# 用上一步拿到的 Ensembl ID 做批量查询
ensembl_ids = [r["ensembl_id"] for r in results if r.get("ensembl_id")]
if ensembl_ids:
    batch_results = batch_lookup_genes(ensembl_ids)
    if batch_results:
        for eid, info in batch_results.items():
            print(f"{eid}: {info.get('display_name')} - {info.get('description', '')[:60]}")

# ============================================
# UniProt REST API：批量查询蛋白质信息
# ============================================
import requests
import time

# UniProt 新版 REST API 基础地址（2022年改版后）
BASE_URL = "https://rest.uniprot.org"

def search_uniprot(query, format="json", size=25):
    """搜索 UniProt 数据库"""
    url = f"{BASE_URL}/uniprotkb/search"
    params = {
        "query": query,  # 搜索语句
        "format": format,  # 返回格式：json/tsv/fasta
        "size": size,  # 每页结果数
        "fields": "accession,id,protein_name,gene_names,organism_name,length,go_f"
        # fields 指定返回哪些字段（减少数据量）
    }
    response = requests.get(url, params=params)

    if response.status_code == 200:
        return response.json()
    else:
        print(f"搜索失败: {response.status_code} - {response.text[:200]}")
        return None

# 搜索 T2D 相关的人类蛋白质
result = search_uniprot(
    query="type 2 diabetes AND organism_id:9606 AND reviewed:true",
    # organism_id:9606 是人类
    # reviewed:true 只搜 Swiss-Prot（人工审核的，质量高）
    size=10
)

if result and "results" in result:
    for entry in result["results"]:
        acc = entry.get("primaryAccession", "N/A")  # UniProt 登录号
        name = entry.get("proteinDescription", {}).get("recommendedName", {}).get("fullName", {}).get("value", "N/A")
        genes = entry.get("genes", [{}])
        gene_name = genes[0].get("geneName", {}).get("value", "N/A") if genes else "N/A"
        organism = entry.get("organism", {}).get("scientificName", "N/A")
        print(f"{acc} | {gene_name} | {name} | {organism}")

# 按 accession 批量获取（直接用 accession 列表）
def get_uniprot_entry(accession):
    """获取单条 UniProt 记录的完整信息"""
    url = f"{BASE_URL}/uniprotkb/{accession}.json"
    response = requests.get(url)
    if response.status_code == 200:
        return response.json()
    return None

# 示例：获取胰岛素受体的详细信息
insulin_receptor = get_uniprot_entry("P06213")
if insulin_receptor:
    print(f"蛋白名: {insulin_receptor['proteinDescription']['recommendedName']['fullName']['value']}")
    print(f"序列长度: {insulin_receptor['sequence']['length']} aa")
    # 打印 GO 功能注释
    for ref in insulin_receptor.get("uniProtKBCrossReferences", []):
        if ref.get("database") == "GO":
            props = {p["key"]: p["value"] for p in ref.get("properties", [])}
            print(f"  GO: {ref['id']} - {props.get('GoTerm', '')}")

# ============================================
# KEGG REST API：查询代谢通路
# ============================================
import requests
import time

# KEGG API 基础地址（http 和 https 均可用，官方文档常写 http）
BASE_URL = "https://rest.kegg.jp"

def kegg_find(database, query):
    """在 KEGG 指定数据库中搜索"""
    url = f"{BASE_URL}/find/{database}/{query}"
    response = requests.get(url)
    if response.status_code == 200:
        return response.text  # KEGG 返回的是纯文本，不是 JSON
    return None

def kegg_get(entry_id):
    """获取 KEGG 条目的详细信息"""
    url = f"{BASE_URL}/get/{entry_id}"
    response = requests.get(url)
    if response.status_code == 200:
        return response.text
    return None

def kegg_link(target_db, source_id):
    """查询两个数据库之间的映射关系"""
    url = f"{BASE_URL}/link/{target_db}/{source_id}"
    response = requests.get(url)
    if response.status_code == 200:
        return response.text
    return None

# 示例 1：搜索 T2D 相关通路
print("=== 搜索 T2D 相关通路 ===")
result = kegg_find("pathway", "type II diabetes")
if result:
    for line in result.strip().split("\n"):
        print(line)
        # 输出类似: path:hsa04930  Type II diabetes mellitus - Homo sapiens

# 示例 2：获取 T2D 通路的详细基因列表
print("\n=== T2D 通路包含的基因 ===")
genes_in_pathway = kegg_link("genes", "hsa04930")  # hsa04930 是人类 T2D 通路
if genes_in_pathway:
    gene_list = []
    for line in genes_in_pathway.strip().split("\n")[:10]:  # 只打印前 10 个
        parts = line.split("\t")
        if len(parts) == 2:
            gene_list.append(parts[1])
            print(parts[1])  # 输出类似: hsa:3630（INS 基因）

# 示例 3：查基因对应的所有通路
print("\n=== INS 基因参与的通路 ===")
ins_pathways = kegg_link("pathway", "hsa:3630")  # hsa:3630 是 INS（胰岛素）
if ins_pathways:
    for line in ins_pathways.strip().split("\n"):
        parts = line.split("\t")
        if len(parts) == 2:
            pathway_id = parts[1]
            # 获取通路名称
            info = kegg_find("pathway", pathway_id.replace("path:", ""))
            if info:
                print(f"  {info.strip().split(chr(10))[0]}")
            time.sleep(0.3)  # KEGG 限速

# ============================================
# STRING API：查询蛋白互作网络
# ============================================
import requests
import time

# STRING API 基础地址
BASE_URL = "https://string-db.org/api"
# STRING 使用 NCBI taxonomy ID 标识物种
SPECIES_TAXID = 9606  # 人类

def get_string_interactions(protein, required_score=400):
    """获取蛋白质的互作伙伴"""
    url = f"{BASE_URL}/json/network"
    params = {
        "identifiers": protein,  # 蛋白名（基因名也行）
        "species": SPECIES_TAXID,  # 物种
        "required_score": required_score,  # 最低置信度（0-1000）
        "limit": 20  # 最多返回 20 个互作
    }
    response = requests.get(url, params=params)
    if response.status_code == 200:
        return response.json()
    return None

def get_enrichment(proteins):
    """对蛋白列表做功能富集分析"""
    url = f"{BASE_URL}/json/enrichment"
    params = {
        "identifiers": "\r".join(proteins),  # 多个蛋白用回车分隔
        "species": SPECIES_TAXID
    }
    response = requests.get(url, params=params)
    if response.status_code == 200:
        return response.json()
    return None

# 查询 INS（胰岛素）的互作伙伴
print("=== INS 的蛋白互作伙伴 ===")
interactions = get_string_interactions("INS")
if interactions:
    for item in interactions[:10]:
        print(f"  {item['preferredName_A']} <-> {item['preferredName_B']} "
              f"(score: {item['score']:.3f})")

# 对 T2D 相关蛋白做富集分析
print("\n=== T2D 相关蛋白的功能富集 ===")
t2d_proteins = ["INS", "INSR", "IRS1", "IRS2", "PIK3CA", "AKT1", "SLC2A4"]
enrichment = get_enrichment(t2d_proteins)
if enrichment:
    # 只打印 KEGG 通路的富集结果
    kegg_terms = [e for e in enrichment if e.get("category") == "KEGG"]
    for term in kegg_terms[:5]:
        print(f"  {term['term']} | {term['description']} | p={term['p_value']:.2e}")

# ============================================
# Biopython Entrez 安装与基础配置
# ============================================
# 安装：conda install -c conda-forge biopython
# 或者：pip install biopython

from Bio import Entrez

# 必须设置！不设置 NCBI 会拒绝请求
Entrez.email = "yourname@example.com"

# 强烈建议注册 API Key（免费）：
# 去 https://www.ncbi.nlm.nih.gov/account/settings/ 注册
# 有 Key：每秒 10 次请求；没有：每秒 3 次
# Entrez.api_key = "你的key"

# ============================================
# 完整 Entrez 工作流：批量获取 T2D 宏基因组的 BioProject 信息
# ============================================
from Bio import Entrez
import time
import json

Entrez.email = "yourname@example.com"

# === 第一步：搜索 BioProject ===
print("Step 1: 搜索 T2D 相关的宏基因组项目...")
handle = Entrez.esearch(
    db="bioproject",  # BioProject 数据库
    term='("type 2 diabetes" OR "T2D") AND metagenome AND gut',
    retmax=100,  # 最多 100 条
    usehistory="y"  # 开启历史服务器（大结果集必开）
)
search_results = Entrez.read(handle)
handle.close()

count = int(search_results["Count"])
webenv = search_results["WebEnv"]  # 服务器端会话 ID
query_key = search_results["QueryKey"]  # 查询标识
print(f"  找到 {count} 个项目")

# === 第二步：分批获取（避免一次请求太多被封）===
print("\nStep 2: 分批获取项目详情...")
batch_size = 20  # 每批 20 条
all_projects = []

for start in range(0, min(count, 60), batch_size):  # 最多取 60 条
    print(f"  获取第 {start+1}-{start+batch_size} 条...")
    handle = Entrez.esummary(
        db="bioproject",
        query_key=query_key,  # 用之前搜索的 query_key
        WebEnv=webenv,  # 用之前的 WebEnv（避免重复搜索）
        retstart=start,  # 从第几条开始
        retmax=batch_size  # 每次取几条
    )
    summaries = Entrez.read(handle)
    handle.close()

    # 解析每条记录
    for doc in summaries.get("DocumentSummarySet", {}).get("DocumentSummary", []):
        project_info = {
            "uid": doc.attributes.get("uid", ""),
            "name": doc.get("Project_Name", ""),
            "title": doc.get("Project_Title", ""),
            "description": doc.get("Project_Description", "")[:200],
            "organism": doc.get("Organism_Name", ""),
            "registration_date": doc.get("Registration_Date", "")
        }
        all_projects.append(project_info)

    time.sleep(0.5)  # 批次间等待

# === 第三步：输出结果 ===
print(f"\nStep 3: 获取完成，共 {len(all_projects)} 个项目")
for i, proj in enumerate(all_projects[:5], 1):
    print(f"\n[{i}] {proj['title'][:80]}")
    print(f"    物种: {proj['organism']}")
    print(f"    日期: {proj['registration_date']}")

# 可以保存为 JSON 供后续分析
with open("t2d_bioprojects.json", "w", encoding="utf-8") as f:
    json.dump(all_projects, f, ensure_ascii=False, indent=2)
print("\n结果已保存到 t2d_bioprojects.json")

# ============================================
# 处理超大结果集（上万条记录）的正确姿势
# ============================================
from Bio import Entrez
import time

Entrez.email = "yourname@example.com"

# 搜索所有人类肠道宏基因组数据（可能有几万条）
handle = Entrez.esearch(
    db="sra",
    term="human gut metagenome[Organism] AND WGS[Strategy]",
    retmax=0,  # 先不返回 ID，只看总数
    usehistory="y"  # 关键！把结果存在 NCBI 服务器上
)
results = Entrez.read(handle)
handle.close()

total = int(results["Count"])
webenv = results["WebEnv"]
query_key = results["QueryKey"]
print(f"总共有 {total} 条记录，用 History Server 分批获取")

# 分批获取（每次 500 条，NCBI 建议不超过 500）
batch_size = 500
for start in range(0, min(total, 2000), batch_size):  # 示例只取前 2000 条
    handle = Entrez.efetch(
        db="sra",
        rettype="runinfo",  # SRA 特有：返回 CSV 格式的运行信息
        retmode="text",
        retstart=start,
        retmax=batch_size,
        webenv=webenv,
        query_key=query_key
    )
    data = handle.read()
    handle.close()

    # 保存每批数据（追加模式）
    mode = "w" if start == 0 else "a"
    with open("sra_runinfo.csv", mode) as f:
        if start == 0:
            f.write(data)  # 第一批包含表头
        else:
            # 跳过表头行（后续批次不需要重复写表头）
            lines = data.split("\n")
            f.write("\n".join(lines[1:]))

    print(f"  已获取 {start + batch_size}/{total}")
    time.sleep(1)  # 大量请求时间隔要更长

print("全部保存到 sra_runinfo.csv")

# conda 安装（推荐，最简单）
conda install -c bioconda sra-tools

# 验证安装
prefetch --version
fasterq-dump --version

# 配置缓存路径（避免下载到默认路径撑爆硬盘）
vdb-config --set /repository/user/main/public/root=/home/user/sra_cache

# ============================================
# 批量下载 SRA 数据：从 accession list 到 FASTQ
# ============================================

# 第一步：准备 accession 列表文件
# 文件内容每行一个 SRA 编号，例如：
# SRR12345678
# SRR12345679
# SRR12345680

# 第二步：批量预下载（prefetch）—— 先下载压缩的 .sra 文件
# -O 指定输出目录
# --max-size 50G 允许单个文件最大 50GB（默认 20GB 会跳过大文件）
prefetch --option-file sra_accession_list.txt \
         -O ./sra_downloads \
         --max-size 50G \
         --progress

# 第三步：批量转换为 FASTQ（fasterq-dump）
# 比旧版 fastq-dump 快 5-10 倍！
# --split-3：双端数据自动分成 _1.fastq 和 _2.fastq
# --threads 8：用 8 个线程加速
# --temp /tmp/sra_tmp：指定临时目录（需要大空间！约为输出的 2 倍）
for sra_file in ./sra_downloads/*/*.sra; do
    fasterq-dump "$sra_file" \
        --outdir ./fastq_output \
        --split-3 \
        --threads 8 \
        --temp /tmp/sra_tmp \
        --progress
done

# 第四步：压缩 FASTQ（节省 70% 空间）
pigz -p 8 ./fastq_output/*.fastq
# pigz 是多线程版 gzip，-p 8 用 8 线程

# ============================================
# Python 自动化 SRA 批量下载
# ============================================
import subprocess  # 调用系统命令
import os
from pathlib import Path  # 路径处理

# 配置
ACCESSION_FILE = "sra_accession_list.txt"  # accession 列表文件
OUTPUT_DIR = Path("./fastq_output")  # FASTQ 输出目录
SRA_CACHE = Path("./sra_cache")  # SRA 缓存目录
THREADS = 8  # 线程数

# 创建目录
OUTPUT_DIR.mkdir(exist_ok=True)
SRA_CACHE.mkdir(exist_ok=True)

# 读取 accession 列表
with open(ACCESSION_FILE) as f:
    accessions = [line.strip() for line in f if line.strip()]

print(f"共 {len(accessions)} 个样本待下载")

# 逐个下载并转换
for i, acc in enumerate(accessions, 1):
    print(f"\n[{i}/{len(accessions)}] 处理 {acc}...")

    # Step 1: prefetch 下载
    cmd_prefetch = [
        "prefetch", acc,
        "-O", str(SRA_CACHE),
        "--max-size", "50G"
    ]
    result = subprocess.run(cmd_prefetch, capture_output=True, text=True)
    if result.returncode != 0:
        print(f"  prefetch 失败: {result.stderr[:200]}")
        continue  # 跳过失败的，继续下一个

    # Step 2: fasterq-dump 转换
    sra_path = SRA_CACHE / acc / f"{acc}.sra"
    cmd_fasterq = [
        "fasterq-dump", str(sra_path),
        "--outdir", str(OUTPUT_DIR),
        "--split-3",
        "--threads", str(THREADS)
    ]
    result = subprocess.run(cmd_fasterq, capture_output=True, text=True)
    if result.returncode != 0:
        print(f"  fasterq-dump 失败: {result.stderr[:200]}")
        continue

    # Step 3: 压缩 FASTQ
    for fq in OUTPUT_DIR.glob(f"{acc}*.fastq"):
        subprocess.run(["gzip", str(fq)])

    print(f"  完成! FASTQ 已保存并压缩")

    # 清理缓存（释放空间）
    sra_file = SRA_CACHE / acc / f"{acc}.sra"
    if sra_file.exists():
        sra_file.unlink()

print(f"\n全部完成! FASTQ 文件在: {OUTPUT_DIR}")

# ============================================
# UniProt 批量查询完整示例
# 场景：从宏基因组功能注释结果中，批量查询蛋白功能
# ============================================
import requests
import time
import csv
from io import StringIO

BASE_URL = "https://rest.uniprot.org"

# === 方法 1：ID Mapping（最常用）===
# 场景：你有一批 UniProt ID，想批量获取注释信息
def uniprot_id_mapping(from_db, to_db, ids):
    """
    UniProt ID 映射服务
    from_db: 来源数据库（如 "UniProtKB_AC-ID"）
    to_db: 目标数据库（如 "UniProtKB"）
    ids: ID 列表
    """
    # 第一步：提交映射任务
    url = f"{BASE_URL}/idmapping/run"
    data = {
        "from": from_db,
        "to": to_db,
        "ids": ",".join(ids)  # 逗号分隔
    }
    response = requests.post(url, data=data)
    job_id = response.json()["jobId"]
    print(f"  任务已提交，Job ID: {job_id}")

    # 第二步：轮询等待结果
    status_url = f"{BASE_URL}/idmapping/status/{job_id}"
    while True:
        status = requests.get(status_url).json()
        if "results" in status or "jobStatus" not in status:
            break
        if status["jobStatus"] == "RUNNING":
            time.sleep(2)  # 每 2 秒检查一次
            print("  等待中...")
        else:
            print(f"  任务状态: {status['jobStatus']}")
            break

    # 第三步：获取结果
    result_url = f"{BASE_URL}/idmapping/results/{job_id}"
    response = requests.get(result_url)
    return response.json()

# 示例：批量查询蛋白
protein_ids = ["P06213", "P01308", "P35568", "Q14524"]  # INSR, INS, IRS1, SCN5A
print("=== ID Mapping 批量查询 ===")
results = uniprot_id_mapping("UniProtKB_AC-ID", "UniProtKB", protein_ids)

if "results" in results:
    for item in results["results"]:
        entry = item.get("to", {})
        acc = entry.get("primaryAccession", "N/A")
        name = entry.get("proteinDescription", {}).get("recommendedName", {}).get("fullName", {}).get("value", "N/A")
        print(f"  {acc}: {name}")

# === 方法 2：流式下载大量结果 ===
def stream_uniprot_search(query, fields, format="tsv"):
    """
    流式搜索，适合下载大量结果
    返回 TSV 格式，可以直接导入 pandas
    """
    url = f"{BASE_URL}/uniprotkb/stream"
    params = {
        "query": query,
        "fields": fields,  # 逗号分隔的字段名
        "format": format
    }
    response = requests.get(url, params=params, stream=True)

    if response.status_code == 200:
        return response.text
    else:
        print(f"查询失败: {response.status_code}")
        return None

# 下载所有与短链脂肪酸（SCFA）代谢相关的细菌蛋白
# （SCFA 在 T2D 肠道菌群研究中非常重要）
print("\n=== 流式下载 SCFA 相关细菌蛋白 ===")
tsv_data = stream_uniprot_search(
    query="short-chain fatty acid AND taxonomy_id:2 AND reviewed:true",
    # taxonomy_id:2 是细菌域
    fields="accession,id,protein_name,gene_names,organism_name,ec,go_f",
    format="tsv"
)

if tsv_data:
    # 解析 TSV
    reader = csv.DictReader(StringIO(tsv_data), delimiter="\t")
    rows = list(reader)
    print(f"  共找到 {len(rows)} 条结果")
    for row in rows[:5]:
        print(f"  {row.get('Entry', '')} | {row.get('Gene Names', '')} | "
              f"{row.get('Protein names', '')[:50]}")

# ============================================
# Ensembl REST API 进阶用法
# ============================================
import requests
import time
import json

BASE_URL = "https://rest.ensembl.org"
HEADERS = {"Content-Type": "application/json"}

# === 获取基因的所有转录本 ===
def get_transcripts(gene_id):
    """获取一个基因的所有转录本信息"""
    url = f"{BASE_URL}/lookup/id/{gene_id}?expand=1"
    # expand=1 会展开子对象（包括转录本列表）
    response = requests.get(url, headers=HEADERS)
    if response.status_code == 200:
        data = response.json()
        transcripts = data.get("Transcript", [])
        return transcripts
    return []

# 查询 TCF7L2（T2D 最强关联基因）的转录本
print("=== TCF7L2 的转录本 ===")
# 先拿到 Ensembl ID
lookup = requests.get(
    f"{BASE_URL}/lookup/symbol/homo_sapiens/TCF7L2",
    headers=HEADERS
).json()
tcf7l2_id = lookup["id"]
print(f"TCF7L2 Ensembl ID: {tcf7l2_id}")

transcripts = get_transcripts(tcf7l2_id)
for t in transcripts[:5]:
    print(f"  {t['id']} | {t['biotype']} | 长度: {t['length']} bp")

# === 获取序列 ===
def get_sequence(seq_id, seq_type="genomic"):
    """
    获取序列
    seq_type: genomic/cds/cdna/protein
    """
    url = f"{BASE_URL}/sequence/id/{seq_id}?type={seq_type}"
    headers = {"Content-Type": "text/plain"}  # 纯文本返回序列
    response = requests.get(url, headers=headers)
    if response.status_code == 200:
        return response.text
    return None

# 获取 TCF7L2 蛋白序列
print(f"\n=== TCF7L2 蛋白序列（前 100 aa）===")
protein_seq = get_sequence(tcf7l2_id, "protein")
if protein_seq:
    print(protein_seq[:100] + "...")

# === 变异信息查询 ===
def get_variants(gene_symbol, species="human"):
    """获取基因区域的已知变异"""
    # 先获取基因坐标
    lookup_url = f"{BASE_URL}/lookup/symbol/{species}/{gene_symbol}"
    gene_info = requests.get(lookup_url, headers=HEADERS).json()

    chrom = gene_info["seq_region_name"]
    start = gene_info["start"]
    end = gene_info["end"]

    # 用 overlap 端点获取变异
    # 注意：区域太大会超时，限制在基因范围内
    region = f"{chrom}:{start}-{min(end, start + 50000)}"  # 最多查 50kb
    url = f"{BASE_URL}/overlap/region/human/{region}?feature=variation"
    response = requests.get(url, headers=HEADERS)

    if response.status_code == 200:
        return response.json()
    return []

print(f"\n=== KCNJ11 基因区域的变异（前 5 个）===")
variants = get_variants("KCNJ11")
if variants:
    print(f"  共找到 {len(variants)} 个变异")
    for v in variants[:5]:
        print(f"  {v.get('id', 'N/A')} | {v.get('consequence_type', 'N/A')} | "
              f"位置: {v.get('start', '')}")
    time.sleep(0.2)

# === 比较基因组学：查同源基因 ===
def get_homologs(gene_id, target_species=None):
    """获取同源基因"""
    url = f"{BASE_URL}/homology/id/{gene_id}?type=orthologues"
    if target_species:
        url += f"&target_species={target_species}"
    response = requests.get(url, headers=HEADERS)
    if response.status_code == 200:
        return response.json()
    return None

print(f"\n=== TCF7L2 在小鼠中的同源基因 ===")
homologs = get_homologs(tcf7l2_id, "mus_musculus")
if homologs:
    for h in homologs.get("data", [{}])[0].get("homologies", [])[:3]:
        target = h.get("target", {})
        print(f"  {target.get('id')} | {target.get('species')} | "
              f"相似度: {h.get('dn_ds', 'N/A')}")

# ============================================
# 实战：把 Kraken2 差异菌的代谢功能查出来
# ============================================
import requests
import time

# 假设你差异分析发现这些菌在 T2D 患者中显著变化
differential_bacteria = [
    "Faecalibacterium prausnitzii",  # 产丁酸菌，T2D 中减少
    "Akkermansia muciniphila",  # 粘液降解菌，T2D 中减少
    "Roseburia intestinalis",  # 产丁酸菌
    "Bacteroides vulgatus",  # 拟杆菌
]

# 在 KEGG 中查这些菌的代谢模块
KEGG_BASE = "https://rest.kegg.jp"

for bacterium in differential_bacteria:
    print(f"\n{'='*50}")
    print(f"查询: {bacterium}")

    # 在 KEGG 的 organism 库中搜索
    result = requests.get(f"{KEGG_BASE}/find/genome/{bacterium}").text

    if result.strip():
        # 解析第一个结果获取 organism code
        first_line = result.strip().split("\n")[0]
        org_code = first_line.split("\t")[0].replace("gn:", "")
        print(f"  KEGG organism: {org_code}")

        # 查该菌的代谢通路
        pathways = requests.get(f"{KEGG_BASE}/list/pathway/{org_code}").text
        pathway_lines = pathways.strip().split("\n")
        print(f"  共有 {len(pathway_lines)} 条代谢通路")

        # 过滤出与碳水化合物/脂肪酸代谢相关的通路（与 T2D 相关）
        keywords = ["carbohydrate", "fatty acid", "butyrate", "propanoate", "insulin"]
        relevant = [l for l in pathway_lines 
                    if any(kw in l.lower() for kw in keywords)]
        if relevant:
            print(f"  T2D 相关通路:")
            for p in relevant[:3]:
                print(f"    {p}")
    else:
        print(f"  KEGG 中未找到该菌")

    time.sleep(0.5)  # 限速

错误信息：requests.exceptions.HTTPError: 429 Client Error: Too Many Requests
原因：请求太频繁，超过 API 的速率限制
解决：
  1. 每次请求之间加 time.sleep(0.5)
  2. NCBI：注册 API Key（https://www.ncbi.nlm.nih.gov/account/settings/）
  3. Ensembl：每秒不超过 15 次，批量用 POST 方式
  4. UniProt：使用 stream 端点替代分页查询

现象：Entrez.read() 返回的 IdList 为空，但网页搜索有结果
原因：
  1. 搜索词语法错误（布尔运算符必须大写：AND/OR/NOT）
  2. 字段限定符拼写错误（如 [Titl] 写成了 [Title] 的缩写形式）
  3. 没设置 Entrez.email（被静默拒绝）
解决：
  - 先在 NCBI 网页上测试搜索词确认有结果
  - 检查布尔运算符大小写：term="T2D[Title] AND metagenome[All Fields]"
  - 确保设置了 Entrez.email

错误信息：400 Client Error: Bad Request
原因：
  1. 新版 API（2022 年后）URL 格式变了
  2. fields 参数中的字段名拼错了
  3. query 语法错误
解决：
  - 旧版 URL（已弃用）：https://www.uniprot.org/uniprot/
  - 新版 URL（用这个）：https://rest.uniprot.org/uniprotkb/
  - 查 UniProt 官方字段名文档：https://www.uniprot.org/help/return_fields
  - 用 Postman 或 curl 先测试单条请求

错误信息：requests.exceptions.ReadTimeout 或 504 Gateway Timeout
原因：请求的数据范围太大（如查一整条染色体的变异）
解决：
  - 缩小查询范围（区域不超过 5Mb）
  - 加 timeout 参数：requests.get(url, timeout=30)
  - 用 POST /lookup/id 批量查询替代多次 GET
  - 对于大规模数据，考虑下载 Ensembl 的 FTP 文件

错误信息：disk is full / no space left on device
原因：fasterq-dump 需要临时空间，约为输出文件的 2-3 倍
解决：
  1. 指定临时目录到大硬盘：fasterq-dump --temp /大硬盘/tmp
  2. 逐个处理，处理完一个压缩一个再下一个
  3. 提前算空间：1个 10GB 的 SRA 约产生 30GB FASTQ + 30GB 临时空间
  4. 下载后立即 pigz 压缩（节省 70%）

错误信息：403 Forbidden
原因：
  1. KEGG 限制商业用途（学术免费，商业需要授权）
  2. 请求频率过高（虽然没有明确文档，一般每秒 3-5 次安全）
  3. 某些端点确实需要付费（如 KEGG MEDICUS）
解决：
  - 确保用于学术研究
  - 加大请求间隔：time.sleep(0.5)
  - 用 https://rest.kegg.jp（http 和 https 均可）
  - 大规模需求考虑下载 KEGG FTP 数据（需授权）

=== NCBI Entrez（通过 Biopython 调用）===
Entrez.esearch(db, term)           # 搜索
Entrez.efetch(db, id)              # 获取完整记录
Entrez.esummary(db, id)            # 获取摘要
Entrez.elink(dbfrom, db, id)       # 跨库链接
支持的 db：pubmed, sra, gds, bioproject, biosample, gene, protein, nucleotide...

=== Ensembl REST ===
GET /lookup/symbol/{species}/{symbol}     # 按基因名查
GET /lookup/id/{id}                        # 按 Ensembl ID 查
POST /lookup/id   (body: {"ids":[...]})   # 批量查
GET /sequence/id/{id}?type=protein        # 获取序列
GET /overlap/region/{species}/{region}    # 区域内的注释/变异
GET /homology/id/{id}                     # 同源基因
GET /xrefs/id/{id}                        # 交叉引用（链接到其他数据库）

=== UniProt REST ===
GET  /uniprotkb/search?query=...&fields=...&format=json    # 搜索
GET  /uniprotkb/{accession}.json                           # 单条获取
GET  /uniprotkb/stream?query=...&format=tsv                # 流式下载
POST /idmapping/run                                         # ID 映射（提交）
GET  /idmapping/status/{jobId}                             # 映射状态
GET  /idmapping/results/{jobId}                            # 映射结果

=== KEGG REST ===
GET /find/{database}/{query}       # 搜索
GET /get/{entry_id}                # 获取条目详情
GET /link/{target_db}/{source_id}  # 数据库间映射
GET /list/{database}               # 列出所有条目
GET /get/{pathway_id}/image        # 获取通路图（PNG）

=== STRING ===
GET /api/json/network?identifiers=...&species=...        # 获取互作网络
GET /api/json/enrichment?identifiers=...&species=...     # 功能富集
GET /api/image/network?identifiers=...&species=...       # 获取网络图（PNG）
GET /api/json/interaction_partners?identifiers=...       # 互作伙伴列表