Python高效下载NCBI生物数据:从Entrez到SRA实用指南213


在生物信息学研究中,美国国家生物技术信息中心(NCBI)是全球最重要、最全面的生物数据宝库之一。它包含了基因组序列、蛋白质序列、基因表达数据、文献信息以及高通量测序原始数据等海量信息。对于生物学家和生物信息学工程师而言,高效、自动化地获取NCBI数据是日常工作不可或缺的一部分。Python,凭借其强大的生态系统和简洁的语法,成为了完成这项任务的理想工具。

本文将深入探讨如何使用Python来下载NCBI的各类生物数据,从基础的Entrez E-utilities查询到复杂的高通量测序数据(SRA)获取,并提供实用的代码示例和最佳实践,帮助读者构建健壮、高效的数据下载流程。

理解NCBI数据及其访问方式

NCBI提供了多种数据类型,每种类型都有其特定的存储方式和访问接口:
序列数据(Sequences): 包括GenBank、RefSeq、dbSNP、蛋白质序列等。主要通过Entrez E-utilities访问。
文献数据(Literature): 如PubMed,收录了生物医学文献摘要和全文。同样通过Entrez E-utilities访问。
基因表达数据(Gene Expression): 如GEO(Gene Expression Omnibus),存储芯片和RNA-seq实验数据。可以通过E-utilities查询,但原始数据通常通过FTP下载。
高通量测序原始数据(SRA): Sequence Read Archive,存储原始测序reads。可以通过Entrez E-utilities查询元数据,但数据文件本身需要专门的SRA Toolkit进行下载和处理。
基因组数据: 完整的基因组序列、注释信息。通常通过FTP或NCBI Genomes下载。

Python主要通过以下几种方式与NCBI交互:
BioPython库的Entrez模块: 这是最常用和推荐的方式,它封装了NCBI的E-utilities Web服务,允许通过HTTP请求进行数据查询和获取。
`requests`库: 对于需要直接进行HTTP或FTP下载的场景,特别是下载大文件或处理非E-utilities接口的数据时,`requests`库提供了更灵活的控制。
`subprocess`模块: 用于在Python脚本中调用外部命令行工具,例如SRA Toolkit。

使用BioPython Entrez模块进行数据检索与下载

BioPython是生物信息学领域最受欢迎的Python库之一,其``模块是访问NCBI E-utilities的核心。

1. 准备工作


首先,确保你已经安装了BioPython库:pip install biopython

在使用Entrez模块之前,你必须设置``,这是NCBI要求用于识别用户身份的。虽然不是强制的,但强烈建议注册一个NCBI账户并获取API Key,这可以显著提高查询限制,减少因请求频繁而被限流的风险。from Bio import Entrez
import time
# 替换为你的邮箱地址
= "@"
# 可选:如果注册了NCBI API Key,请在此处设置,可以提高查询限制
# Entrez.api_key = "YOUR_NCBI_API_KEY"

2. E-utilities核心概念:esearch 和 efetch


Entrez E-utilities的核心是`esearch`和`efetch`:
`esearch`:用于在指定NCBI数据库中执行搜索查询,返回匹配的记录ID列表。
`efetch`:根据`esearch`返回的ID列表,获取实际的数据内容(如序列、文献摘要等),并可指定输出格式。

3. 实例:下载GenBank序列数据


假设我们要下载某个基因的GenBank序列,我们可以通过基因名或Accession ID进行搜索。from Bio import Entrez
from Bio import SeqIO
import os
= "@"
# Entrez.api_key = "YOUR_NCBI_API_KEY"
def download_genbank_sequence(query, db="nuccore", retmode="xml", rettype="gb", output_dir="ncbi_data"):
"""
通过基因名或Accession ID下载GenBank序列。
Args:
query (str): 搜索关键词或Accession ID。
db (str): NCBI数据库名,例如"nuccore"(核苷酸)、"protein"(蛋白质)。
retmode (str): 返回模式,例如"xml"、"text"。
rettype (str): 返回类型,例如"gb"(GenBank格式)、"fasta"。
output_dir (str): 保存文件的目录。
"""
(output_dir, exist_ok=True)

print(f"Searching for '{query}' in {db} database...")
handle = (db=db, term=query, retmax="10") # retmax限制返回条目数量
record = (handle)
()

