利用Biopython来解析KEGG记录

Biopython

更好的阅读体验请移步这里。

简介

Biopython项目是旨在减少计算生物学中代码重复的开源项目之一，由国际开发人员协会创建。它包含表示生物序列和序列注释的类，并且能够读取和写入各种文件格式（FASTA，FASTQ，GenBank和Clustal等），支持以程序化方式访问生物信息的在线数据库（例如，NCBI）。独立的模块扩展了Biopython的序列比对，蛋白质结构，群体遗传学，系统发育，序列基序和机器学习等功能。

Biopython中的Bio.KEGG模块可以对KEGG数据库中的条目记录进行解析和对KEGG的API进行查询，这对于批量处理KEGG条目有很大帮助。

解析KEGG记录

解析KEGG记录就像利用Biopython解析任何其他文件格式类似。（在运行以下代码之前，请首先将http://rest.kegg.jp/get/ec:5.4.2.2处的文件另存为ec_5.4.2.2.txt。）

>>> from Bio.KEGG import Enzyme
>>> records = Enzyme.parse(open("ec_5.4.2.2.txt"))
>>> record = list(records)[0]
>>> record.classname
['Isomerases;', 'Intramolecular transferases;', 'Phosphotransferases (phosphomutases)']
>>> record.entry
'5.4.2.2'

如果输入的KEGG记录只有一条，也可以使用read方法：

>>> from Bio.KEGG import Enzyme
>>> record = Enzyme.read(open("ec_5.4.2.2.txt"))
>>> record.classname
['Isomerases;', 'Intramolecular transferases;', 'Phosphotransferases (phosphomutases)']
>>> record.entry
'5.4.2.2'

查询KEGG的API

Biopython支持查询KEGG api。支持查询所有KEGG端口；支持KEGG（https://www.kegg.jp/kegg/rest/keggapi.html）条目记录查询的所有方法。该界面对查询进行了一些验证，这些查询遵循KEGG网站上的规则。

下面代码展示了如何通过下载上述相关的酶来扩展示例数据，并将结果传递给解析函数：

>>> from Bio.KEGG import REST
>>> from Bio.KEGG import Enzyme
>>> request = REST.kegg_get("ec:5.4.2.2")
>>> open("ec_5.4.2.2.txt", "w").write(request.read())
>>> records = Enzyme.parse(open("ec_5.4.2.2.txt"))
>>> record = list(records)[0]
>>> record.classname
['Isomerases;', 'Intramolecular transferases;', 'Phosphotransferases (phosphomutases)']
>>> record.entry
'5.4.2.2'

再来一个更实际的例子：

from Bio.KEGG import REST

human_pathways = REST.kegg_list("pathway", "hsa").read()

# 从所有人类相关的通路中过滤出与修复相关的通路
repair_pathways = []
for line in human_pathways.rstrip().split("\n"):
    entry, description = line.split("\t")
    if "repair" in description:
        repair_pathways.append(entry)

# 得到与修复通路相关的基因并将其加入列表
repair_genes = [] 
for pathway in repair_pathways:
    pathway_file = REST.kegg_get(pathway).read()  # 查询并读取每一条通路

    # 遍历每个KEGG通路文件，只读取每个通路中的基因
    current_section = None
    for line in pathway_file.rstrip().split("\n"):
        section = line[:12].strip()  # 12列以内部分的名字
        if not section == "":
            current_section = section
        
        if current_section == "GENE":
            gene_identifiers, gene_description = line[12:].split("; ")
            gene_id, gene_symbol = gene_identifiers.split()

            if not gene_symbol in repair_genes:
                repair_genes.append(gene_symbol)

print("There are %d repair pathways and %d repair genes. The genes are:" % \
      (len(repair_pathways), len(repair_genes)))
print(", ".join(repair_genes))

关于KEGG API的详细说明可以移步其官方说明文档。

参考

Biopython Tutorial and Cookbook

其他资源

Biopython 教程与手册