生信数据库批量查询¶

彭文强 | 2026届 | 生信分析工程师求职

1. 一句话说明¶

生信数据库批量查询就是用编程方式（API/命令行）一次性从 NCBI、UniProt、Ensembl 等公共数据库中获取成百上千条数据，替代手动逐条搜索下载的低效操作。

2. 为什么要学¶

你目前的情况¶

你的 T2D 项目中已经用过 prefetch + fasterq-dump 从 NCBI SRA 下载宏基因组测序数据，但那只是"知道下载命令"。系统学习批量查询能让你：

场景	手动查（网页）	批量查（API）
下载 100 个 SRA 样本的元数据	点 100 次，复制粘贴半天	一个脚本 30 秒搞定
查 T2D 相关基因在 KEGG 中的通路	一个个搜，容易漏	批量提交基因列表，全部返回
获取差异菌对应蛋白的功能注释	UniProt 网页一条条查	REST API 一次查完
做文献综述（T2D + gut microbiome）	PubMed 手动翻页	Entrez 自动拉取摘要
查宏基因组 MAG 中基因的同源信息	Ensembl 逐个搜	REST API 批量获取

白话总结：手动查数据库像去超市一个个挑菜，API 批量查像直接给超市发采购单，效率差 100 倍。面试中说"我能用 Python 批量查询各大数据库做自动化分析"，直接加分。

3. 核心概念白话版¶

3.1 REST API（表述性状态转移应用编程接口）¶

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

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

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

3.2 Entrez（NCBI 统一检索系统）¶

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

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

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

3.3 批量下载 vs 批量查询¶

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

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

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

3.4 元数据（Metadata）¶

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

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

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

4. 主要数据库详解¶

4.1 NCBI SRA（Sequence Read Archive）¶

是什么：全球最大的高通量测序原始数据仓库，存储了几千万个样本的测序文件。

能查什么： - 测序原始数据（FASTQ） - 样本元数据（物种、疾病、测序平台等） - 项目信息（BioProject）、样本信息（BioSample）

API 用法：通过 Entrez 的 esearch + efetch 查询

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

# 必须设置邮箱，NCBI 要求（不设会被封 IP）
Entrez.email = "pengwenqiang@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 秒，避免被封

4.2 NCBI GEO（Gene Expression Omnibus）¶

是什么：基因表达数据仓库，存储芯片和 RNA-seq 的表达矩阵和实验设计。

能查什么： - 基因表达矩阵（做差异分析用的） - 实验设计（哪些是实验组、对照组） - 数据系列（GSE）、样本（GSM）、平台（GPL）

API 用法：

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

Entrez.email = "pengwenqiang@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)  # 限速

4.3 NCBI PubMed¶

是什么：生物医学文献数据库，收录超 3600 万篇论文的摘要和引用信息。

能查什么： - 论文标题、摘要、作者、期刊 - 关键词、MeSH 术语 - 引用关系

API 用法：

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

Entrez.email = "pengwenqiang@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()

4.4 Ensembl¶

是什么：欧洲生物信息学研究所（EBI）维护的基因组注释数据库，提供基因结构、变异、比较基因组等信息。

能查什么： - 基因坐标、转录本、外显子结构 - 基因功能注释（GO terms） - 物种间同源基因 - 遗传变异（SNP）

API 用法（REST API，不需要额外安装包）：

# ============================================
# 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]}")

4.5 UniProt¶

是什么：全球最权威的蛋白质序列和功能注释数据库，分为 Swiss-Prot（人工审核）和 TrEMBL（自动注释）。

能查什么： - 蛋白质序列、结构域、修饰位点 - 功能注释（GO terms、EC number） - 蛋白互作、亚细胞定位 - 疾病关联

API 用法（2022 年后新 REST API）：

# ============================================
# 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', '')}")

4.6 KEGG（Kyoto Encyclopedia of Genes and Genomes）¶

是什么：京都基因与基因组百科全书，以代谢通路图著称，把基因、酶、代谢物、疾病串联起来。

