Biopython:如何在NCBI中下载所有肽序列(或与特定术语相关的所有记录)?

  • 本文关键字:术语 记录 NCBI Biopython 下载 biopython
  • 更新时间 :
  • 英文 :


我想以fasta格式下载NCBI蛋白质数据库中的所有肽序列(即>和肽名称,后跟肽序列(,我看到这里有一个MESH术语描述了肽是什么,但我无法弄清楚如何整合它。

我写了这个:

import Bio
from Bio import Entrez
Entrez.email = 'test@gmail.com'
handle = Entrez.esearch(db="protein", term="peptide")
record = handle.read()
out_handle = open('myfasta.fasta', 'w')
out_handle.write(record.rstrip('n'))

但它只打印出 995 个 ID,没有要归档的序列,我想知道是否有人可以证明我哪里出错了。

请注意,在NCBI蛋白质数据库中搜索术语肽会返回8187908次命中,因此请确保您确实希望将所有这些命中的肽序列下载到一个大的fasta文件中。

>>> from Bio import Entrez
>>> Entrez.email = 'test@gmail.com'
>>> handle = Entrez.esearch(db="protein", term="peptide")
>>> record = Entrez.read(handle)
>>> record["Count"]
'8187908'

Entrez.esearch 返回的默认记录数为 20。这是为了防止NCBI的服务器过载。

>>> len(record["IdList"])
20

若要获取记录的完整列表,请更改retmax参数:

>>> count = record["Count"]
>>> handle = Entrez.esearch(db="protein", term="peptide", retmax=count)
>>> record = Entrez.read(handle)
>>> len(record['IdList'])
8187908

下载所有记录的方法是使用Entrez.epost

从BioPython教程的第9.4章开始:

EPost 上传 UI 列表,用于后续搜索策略;有关更多信息,请参阅 EPost 帮助页面。它可以通过Bio.Entrez.epost((函数从Biopython获得。

举例说明这何时有用,假设您有一长串要使用 EFetch 下载的 ID(可能是序列,也许是引文 - 任何东西(。当您使用 EFetch 发出请求时,您的 ID 列表、数据库等都会变成发送到服务器的长 URL。如果您的 ID 列表很长,则此 URL 会变长,并且长 URL 可能会中断(例如,某些代理无法很好地处理(。

相反,您可以将其分为两个步骤,首先使用 EPost 上传 ID 列表(这在内部使用"HTML post",而不是"HTML get",从而绕过长 URL 问题(。借助历史记录支持,您可以参考这一长串 ID,并使用 EFetch 下载相关数据。

[...]返回的 XML 包括两个重要的字符串,QueryKey 和 WebEnv,它们共同定义了您的历史记录会话。您可以提取这些值以用于另一个 Entrez 调用,例如 EFetch。

阅读[第9.15章:使用历史记录搜索和下载序列][3]以了解如何使用QueryKeyWebEnv

一个完整的工作示例是:

from Bio import Entrez
import time
from urllib.error import HTTPError
DB = "protein"
QUERY = "peptide"
Entrez.email = 'test@gmail.com'
handle = Entrez.esearch(db=DB, term=QUERY, rettype='fasta')
record = Entrez.read(handle)
count = record['Count']
handle = Entrez.esearch(db=DB, term=QUERY, retmax=count, rettype='fasta')
record = Entrez.read(handle)
id_list = record['IdList']
post_xml = Entrez.epost(DB, id=",".join(id_list))
search_results = Entrez.read(post_xml)
webenv = search_results['WebEnv']
query_key = search_results['QueryKey']

batch_size = 200
with open('peptides.fasta', 'w') as out_handle:
for start in range(0, count, batch_size):
end = min(count, start+batch_size)
print(f"Going to download record {start+1} to {end}")
attempt = 0
success = False
while attempt < 3 and not success:
attempt += 1
try:
fetch_handle = Entrez.efetch(db=DB, rettype='fasta',
retstart=start, retmax=batch_size,
webenv=webenv, query_key=query_key)
success = True
except HTTPError as err:
if 500 <= err.code <= 599:
print(f"Received error from server {err}")
print(f"Attempt {attempt} of 3")
time.sleep(15)
else:
raise
data = fetch_handle.read()
fetch_handle.close()
out_handle.write(data)

peptides.fasta的前几行如下所示:

>QGT67293.1 RepA leader peptide Tap (plasmid) [Klebsiella pneumoniae]
MLRKLQAQFLCHSLLLCNISAGSGD
>QGT67288.1 RepA leader peptide Tap (plasmid) [Klebsiella pneumoniae]
MLRKLQAQFLCHSLLLCNISAGSGD
>QGT67085.1 thr operon leader peptide [Klebsiella pneumoniae]
MNRIGMITTIITTTITTGNGAG
>QGT67083.1 leu operon leader peptide [Klebsiella pneumoniae]
MIRTARITSLLLLNACHLRGRLLGDVQR
>QGT67062.1 peptide antibiotic transporter SbmA [Klebsiella pneumoniae]
MFKSFFPKPGPFFISAFIWSMLAVIFWQAGGGDWLLRVTGASQNVAISAARFWSLNYLVFYAYYLFCVGV
FALFWFVYCPHRWQYWSILGTSLIIFVTWFLVEVGVAINAWYAPFYDLIQSALATPHKVSINQFYQEIGV
FLGIAIIAVIIGVMNNFFVSHYVFRWRTAMNEHYMAHWQHLRHIEGAAQRVQEDTMRFASTLEDMGVSFI
NAVMTLIAFLPVLVTLSEHVPDLPIVGHLPYGLVIAAIVWSLMGTGLLAVVGIKLPGLEFKNQRVEAAYR
KELVYGEDDETRATPPTVRELFRAVRRNYFRLYFHYMYFNIARILYLQVDNVFGLFLLFPSIVAGTITLG
LMTQITNVFGQVRGSFQYLISSWTTLVELMSIYKRLRSFERELDGKPLQEAIPTLR

最新更新