尧图网站设计 尧图网站设计YAOTU DESIGN
ARTICLE DETAIL

资讯详情

深耕网站设计与一线实操的经验洞察。

使用Biopython与NCBI Entrez批量获取RefSeq标准转录本(NM_)注释信息

使用Biopython与NCBI Entrez批量获取RefSeq标准转录本(NM_)注释信息 1. 项目概述从“NM_”说起我们到底在找什么如果你在生物信息学或者分子生物学领域摸爬滚打过一阵子肯定对“NM_”这个前缀不陌生。它就像基因世界里一个特定产品的“标准型号编码”。简单来说NM_开头的编号是NCBI美国国家生物技术信息中心RefSeq数据库中为一条经过人工或计算审校、具有明确生物学意义的“标准”mRNA转录本所分配的标识符。它代表的是一个基因可能产生的、被广泛认可和引用的“模范”转录本序列。所以当有人问“如何下载关于人类或者其它物种的全部转录本名称NM_”他真正的需求往往不止是拿到一串以“NM_”开头的ID列表。更深层的需求可能是我需要一份高质量的、权威的基因注释参考列表用于设计引物、进行表达量分析比如RNA-seq、构建基因集做富集分析或者仅仅是整理一个可靠的基因目录。这个需求听起来简单但实操起来有几个关键点容易被忽略。首先“全部”是一个动态概念基因组注释在不断更新新的转录本会被发现和验证。其次“人类或者其它物种”意味着方法需要具备通用性能适应不同物种的数据库。最后直接拿到“名称”往往不够我们通常还需要与之关联的基因符号、基因描述、染色体位置等信息形成一个可用的数据表格。因此这个项目的核心是掌握从权威生物数据库主要是NCBI中以编程化、可重复的方式批量获取并处理标准转录本注释信息的能力。下面我就以一个多年“数据矿工”的经验带你走通这条路并分享那些手册上不会写的细节和坑。2. 核心思路与工具选型为什么是Entrez和Biopython面对从海量数据库中获取特定数据的需求我们首先得摒弃手动在网站上点击下载的念头。对于人类约2万多个基因每个基因可能有多个转录本或其他物种如小鼠、斑马鱼数据量可能从几万到几十万条记录自动化是唯一可行的路径。2.1 为什么选择NCBI Entrez系统在生物数据领域NCBI是不折不扣的“中央数据库”。而Entrez是其提供的一套用于检索和获取数据的编程接口API它就像数据库的“总开关”。通过E-utilitiesEntrez工具集我们可以用URL链接或专门的库来查询和下载数据。选择它的理由很充分权威性与完整性RefSeqNM_所在库是NCBI维护的黄金标准参考序列数据库数据质量高注释可靠。统一入口Entrez可以跨多个NCBI数据库如Gene, Nucleotide, Protein进行联合查询方便我们获取关联信息。稳定与免费作为政府资助的公共资源其API稳定且对学术用途完全免费。标准化输出支持多种格式如XML ASN.1, FASTA便于程序解析。2.2 为什么选择Biopython进行本地操作虽然可以直接用requests库构造URL调用E-utilities但我强烈推荐使用Biopython的Bio.Entrez模块。它是一个专门为生物信息学设计的Python工具包封装了Entrez的复杂参数让操作变得异常简单。它的优势在于简化流程自动处理API密钥、请求频率限制NCBI要求每秒不超过3次请求、错误重试等琐事。直接解析返回的数据可以直接被解析成Python对象比如BeautifulSoup解析XML或特定的记录对象省去自己写解析器的麻烦。生态丰富Biopython还包含大量其他生物信息学工具后续数据处理如序列操作、格式转换可以无缝衔接。注意使用NCBI的E-utilities时务必遵守其使用规范。强烈建议在代码中设置你的邮箱地址Entrez.email “your.emailexample.com“这样当你的请求出现问题时NCBI管理员可以联系你。如果进行大规模批量请求最好申请一个API密钥Entrez.api_key “your_api_key“这样可以将请求频率限制从每秒3次提升到每秒10次。3. 实战演练分步获取人类全部NM_转录本理论说完我们直接上代码。我们的目标是获取智人Homo sapiens所有RefSeq数据库中状态为“reviewed”已审阅的、NM_开头的mRNA转录本ID并附带其对应的基因符号和完整描述。3.1 环境准备与Biopython安装首先确保你的Python环境建议3.7以上已经就绪。安装Biopython非常简单pip install biopython如果下载速度慢可以使用国内镜像源例如pip install biopython -i https://pypi.tuna.tsinghua.edu.cn/simple3.2 第一步在Nucleotide库中定位目标数据集我们需要的NM_记录存放在NCBI的“Nucleotide”数据库中。第一步是进行搜索以获取符合条件的所有记录的ID列表即UID列表。这里的关键是构建正确的搜索词Term。from Bio import Entrez import time # 1. 设置你的邮箱这是礼貌也是规范 Entrez.email “your_real_emailinstitution.com“ # 请务必替换成你的真实邮箱 # 2. 定义搜索参数 db “nucleotide“ # 数据库核苷酸序列库 # 核心搜索词物种为“Homo sapiens”[Organism]且序列版本前缀为“NM_”[Filter] # 加上“AND srcdb_refseq[prop]”可以更精确地限定在RefSeq数据库内 # 加上“AND alive[prop]”确保获取当前有效的记录非撤回版本 term ‘(“Homo sapiens“[Organism]) AND (NM_[Filter]) AND srcdb_refseq[prop] AND alive[prop]‘ # 3. 使用esearch进行搜索获取查询结果句柄和ID列表 print(“正在搜索并获取ID列表...“) search_handle Entrez.esearch(dbdb, termterm, retmax100000, usehistory“y“) search_results Entrez.read(search_handle) search_handle.close() # 提取关键信息 count int(search_results[“Count“]) webenv search_results[“WebEnv“] query_key search_results[“QueryKey“] id_list search_results[“IdList“] print(f“找到总计 {count} 条记录。“) print(f“前10个ID示例{id_list[:10]}“)关键参数解析retmax单次返回的最大记录数。这里设为10万对于人类转录本通常足够。如果超过这个数需要分批次获取。usehistory“y“这是核心技巧对于大型结果集NCBI推荐使用“历史会话”模式。它不会一次性返回所有ID可能超时或数据包过大而是返回一个WebEnv网络环境句柄和QueryKey查询键。后续我们可以用这两个参数以分批batch的方式安全地获取数据这是处理大数据集的标准做法。3.3 第二步分批获取完整的记录详情有了WebEnv和QueryKey我们就可以分批获取每条记录的详细XML信息。这里我们使用efetch函数并指定返回格式为“xml”。# 4. 配置分批获取参数 batch_size 500 # 每批获取500条记录。NCBI建议单次不超过500条以防服务器过载。 all_records [] print(“开始分批获取记录详情...“) for start in range(0, count, batch_size): end min(count, start batch_size) print(f“正在获取第 {start1} 到 {end} 条记录...“) # 使用历史会话参数进行获取 fetch_handle Entrez.efetch( dbdb, rettype“gb“, # 返回GenBank格式信息最全。也可用“fasta”只取序列。 retmode“xml“, # 以XML格式返回便于用Entrez.read解析 retstartstart, retmaxbatch_size, webenvwebenv, query_keyquery_key ) # 解析XML数据 batch_records Entrez.read(fetch_handle) all_records.extend(batch_records) fetch_handle.close() # 遵守NCBI的请求频率限制每批请求后暂停 time.sleep(0.34) # 略高于1/3秒确保平均每秒请求不超过3次 print(f“所有批次获取完成共解析 {len(all_records)} 条记录。“)实操心得rettype“gb“和retmode“xml“的组合能获取到包含基因注释、来源、参考文献等最丰富信息的结构化数据。如果你只需要序列可以改用rettype“fasta“,retmode“text“。time.sleep(0.34)这个暂停至关重要。不加的话频繁请求会被NCBI服务器暂时封禁。如果你申请了API密钥可以适当缩短间隔如0.11秒。3.4 第三步从XML记录中提取目标信息获取到的all_records是一个复杂的嵌套数据结构Python列表和字典。我们需要从中精准地提取出NM_编号、基因符号、描述等信息。这需要你对GenBank XML的结构有一定了解。# 5. 解析并提取信息 transcript_info_list [] for record in all_records: # 每条record是一个字典对应一条序列记录 accession record.get(“GBSeq_primary-accession“, “N/A“) # 主登录号如NM_001370658.2 version record.get(“GBSeq_accession-version“, “N/A“) # 带版本号的登录号更精确 # 提取基因符号和完整名称描述 gene_symbol “N/A“ full_name “N/A“ # 注释信息在 GBSeq_feature-table 中 features record.get(“GBSeq_feature-table“, []) for feature in features: if feature.get(“GBFeature_key“) “gene“: # 在gene特征的限定符(qualifiers)里找gene符号 quals feature.get(“GBFeature_quals“, []) for qual in quals: if qual.get(“GBQualifier_name“) “gene“: gene_symbol qual.get(“GBQualifier_value“, “N/A“) break elif feature.get(“GBFeature_key“) “CDS“: # CDS特征里的“product”限定符通常是蛋白产物名可作为描述的一部分 quals feature.get(“GBFeature_quals“, []) for qual in quals: if qual.get(“GBQualifier_name“) “product“: full_name qual.get(“GBQualifier_value“, “N/A“) # 有时更完整的描述在记录的GBSeq_definition字段 if full_name “N/A“: full_name record.get(“GBSeq_definition“, “N/A“) break # 如果gene特征里没找到可以尝试从记录的GBSeq_definition中解析 if gene_symbol “N/A“: def_line record.get(“GBSeq_definition“, ““) # 一个简单的启发式解析定义行通常格式为“基因符号 蛋白/产品描述” parts def_line.split() if parts: # 假设第一个单词是基因符号这不总是成立但常见 potential_symbol parts[0] # 过滤掉一些明显不是基因符号的常见词 if not potential_symbol.startswith((“PREDICTED:“, “SIMILAR“, “Homo“)): gene_symbol potential_symbol transcript_info_list.append({ “Accession“: accession, “Accession_Version“: version, “Gene_Symbol“: gene_symbol, “Description“: full_name }) print(f“信息提取完成共处理 {len(transcript_info_list)} 条转录本记录。“)注意事项 XML解析是这一步最易出错的地方。NCBI的XML结构并非一成不变不同记录间可能存在字段缺失或结构差异。上述代码提供了基本的提取逻辑但为了鲁棒性在实际生产环境中你需要添加更多的错误处理try...except和备用提取路径。例如基因符号也可能存在于“source“特征的“gene“限定符中。3.5 第四步数据输出与保存最后我们将提取的信息保存为最通用的CSV格式方便用Excel、R或Pandas进行后续分析。import pandas as pd # 6. 转换为DataFrame并保存 df pd.DataFrame(transcript_info_list) # 去重理论上不应该有重复但操作一下更安全 df.drop_duplicates(subset[“Accession_Version“], inplaceTrue) output_file “human_refseq_nm_transcripts.csv“ df.to_csv(output_file, indexFalse, encoding‘utf-8-sig‘) # utf-8-sig确保Excel打开中文不乱码 print(f“数据已成功保存至 {output_file}“) print(f“数据预览\n“, df.head())至此一个包含人类全部当前数据库中的NM_转录本ID、基因符号和描述的表格就生成了。你可以用文本编辑器或Excel打开这个CSV文件进行查看。4. 方案扩展适配其它物种与获取更多信息上面的流程是针对人类的。如何适配其他物种呢关键在于修改搜索词term中的物种限定部分。4.1 修改物种例如要获取小鼠Mus musculus的数据term_mouse ‘(“Mus musculus“[Organism]) AND (NM_[Filter]) AND srcdb_refseq[prop] AND alive[prop]‘要获取大鼠Rattus norvegicus的数据term_rat ‘(“Rattus norvegicus“[Organism]) AND (NM_[Filter]) AND srcdb_refseq[prop] AND alive[prop]‘你可以在NCBI的网站上先搜索目标物种的学名确保拼写正确。通常使用双引号包裹完整的物种学名是最精确的。4.2 获取染色体位置信息有时我们还需要转录本的染色体位置chr, start, end, strand。这些信息也包含在XML记录中但位于“GBSeq_feature-table“里key为“source“的特征中。def extract_genomic_location(feature_table): “““从特征表中提取基因组位置信息。“““ chrom, start, end, strand “N/A“, “N/A“, “N/A“, “N/A“ for feature in feature_table: if feature.get(“GBFeature_key“) “source“: location feature.get(“GBFeature_location“, ““) # 位置字符串格式可能为“1..1000”或“complement(1..1000)” if location: # 简单解析示例实际可能需要更复杂的正则表达式 if “complement“ in location: strand “-“ location location.replace(“complement(“, ““).replace(“)“, ““) else: strand ““ if “..“ in location: start, end location.split(“..“)[:2] # 尝试从source的限定符中找染色体号 quals feature.get(“GBFeature_quals“, []) for qual in quals: if qual.get(“GBQualifier_name“) “chromosome“: chrom qual.get(“GBQualifier_value“, “N/A“) elif qual.get(“GBQualifier_name“) “map“: # 旧记录可能用map字段 map_info qual.get(“GBQualifier_value“, ““) if map_info.startswith(“chr“): chrom map_info break return chrom, start, end, strand # 在之前的循环中可以调用此函数 for record in all_records: # ... 之前的提取代码 ... features record.get(“GBSeq_feature-table“, []) chrom, start, end, strand extract_genomic_location(features) # 将chrom, start, end, strand加入到transcript_info_list的字典中4.3 使用NCBI Gene数据库作为补充或替代途径除了从Nucleotide库抓取另一个更直接的获取基因-转录本对应关系的方法是查询NCBI Gene数据库。Gene库整合了基因的综合信息通常能更清晰地列出该基因的所有RefSeq转录本包括NM_和NR_等。# 搜索人类所有基因 term_gene ‘(“Homo sapiens“[Organism]) AND alive[property]‘ search_handle_gene Entrez.esearch(db“gene“, termterm_gene, retmax5000, usehistory“y“) gene_results Entrez.read(search_handle_gene) search_handle_gene.close() gene_count int(gene_results[“Count“]) print(f“在Gene库中找到 {gene_count} 个基因记录。“) # 然后可以获取每个基因的摘要(summary)其中包含“Genomic context”和“RefSeq”部分 # 或者使用elink获取基因关联的核苷酸记录 # 这种方法更适合获取基因层面的汇总信息对于精确抓取所有NM_ ID可能不如直接查询Nucleotide库全面。5. 常见问题、排查技巧与优化建议在实际操作中你几乎一定会遇到下面这些问题。这里是我踩过坑后总结的解决方案。5.1 网络连接与超时问题问题HTTPError或URLError连接被重置或超时。排查检查网络确保可以正常访问NCBI官网https://www.ncbi.nlm.nih.gov/。某些网络环境可能需要配置代理。添加超时和重试在Entrez.efetch或esearch调用中增加timeout参数如timeout30。并可以用try...except包裹请求实现失败重试。import urllib.error max_retries 3 for attempt in range(max_retries): try: fetch_handle Entrez.efetch(..., timeout30) # ... 处理逻辑 ... break # 成功则跳出循环 except (urllib.error.URLError, urllib.error.HTTPError) as e: print(f“请求失败第{attempt1}次重试。错误{e}“) time.sleep(5) # 重试前等待更久 else: print(“多次重试后仍失败请检查网络或稍后再试。“)5.2 数据不完整或缺失字段问题提取到的基因符号或描述大量为“N/A”。排查检查XML结构打印一条完整记录的XMLprint(record)仔细查看目标信息到底藏在哪个路径下。不同生物类型或不同提交日期的记录结构可能有细微差别。使用备用字段如果“gene“特征里没有尝试从“CDS“特征的“gene“限定符或“source“特征的“gene“限定符中提取。描述信息也可以优先使用“GBSeq_definition“。考虑使用其他工具如果对数据质量要求极高可以考虑使用NCBI Datasets命令行工具datasets或BioMart通过biomartR包或Ensembl网站。这些工具提供了更稳定、格式更统一的基因注释文件下载。5.3 处理大规模数据数十万条记录问题获取全部记录耗时极长甚至中途失败。优化建议善用历史会话我们之前已经用了usehistory“y“这是必须的。调整批次大小batch_size可以调小如200或250减少单次请求的数据量提高成功率。使用API密钥前往NCBI账户设置申请API密钥并在代码中设置Entrez.api_key将速率限制从3次/秒提升到10次/秒能显著缩短总时间。保存中间状态对于超大规模任务可以将每批获取并解析后的数据立即追加保存到文件如JSON Lines格式.jsonl。这样即使程序中断也可以从断点恢复避免前功尽弃。考虑离线数据对于人类、小鼠等模式生物NCBI提供预编译的注释文件如refseq_*.gff.gz或*.gtf.gz。直接从FTPftp://ftp.ncbi.nlm.nih.gov/下载这些文件然后用gffread或编程方式解析通常是最快最稳定的方法。这尤其适合需要频繁获取相同物种数据的情况。5.4 结果验证与去重问题如何确保我拿到的是“全部”且“最新”的NM_记录验证步骤数量级核对人类蛋白质编码基因约2万个但许多基因有多个转录本所以NM_总数通常在8-10万左右。如果你的结果远小于这个数可能是搜索词太严格比如漏掉了某些变体或解析过程丢失了大量记录。交叉检查随机抽取几十个你熟悉的基因检查其主要的NM_转录本是否都在你的列表中。版本号注意Accession_Version如NM_001370658.2中的版本号。同一个NM_登录号版本号更新意味着序列有修订。我们的搜索词alive[prop]确保了获取的是当前有效版本。去重务必基于Accession_Version登录号版本进行去重因为这是记录的唯一标识。仅用Accession不带版本去重可能会遗漏不同版本的信息。最后记住生物数据库是动态的。今天下载的“全量”数据三个月后可能就不全了。对于长期项目建议将数据获取脚本化并记录下载日期和数据库版本或者在关键分析前重新运行脚本获取最新数据。这套基于Entrez和Biopython的流程为你提供了一个强大、灵活且可编程的解决方案足以应对绝大多数从NCBI获取标准转录本列表的需求。
返回列表