初尝BioPython小记

NCBI官网：https://www.ncbi.nlm.nih.gov/
BioPython官方文档：https://biopython-cn.readthedocs.io/zh_CN/latest/

bioPython是生物信息领域用于基因序列比对的一个工具，在搭建本地化blast的时候很实用。
以下笔记内容是在开发过程中参考官方文档 第七章：BLAST 后的相关问题及解决记录。

首先，使用NCBIApplications编写blast命令行

blastx_cline = NcbiblastnCommandline(cmd='blastn',query='file.txt',db='nt',evalue=0.01,outfmt=5,out='output.xml')
blastx_cline.set_parameter("task",'megablast')
blastx_cline.set_parameter("num_threads",100)
blastx_cline.set_parameter("num_alignments",100)

其中有个小问题，-task参数我一直加不进去，这个参数用于不同的查询方式，即-task megablast/dc-megablast/blastn，运行后根据报错信息，我进了源码查看，结果发现它与众不同……相关的set方法没找到正确的将其置空的位置，于是我心一横，直接改源码（。

改完后正常了，由于目前应用有限，暂时没有bug，这个等待后续更新了。

blast跑一次大概20s，我掐着计时器算的，同一个序列文件不论是megablast还是dc-megablast或者是blastn都是差不多时间，这块大概也需要后续更多的测试。

在跑完blast数据后，下一步就是解析了。从生成的xml文件提取数据，其实一开始没啥头绪，看着那HTML标签我甚至以为要搞爬虫。。不过并不用，事实上，bioPython已经提供了相应的解析工具，而这里我用到了两个：SearchIO和NCBIXML。

results_io = SearchIO.parse('output.xml',"blast-xml")
for result in results_io:
    seq_map_id = "".join([result.id ,";"])
    for hsp in result.hits:
        hsp_id = hsp.id
        hsp_description = hsp.description

之所以使用两个解析器，是因为NCBIXML解析出来的list并没有按照query id罗列。在使用blast比对序列的时候使用的fastq文件内其实很经常会有多个序列，如果没有query id帮忙区分，会造成很大的困扰。相当于一排没有标签的书架，我知道我要的书就在里面，但我没有指引，只能一排排书架去找，这很浪费时间。所以我先用SearchIO取出query id和description。其实可以的话，我应该也能够在SearchIO里找到我需要的属性元素，然鹅大脑当机没看源码。。也没找到详细的UML图给我取用他的元素。这个后面再研究下。

然后再用NCBIXML解析，前面说到的UML图，其实BioPython官方文档里的也有点小坑，不过不妨碍事，源码一看即得x（所以说上面的SearchIO也完全可以看源码取属性的。。大意了。。

result_handle = open('output.xml')
blast_records = NCBIXML.parse(result_handle)
for blast_record in blast_records:
    for alignment in blast_record.alignments:
        title = alignment.title.split(" ")[0]
        length = alignment.length
        accession = alignment.accession
        hit_def = alignment.hit_def
        hit_id = alignment.hit_id
        for hsp in alignment.hsps:
            bits = str('%.1f' % hsp.bits)
            score = str('%.1f' % hsp.score)
            identities = hsp.identities
            positives = str(hsp.positives)
            gaps = str(hsp.gaps)
            strand = str(hsp.strand)
            frame = str(hsp.frame)
            query = str(hsp.query)
            query_start = str(hsp.query_start)
            query_end = str(hsp.query_end)
            match = hsp.match
            sbjct = str(hsp.sbjct)
            sbjct_start = str(hsp.sbjct_start)
            sbjct_end = str(hsp.sbjct_end)
            num_alignment = str(hsp.num_alignments)
            align_len = hsp.align_length
            # 这个percentage对应NCBI上的identities百分比
            percentage = '%.2f' % (float(identities)/align_len*100)