能查什么： - 代谢通路（Pathway） - 基因 → 酶 → 反应 → 通路的映射关系 - 疾病相关通路 - 微生物功能模块

API 用法（KEGG REST API，免费但有限制）：

# ============================================
# 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 限速

4.7 STRING（蛋白互作网络）¶

是什么：已知和预测的蛋白-蛋白相互作用网络数据库，整合了实验证据、文本挖掘、共表达等多种证据。

能查什么： - 蛋白互作伙伴 - 互作置信度评分 - 功能富集分析 - 互作网络图

API 用法：

# ============================================
# 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}")

5. Entrez 编程接口（Biopython）完整教程¶

5.1 安装和配置¶

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

from Bio import Entrez

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

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

5.2 Entrez 八大函数速查¶

函数	作用	白话解释
`esearch`	搜索，返回 ID 列表	"有哪些符合条件的？"
`efetch`	根据 ID 获取完整记录	"把这几条的详细信息给我"
`esummary`	根据 ID 获取摘要	"简单说说这几条是什么"
`einfo`	查数据库信息	"这个库有哪些字段？"
`elink`	查记录间关联	"这篇文章引用了哪些？"
`epost`	上传 ID 列表到服务器	"先存着，等会儿我要批量查"
`egquery`	跨库搜索计数	"这个词在各个库里各有多少条？"
`espell`	拼写检查	"我是不是拼错了？"

5.3 完整工作流：搜索 → 获取 → 解析¶

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

Entrez.email = "pengwenqiang@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")

5.4 使用 History Server 处理海量结果¶

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

Entrez.email = "pengwenqiang@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")

6. 批量下载 SRA 数据（prefetch + fasterq-dump）¶

6.1 安装 SRA Toolkit¶

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

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

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

6.2 批量下载流程¶

# ============================================
# 批量下载 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 线程

6.3 Python 自动化完整脚本¶

# ============================================
# 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}")

7. 用 Python 批量查询 UniProt REST API¶

# ============================================
# 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]}")

8. 用 Ensembl REST API 获取基因信息¶

# ============================================
# 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')}")

9. 和你 T2D 项目的关联¶

9.1 你的项目能直接用到的查询场景¶

你项目中的步骤	需要批量查询的内容	用哪个库/API
Kraken2 物种注释后	查各菌种的代谢功能注释	KEGG + UniProt
差异丰度分析后	查差异菌的已知 T2D 关联	PubMed Entrez
功能注释（HUMAnN3）后	查 KO/EC 号对应的通路	KEGG REST API
MAGs 基因预测后	批量获取蛋白功能注释	UniProt ID Mapping
写论文讨论部分	查 T2D 菌群相关最新文献	PubMed + Entrez
想扩大样本量	搜索更多 T2D 宏基因组 SRA	Entrez esearch
分析菌群代谢互作	查蛋白互作网络	STRING API

9.2 实战示例：从你的 Kraken2 结果出发¶

# ============================================
# 实战：把 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)  # 限速

10. 常见报错与解决¶

报错 1：HTTP 429 Too Many Requests¶

错误信息：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 端点替代分页查询

报错 2：NCBI Entrez 返回空结果¶

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

报错 3：UniProt 400 Bad Request¶

错误信息：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 先测试单条请求

报错 4：Ensembl 请求超时¶

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

报错 5：fasterq-dump 磁盘空间不足¶

错误信息：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%）

报错 6：KEGG API 返回 403 Forbidden¶

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

11. 速查表：各数据库 API 端点¶

数据库	基础 URL	认证方式	限速	返回格式
NCBI Entrez	`https://eutils.ncbi.nlm.nih.gov/entrez/eutils/`	Email + API Key（可选）	3次/秒（无Key）、10次/秒（有Key）	XML/JSON
Ensembl REST	`https://rest.ensembl.org`	无需认证	15次/秒	JSON/XML/FASTA
UniProt REST	`https://rest.uniprot.org`	无需认证	无明确限制，建议<20次/秒	JSON/TSV/FASTA/XML
KEGG REST	`https://rest.kegg.jp`	无需认证	建议<5次/秒	纯文本/图片
STRING	`https://string-db.org/api`	无需认证	建议<1次/秒	JSON/TSV/图片
InterPro	`https://www.ebi.ac.uk/interpro/api`	无需认证	建议<10次/秒	JSON