id_list = record["IdList"]
if not id_list:
print(f"No records found for '{query}'.")
return
print(f"Found {len(id_list)} IDs: {', '.join(id_list)}")
for i, seq_id in enumerate(id_list):
print(f"Fetching sequence for ID: {seq_id}...")
try:
fetch_handle = (db=db, id=seq_id, rettype=rettype, retmode=retmode)
if retmode == "xml":
# 对于XML格式,需要用解析
with open((output_dir, f"{seq_id}.gb"), "w") as out_handle:
(())
# 若要解析为SeqRecord对象:
# records = list((fetch_handle, "gb"))
# if records: print(f"First record length: {len(records[0].seq)}")
else: # 例如rettype="fasta", retmode="text"
with open((output_dir, f"{seq_id}.fasta"), "w") as out_handle:
(())
print(f"Successfully downloaded {seq_id} to {output_dir}.")
except Exception as e:
print(f"Error downloading {seq_id}: {e}")
finally:
()

(0.5) # 遵守NCBI的请求频率限制
# 示例:下载大肠杆菌O157:H7的16S rRNA基因
download_genbank_sequence("Escherichia coli O157:H7 16S rRNA gene", db="nuccore", rettype="fasta", retmode="text")
# 示例:下载特定的Accession ID
download_genbank_sequence("NC_000913", db="nuccore", rettype="gb", retmode="text") # E. coli K-12 MG1655 genome

4. 实例:下载PubMed文献摘要


除了序列,我们也可以下载PubMed文献的摘要信息。from Bio import Entrez
import json # 用于保存解析后的数据
= "@"
# Entrez.api_key = "YOUR_NCBI_API_KEY"
def download_pubmed_abstracts(query, output_file="", retmax="50"):
"""
下载PubMed文献摘要。
Args:
query (str): 搜索关键词,例如"CRISPR-Cas9 AND cancer"。
output_file (str): 保存摘要的JSON文件名。
retmax (str): 返回的最大条目数。
"""
print(f"Searching for '{query}' in PubMed...")
handle = (db="pubmed", term=query, retmax=retmax)
record = (handle)
()
id_list = record["IdList"]
if not id_list:
print(f"No PubMed articles found for '{query}'.")
return
print(f"Found {len(id_list)} PubMed IDs. Fetching abstracts...")
abstracts = []
# 分批下载以避免一次性请求过多数据
batch_size = 10
for i in range(0, len(id_list), batch_size):
batch_ids = id_list[i : i + batch_size]
try:
fetch_handle = (db="pubmed", id=",".join(batch_ids), retmode="xml")
pubmed_records = (fetch_handle)
()
for pubmed_article in pubmed_records["PubmedArticle"]:
article = pubmed_article["MedlineCitation"]["Article"]
title = ("ArticleTitle", "No Title")
# 尝试从Abstract标签中提取文本
abstract_text = ""
if "Abstract" in article and "AbstractText" in article["Abstract"]:
# AbstractText可能是一个列表或单个字符串
if isinstance(article["Abstract"]["AbstractText"], list):
abstract_text = " ".join(article["Abstract"]["AbstractText"])
else:
abstract_text = article["Abstract"]["AbstractText"]
({
"PMID": pubmed_article["MedlineCitation"]["PMID"],
"Title": str(title), # 确保是字符串
"Abstract": str(abstract_text) # 确保是字符串
})
print(f"Fetched {len(batch_ids)} abstracts in this batch.")
except Exception as e:
print(f"Error fetching batch with IDs {batch_ids}: {e}")
finally:
(1) # 遵守NCBI请求频率限制
with open(output_file, "w", encoding="utf-8") as f:
(abstracts, f, indent=4, ensure_ascii=False)
print(f"Successfully downloaded {len(abstracts)} PubMed abstracts to {output_file}.")
# 示例:下载关于"CRISPR-Cas9"的PubMed摘要
download_pubmed_abstracts("CRISPR-Cas9", output_file="", retmax="100")

处理SRA数据:使用BioPython与SRA Toolkit

SRA (Sequence Read Archive) 存储了高通量测序的原始数据,如Illumina、PacBio等平台的Reads。SRA数据通常非常庞大,其下载和处理流程与GenBank序列不同。

1. SRA数据特点与下载工具


SRA数据文件格式是`.sra`,这是一种压缩且带有元数据的特殊格式。NCBI提供了一个专门的命令行工具集——SRA Toolkit——来处理这些文件,包括下载(`prefetch`)和转换为FASTQ格式(`fastq-dump`)。

