相信大家在上一文中下载fasta的时候还没有感觉到下载是多么复杂,但是对于分析比对多个序列文件时,这个工作量说多了都是泪。比如,老板让你比对自己测定序列与 NCBI 库中序列,并构建相应的进化树,而这个序列需要大于100条。我想你的心情不会和下载一条序列时那么平静,那么,接下来通过BioPython提供的接口来实现快速的自动化序列下载。
一、Entrez 库
1.1 Entrez 介绍
Entrez 在线资源检索器是一组服务器端程序,为国家生物技术信息中心(NCBI)的Entrez查询和数据库系统提供稳定的接口。使用固定的URL语法,将一组标准输入参数转换为各种NCBI软件组件搜索和检索所请求数据所需的值。目前包括38个数据库,涵盖各种生物医学数据,包括核苷酸和蛋白质序列,基因记录,三维分子结构和生物医学文献。该在线资源检索器可以使用任何计算机语言(Perl,Python,Java和C ++等)将URL发送到应用程序服务器并解析响应。
1.2 注意事项
-
最小化请求数
- 如果任务需要搜索和/或下载大量记录,则使用Entrez历史记录批量上载和/或检索这些记录而不是对每条记录使用单独的请求会更有效
- 可以使用单个EPost请求上载数千个ID
- 可以使用一个EFetch请求下载数百个记录
-
访问限制
- 为了不使服务器过载,NCBI建议用户每秒发布不超过三个URL请求
- 将大型作业限制在工作日的周末或东部时间晚上9:00到凌晨5:00之间
-
设置邮箱
- 使用email参数,这样如果遇到什么问题,NCBI可以通过邮件联系到你
- 邮件的参数从2010年6月1日是强制的参数,所以每次必须告诉 NCBI 是谁在访问
-
URL字符处理
- 所有参数使用小写字符
- 参数没有必需的顺序,通常会忽略空值或不适当的参数
- 避免在URL中放置空格,尤其是在查询中。如果需要空格,请使用加号(+)代替空格
- 其他特殊字符(例如引号(“)或用于引用历史记录服务器上的查询键的#符号)应由其URL编码表示(%22表示”;%23表示#)
二、基本操作
2.1 参数设置
# =====一般参数设置=====
# 设置 email 参数,为了方便 NCBI 的工作人员可以联系到你
# 邮件的参数从2010年6月1日是强制的参数,所以每次必须告诉 NCBI 是谁在访问
Entrez.email = "example@163.com"
# 如果你是通过其他脚本调用,可以设定 tool 的名字,默认为 `Biopython`
Entrez.tool = "exampleScript"
# 可选参数,使用代理,一般在无法正常访问时设置
os.environ["http_proxy"] = "http://proxyhost.example.com:8080"
2.2 查看概况
2.2.1 查看目前 NCBI 所有数据库
from Bio import Entrez
# =====查看数据库概况=====
# 获取 Entrez 所有数据库的句柄
hd_info = Entrez.einfo()
# 获取所有数据库列表
read_info = Entrez.read(hd_info)
for db in read_info['DbList']:
print (db)
2.2.2 查看单个数据库概况
from Bio import Entrez
# 获取 Entrez 的 gene 数据库句柄
hd_info_gene = Entrez.einfo(db="gene")
read_info_gene = Entrez.read(hd_info_gene)
# 数据库名
print ("DbName : ", read_info_gene["DbInfo"]["DbName"])
# 在 NCBI 首页顶部下拉菜单栏中的命名
print ("MenuName : ", read_info_gene["DbInfo"]["MenuName"])
# 数据库描述
print ("Description: ", read_info_gene["DbInfo"]["Description"])
# 数据库收录总数
print ("Count : ", read_info_gene["DbInfo"]["Count"])
# 最新更新时间
print ("LastUpdate : ", read_info_gene["DbInfo"]["LastUpdate"])
# Gene 数据库中可用的搜索关键字列表
print ("FieldList : ", read_info_gene["DbInfo"]["FieldList"])
# 我们把它遍历下
for field in read_info_gene["DbInfo"]["FieldList"]:
print("%(Name)s\t %(FullName)s\t %(Description)s" % field)
2.3 查询
2.3.1 全局搜索 | EGQuery
EGQuery:https://www.ncbi.nlm.nih.gov/books/NBK25499/#chapter4.ESummary
这里只关注搜索关键字在数据库中所有的个数,而不关注它的具体内容
在实际使用中我们可以通过在这里得到的数字来确定下载策略
from Bio import Entrez
# =====全局搜索=====
hd_egquery = Entrez.egquery(term="oct4")
read_egquery = Entrez.read(hd_egquery)
print(read_egquery)
for ele in read_egquery["eGQueryResult"]:
print (ele["DbName"], ele["Count"], ele["Status"])
2.3.2 查询单个数据库中的基因 | ESearch
关于 ESearch 的官方文档 https://www.ncbi.nlm.nih.gov/books/NBK25499/#chapter4.ESearch
from Bio import Entrez
# =====在数据库搜索基因=====
# 搜索 Xenopus laevis 物种中名为 oct4 的基因
handle = Entrez.esearch(db="gene", term="oct4[Gene] AND Xenopus laevis[ORGN]")
read_gene = Entrez.read(handle)
print(read_gene)
2.3.3 查询基因详细描述信息 | Esummary
ESummary :https://www.ncbi.nlm.nih.gov/books/NBK25499/#chapter4.ESummary
from Bio import Entrez
# =====获取摘要=====
# 通过 id 来获取 item 的详细信息
hd_esummary = Entrez.esummary(db="gene", id="397784")
read_esummary = Entrez.read(hd_esummary)
# 获取该基因的详细描述
for key, value in read_esummary['DocumentSummarySet']['DocumentSummary'][0].items():
print(key, value)
2.3.4 查询交叉引用条目 | Elink
Elink:https://www.ncbi.nlm.nih.gov/books/NBK25499/#chapter4.ELink
from Bio import Entrez
# =====搜索交叉引用条目=====
# 接下来我们看看 id 为 5460 的基因相关的文献资料
read_elink = Entrez.read(Entrez.elink(dbfrom="gene", db="pubmed", id="5460"))
print ("LinkSetDb: ", read_elink[0]["LinkSetDb"])
# 查看所有相关的目标库
for lsd in read_elink[0]["LinkSetDb"]:
print (lsd["DbTo"], lsd["LinkName"], len(lsd["Link"]))
# 查看相关的所有文献 Id
for link in read_elink[0]["LinkSetDb"][0]["Link"]:
print (link["Id"])
2.4 其他操作
2.4.1 上传id列表到服务器 | EPost
EPost :https://www.ncbi.nlm.nih.gov/books/NBK25499/#chapter4.EPost
为什么要上传列表到服务器?
你要上传的 id 的列表会以 url 的形式上传到服务器,这里有一个问题,如果 id 很多,就会导致url很长。但是在 HTTP 的协议中,上传一般以 GET 形式,这种方式会限制 url 的长度,也就是说如果用户上传的 URL 太长就会只能局限在一定的长度内,而不能完整的上传到服务器。 为了解决这个问题,只能使用 POST 方式上传,它没有限制文本长度,随后以 HTTP 头文件的形式上传服务器,并以历史记录的形式存储在服务器
from Bio import Entrez
# =====上传历史记录=====
# EPost :https://www.ncbi.nlm.nih.gov/books/NBK25499/#chapter4.EPost
#
id_list = ['379522', '397784', '398336']
hd_epost = Entrez.epost("gene", id=",".join(id_list))
read_epost = Entrez.read(hd_epost)
print ("web_env: ", read_epost["WebEnv"])
print ("query_key: ", read_epost["QueryKey"])
2.4.2 拼写建议与纠正 | Espell
这个模块用来纠正输入的查询词条
from Bio import Entrez
# =====拼写建议=====
hd_espell = Entrez.espell(term="steem cell")
read_espell = Entrez.read(hd_espell)
print ("Query: ", read_espell["Query"])
print ("CorrectedQuery: ", read_espell["CorrectedQuery"])
2.4.3 下载 | EFetch
EFetch:https://www.ncbi.nlm.nih.gov/books/NBK25499/#chapter4.EFetch
所有参数组合:https://www.ncbi.nlm.nih.gov/books/NBK25499/table/chapter4.T._valid_values_of__retmode_and/?report=objectonly
from Bio import Entrez
# =====下载=====
hd_efetch_gb = Entrez.efetch(db='nucleotide', id="397784", rettype='gb', retmode='text')
hd_efetch_fa = Entrez.efetch(db='nucleotide', id="397784", rettype='fasta')
print (hd_efetch_gb.read())
print (hd_efetch_fa.read())
with open("res/397784.txt", "w") as file:
hd_efetch_ml = Entrez.efetch(db='pubmed', id="397784", rettype='medline', retmode='text')
file.write(hd_efetch_ml.read())
with open("res/397784.txt") as file:
read_medline = Medline.read(file)
print ("PMID", read_medline["PMID"])
print ("TI", read_medline["TI"])
2.4.4 解析大文件| parse
一般在 NCBI 中的资源会有较大的内存占用,
这里的parse使用迭代器的方式,而不是像列表全部加载,因此了避免了大文件读取时占满内存
- Linux 系统下准备工作
下载实例文件:ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/ASN_BINARY/Mammalia/Homo_sapiens.ags.gz
下载格式转换工具:ftp://ftp.ncbi.nlm.nih.gov/toolbox/ncbi_tools/converters/by_program/gene2xml/linux64.gene2xml.gz
- 在终端依次运行下列命令
mkdir ncbi
cd ncbi
mkdir ags
mkdir tool
cd tool
wget ftp://ftp.ncbi.nlm.nih.gov/toolbox/ncbi_tools/converters/by_program/gene2xml/linux64.gene2xml.gz
gunzip linux64.gene2xml.gz
mv linux64.gene2xml gene2xml
cd ../ags
wget ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/ASN_BINARY/Mammalia/Homo_sapiens.ags.gz
gunzip Homo_sapiens.ags.gz
../tool/gene2xml -b T -i Homo_sapiens.ags -o Homo_sapiens.xml
- 下载你的目录结构类似这样,这里的Homo_sapiens.xml 大约有15G(2018.09)
- 使用 BioPython 解析
from Bio import Entrez
# =====解析大文件=====
hd_parse = open("Homo_sapiens.xml")
res_parse = Entrez.parse(hd_parse)
for record in res_parse:
status = record['Entrezgene_track-info']['Gene-track']['Gene-track_status']
if status.attributes['value']=='discontinued':
continue
geneid = record['Entrezgene_track-info']['Gene-track']['Gene-track_geneid']
genename = record['Entrezgene_gene']['Gene-ref']['Gene-ref_locus']
print (geneid, genename)
网友评论