常用端点速查¶

=== 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=...       # 互作伙伴列表

12. 延伸资源¶

官方文档（遇到问题先查这些）¶

NCBI Entrez 编程文档：https://www.ncbi.nlm.nih.gov/books/NBK25501/
Biopython Entrez 教程：https://biopython.org/wiki/Entrez
Ensembl REST API 文档：https://rest.ensembl.org/documentation
UniProt REST API 文档：https://www.uniprot.org/help/api
KEGG REST API 说明：https://www.kegg.jp/kegg/rest/keggapi.html
STRING API 文档：https://string-db.org/cgi/help?subpage=api

面试高频问题¶

"你怎么从 SRA 批量下载数据？" → prefetch 下载 .sra，fasterq-dump 转 FASTQ，用 accession list 批量处理
"你会用 API 查询生物数据库吗？" → 用 Biopython Entrez 查 NCBI 全家桶，requests + REST API 查 UniProt/Ensembl/KEGG
"处理大规模数据时怎么避免被封 IP？" → 设 email、注册 API Key、time.sleep 限速、用 History Server 分批获取
"你的 T2D 项目用到了哪些数据库？" → SRA（下载原始数据）、KEGG（通路注释）、UniProt（蛋白功能）、PubMed（文献调研）

最后更新：2026-05-03 | 配合 T2D 宏基因组项目使用

包名	用途	安装
`biopython`	NCBI Entrez + 序列处理	`conda install biopython`
`requests`	通用 HTTP 请求	`conda install requests`
`bioservices`	集成多个生物数据库 API	`pip install bioservices`
`mygene`	批量基因注释查询	`pip install mygene`
`pyensembl`	Ensembl 本地离线查询	`pip install pyensembl`

生信数据库批量查询¶

1. 一句话说明¶

2. 为什么要学¶

你目前的情况¶

3. 核心概念白话版¶

3.1 REST API（表述性状态转移应用编程接口）¶

3.2 Entrez（NCBI 统一检索系统）¶

3.3 批量下载 vs 批量查询¶

3.4 元数据（Metadata）¶

4. 主要数据库详解¶

4.1 NCBI SRA（Sequence Read Archive）¶

4.2 NCBI GEO（Gene Expression Omnibus）¶

4.3 NCBI PubMed¶

4.4 Ensembl¶

4.5 UniProt¶

4.6 KEGG（Kyoto Encyclopedia of Genes and Genomes）¶

4.7 STRING（蛋白互作网络）¶

5. Entrez 编程接口（Biopython）完整教程¶

5.1 安装和配置¶

5.2 Entrez 八大函数速查¶

5.3 完整工作流：搜索 → 获取 → 解析¶

5.4 使用 History Server 处理海量结果¶

6. 批量下载 SRA 数据（prefetch + fasterq-dump）¶

6.1 安装 SRA Toolkit¶

6.2 批量下载流程¶

6.3 Python 自动化完整脚本¶

7. 用 Python 批量查询 UniProt REST API¶

8. 用 Ensembl REST API 获取基因信息¶

9. 和你 T2D 项目的关联¶

9.1 你的项目能直接用到的查询场景¶

9.2 实战示例：从你的 Kraken2 结果出发¶

10. 常见报错与解决¶

报错 1：HTTP 429 Too Many Requests¶

报错 2：NCBI Entrez 返回空结果¶

报错 3：UniProt 400 Bad Request¶

报错 4：Ensembl 请求超时¶

报错 5：fasterq-dump 磁盘空间不足¶

报错 6：KEGG API 返回 403 Forbidden¶

11. 速查表：各数据库 API 端点¶

常用端点速查¶

12. 延伸资源¶

官方文档（遇到问题先查这些）¶

推荐 Python 包¶

面试高频问题¶