在使用Python下载SRA数据前,你需要先安装SRA Toolkit并确保其可执行文件(如`prefetch`, `fastq-dump`)在系统的PATH中。安装方法请参考NCBI官方文档。

2. 通过BioPython获取SRA Accession信息


我们可以使用BioPython的`esearch`模块在`sra`数据库中搜索SRA Accession IDs。from Bio import Entrez
import subprocess
import os
import time
= "@"
# Entrez.api_key = "YOUR_NCBI_API_KEY"
def get_sra_accessions(query, retmax="10"):
"""
在SRA数据库中搜索SRA Accession IDs。
Args:
query (str): 搜索关键词,例如"SRR"、"GSE"、"ERP"等。
retmax (str): 返回的最大条目数。
Returns:
list: 匹配的SRA Accession ID列表。
"""
print(f"Searching for '{query}' in SRA database...")
handle = (db="sra", term=query, retmax=retmax)
record = (handle)
()

id_list = record["IdList"]
if not id_list:
print(f"No SRA records found for '{query}'.")
return []
# SRA的esearch返回的ID是Entrez内部ID,需要efetch来获取SRA accession号
print(f"Found {len(id_list)} internal SRA IDs. Fetching SRA accessions...")
try:
fetch_handle = (db="sra", id=",".join(id_list), retmode="xml")
sra_records_xml = ()
()

# 解析XML,提取SRA Accession
# 注意:可能无法直接解析所有SRA XML格式,可能需要手动解析
# 这里用一个简化的方式,实际情况可能需要更复杂的XML解析
sra_accessions = []
from import ElementTree as ET
root = (sra_records_xml)
for experiment_package in (".//Run"): # 查找所有的标签
accession = ("accession")
if accession:
(accession)

print(f"Extracted {len(sra_accessions)} SRA accessions.")
return sra_accessions
except Exception as e:
print(f"Error fetching SRA accessions: {e}")
return []
# 示例:搜索与"SARS-CoV-2"相关的SRA数据
# sra_ids = get_sra_accessions("SARS-CoV-2", retmax="5")
# print(f"SRA Accession IDs: {sra_ids}")

3. 调用SRA Toolkit下载和转换数据


一旦获取了SRA Accession ID(如SRRXXXXXXX),我们就可以使用`subprocess`模块来调用SRA Toolkit的`prefetch`和`fastq-dump`命令。# 接续上面的代码
def download_and_convert_sra(sra_accession_ids, output_dir="sra_data"):
"""
使用SRA Toolkit下载SRA数据并转换为FASTQ格式。
Args:
sra_accession_ids (list): SRA Accession ID列表。
output_dir (str): 保存SRA文件和FASTQ文件的目录。
"""
(output_dir, exist_ok=True)

for sra_id in sra_accession_ids:
sra_file_path = (output_dir, f"{sra_id}.sra")

# 1. 使用prefetch下载SRA文件
print(f"Downloading SRA {sra_id} with prefetch...")
try:
# -O 指定输出目录
prefetch_cmd = ["prefetch", sra_id, "-O", output_dir]
result = (prefetch_cmd, capture_output=True, text=True, check=True)
print(f"Prefetch stdout:{}")
if :
print(f"Prefetch stderr:{}")
print(f"Successfully prefetched {sra_id}.")
except as e:
print(f"Error during prefetch for {sra_id}: {e}")
print(f"Stdout: {}")
print(f"Stderr: {}")
continue

# 2. 使用fastq-dump将SRA文件转换为FASTQ
# --split-files 用于双端测序数据(生成和)
# --gzip 压缩FASTQ文件
print(f"Converting {sra_id}.sra to FASTQ with fastq-dump...")
try:
# -O 指定输出目录,fastq-dump默认输出到当前目录或-O指定的目录
fastqdump_cmd = ["fastq-dump", "--split-files", "--gzip", sra_file_path, "-O", output_dir]
result = (fastqdump_cmd, capture_output=True, text=True, check=True)
print(f"Fastq-dump stdout:{}")
if :
print(f"Fastq-dump stderr:{}")
print(f"Successfully converted {sra_id} to FASTQ.")
except as e:
print(f"Error during fastq-dump for {sra_id}: {e}")
print(f"Stdout: {}")
print(f"Stderr: {}")
continue

