生物信息中的Python 03 | 自动化操作NCBI

相信大家在上一文中下载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"] = # =====一般参数设置=====
# 设置 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']:
    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(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[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)
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():
    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[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[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[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[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）

mark

使用 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']
     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)