这两天处理数据的时候,需要解析Blast的结果。经管之前写过解析blast的xml结果的程序,用起来还是不够优雅。刚好,借这次机会来接触一下biopython。

一、简介

biopython是基于python语言来帮助生物信息学工作这解决问题的工具。它可以做:

1.将生物信息学文件解析为Python可用的数据结构,包含以下支持的格式:

    •  Blast输出结果 – standalone和在线Blast
    • Clustalw
    • FASTA
    • GenBank
    • PubMed和Medline
    • ExPASy文件, 如Enzyme和Prosite
    • SCOP, 包括‘dom’和‘lin’文件
    • UniGene
    • SwissProt

2.被支持格式的文件可以通过记录来重复或者通过字典界面来索引。
3.处理常见的生物信息学在线数据库的代码:

  • NCBI – Blast, Entrez和PubMed服务
  • ExPASy – Swiss-Prot和Prosite条目, 包括Prosite搜索

4.常见生物信息学程序的接口,例如:

  •  NCBI的Standalone
  • Clustalw比对
  • EMBOSS命令行工具

5.一个能处理序列、ID和序列特征的标准序列类。
6.对序列实现常规操作的工具,如翻译,转录和权重计算。
7.利用k最近邻接、Bayes或SVM对数据进行分类的代码。
8.处理比对的代码,包括创建和处理替换矩阵的标准方法。
9.分发并行任务到不同进程的代码。
10.实现序列的基本操作,翻译以及BLAST等功能的GUI程序。
11.使用这些模块的详细文档和帮助,包括此文件,在线的wiki文档,网站和邮件列表。
12.整合BioSQL,一个也被BioPerl和BioJava支持的数据库架构。

感觉好厉害的样子,还是从解析xml开始学起来吧。

biopython下载地址:http://biopython.org/wiki/Download

biopython官网说明:http://biopython-cn.readthedocs.io/zh_CN/latest/index.html

二、实战

1.解析xml文件
blastrecord

from Bio.Blast import NCBIXML
def test_bio_python():
    result_file = "tblastn_100000_2.xml"
    output_file = "tblastn_100000_gene.tsv"

    # with open(output_file,"w") as data1:
    data2 = open(output_file, "w")
    data2.write("seq\tlength\tgen\tmap_len\tmismatch\n")

    result_handle = open(result_file)
    blast_records = NCBIXML.parse(result_handle)
    blast_records = list(blast_records)

    for blast_record in blast_records:
        # blast_record为一条序列的比对结果
        # for description in blast_record.descriptions:
        #     print description.title
        #     print description.score
        #     print description.e
        #     print description.num_alignments
        # break
        # print blast_record.
        # print blast_record.description.title

        # print dir(blast_record) 可以看到每条blast_record的基本属性
        query_id = blast_record.query  # query的名字
        query_length = blast_record.query_length
        # print blast_record.reference
        # print blast_record.query_letters
        # print blast_record.query_id

        for aligment in blast_record.alignments: 
        #aligment为一条序列比对到一个基因,其中一个基因有可能被多次Hits
            for hsp in aligment.hsps:
                # print hsp.__getattr__

                # print aligment.title #被Hit的序列的名字,相当于xml中的
                #<Hit_def>Homo sapiens ER membrane protein complex subunit 1 (EMC1), transcript variant 2, mRNA</Hit_def>

                #获得()中的基因名字
                if "(" in aligment.title:
                    if aligment.title.count("(") == 1:
                        gene = aligment.title.split("(")[1].split(")")[0]
                    elif aligment.title.count("(") == 2:
                        gene = aligment.title.split(")")[1].split("(")[1]
                else:
                    gene = aligment.title

                # print aligment.length   #比对序列的长度


                # print hsp.score
                # print hsp.bits
                # print hsp.expect
                # print hsp.num_alignments
                # print hsp.identities
                identity_len = hsp.identities  #比对到的碱基个数,相当于<Hsp_identity>4</Hsp_identity>
             
                # print hsp.positives  #相当于<Hsp_positive>6</Hsp_positive>
                # print hsp.gaps  #相当于<Hsp_gaps>0</Hsp_gaps>
                # print hsp.strand
                # print hsp.frame
                # print hsp.query #相当于<Hsp_qseq>LLVDDQ</Hsp_qseq>
                # print hsp.query_start
                # print hsp.match  #相当于<Hsp_midline>LL+DD+</Hsp_midline>
                # print hsp.sbjct  #相当于<Hsp_hseq>LLIDDE</Hsp_hseq>
                # print hsp.sbjct_start

                gaps = hsp.gaps
                identity_len = hsp.identities

                align_len = len(query_seq) - int(gaps)
                identity_1 = int(identity_len) / int(align_len) * 100 #比对上的序列的相似性
                identity_2 = int(identity_len) / int(query_length) * 100 #全长比对的百分比
                mismatch = int(query_length) - int(identity_len)
                result = [query_id, str(query_length), gene, str(identity_len),
                          str(mismatch)]
                data2.write("%s\n" % "\t".join(result))
                # break
    data2.close()

后面用到,再来接着补充吧

 

One thought on “biopython实践

发表评论

电子邮件地址不会被公开。 必填项已用*标注