# 可选:转换完成后删除.sra文件
# if (sra_file_path):
# (sra_file_path)
# print(f"Removed {sra_file_path}.")
# 示例:下载并转换一些SARS-CoV-2相关的SRA数据
# 注意:SRA文件通常很大,下载耗时较长
sra_ids_to_download = ["SRR11407008", "SRR11407009"] # 替换为你想要下载的SRA Accession ID
download_and_convert_sra(sra_ids_to_download, output_dir="sars_cov_2_sra_data")

请注意,SRA数据文件通常非常大(几十GB甚至几百GB),下载需要较长时间和大量磁盘空间。确保你的网络连接稳定,并有足够的存储空间。

批量下载与高级技巧

1. 速率限制与API Key


NCBI对E-utilities的请求频率有限制。不使用API Key时,每秒不能超过3个请求。使用API Key后,限制会放宽到每秒10个请求。务必在代码中加入`()`来控制请求频率,避免被NCBI暂时封锁IP。

2. 分批处理(Batch Processing)


当需要下载大量ID时,最好将ID列表分成小批次(例如每次100到500个ID)进行处理,而不是一次性提交所有ID。这有助于减少服务器负载,并简化错误处理。

3. 错误处理与重试机制


网络请求和文件操作都可能失败。使用`try-except`块捕获异常,并考虑实现简单的重试机制,例如在失败后等待一段时间再尝试几次。import time
def safe_efetch(db, id_list, rettype, retmode, max_retries=3, delay=5):
"""带重试机制的efetch函数。"""
for attempt in range(max_retries):
try:
handle = (db=db, id=",".join(id_list), rettype=rettype, retmode=retmode)
return handle
except Exception as e:
print(f"Attempt {attempt + 1} failed: {e}. Retrying in {delay} seconds...")
(delay)
raise Exception(f"Failed to fetch after {max_retries} attempts for IDs: {id_list}")

4. 大文件下载


对于通过URL直接下载的大文件(如GEO的原始数据),可以使用`requests`库的流式下载功能,避免一次性将整个文件载入内存:import requests
def download_large_file(url, local_filename, chunk_size=8192):
"""
流式下载大文件。
"""
print(f"Downloading {url} to {local_filename}...")
try:
with (url, stream=True) as r:
r.raise_for_status() # 检查HTTP请求是否成功
with open(local_filename, 'wb') as f:
for chunk in r.iter_content(chunk_size=chunk_size):
(chunk)
print(f"Successfully downloaded {local_filename}.")
except as e:
print(f"Error downloading {url}: {e}")
# 示例:下载一个虚构的大文件URL
# download_large_file("/geo/series/GSEXXXXX/", "")

5. 进度条显示


对于长时间运行的下载任务,使用`tqdm`库可以提供友好的进度条,提高用户体验:pip install tqdm
from tqdm import tqdm
# ... (在 download_large_file 函数中集成 tqdm)
def download_large_file_with_progress(url, local_filename, chunk_size=8192):
print(f"Downloading {url} to {local_filename}...")
try:
with (url, stream=True) as r:
r.raise_for_status()
total_size = int(('content-length', 0))
with open(local_filename, 'wb') as f, tqdm(
total=total_size, unit='B', unit_scale=True, unit_divisor=1024,
desc=('/')[-1]
) as bar:
for chunk in r.iter_content(chunk_size=chunk_size):
(chunk)
(len(chunk))
print(f"Successfully downloaded {local_filename}.")
except as e:
print(f"Error downloading {url}: {e}")

总结与展望

Python在NCBI生物数据下载方面展现出卓越的灵活性和强大功能。通过BioPython的Entrez模块,我们可以方便地查询和下载各种数据库(如GenBank、PubMed)的文本或XML格式数据。对于高通量测序原始数据(SRA),结合`subprocess`模块调用SRA Toolkit,可以实现从元数据查询到原始文件下载和格式转换的完整自动化流程。

掌握这些技巧,生物信息学研究人员将能够更高效地获取所需数据,为后续的数据分析和科学研究打下坚实基础。在实际应用中,务必关注NCBI的数据使用政策和请求频率限制,设计健壮的错误处理和重试机制,以确保数据下载过程的稳定性和可靠性。

2025-10-11


上一篇:Python多线程数据并发处理:深入理解与最佳实践

下一篇:Python读取Excel数据:高效、灵活的数据处理实战指南