Blast结果处理

1. 说明

该脚本旨在处理由本地版blast序列比对后结果，具体如下所述：
blast输出格式如下所示(-outfmt = "7 qacc sacc score evalue pident qcovs ppos")：

# BLASTP 2.7.1+
# Query: NP_040533.1
# Database: blastdatabase/ElmMotif.fasta
# 0 hits found
# BLASTP 2.7.1+
# Query: NP_056788.1
# Database: blastdatabase/ElmMotif.fasta
# Fields: query acc., subject acc., score, evalue, % identity, % query coverage per subject, % positives
# 2 hits found
NP_056788.1	PF01115-1023-1041-1666	38	2.4	57.143	2	71.43
NP_056788.1	PF01267-1023-1041-1665	38	2.4	57.143	2	71.43
# BLAST processed 2 queries

2. 目的：

从blast的结果中将有Fields的部分对应的第一行的subject及相关的score，evalue等信息提取出来，最终结果如下图所示：

3. Python3脚本：

该脚本的基本思路是：
首先根据[# BLAST ... # BLAST)将文本划分为多个由# BLAST起始的部分，在根据是否有# Fidelds来决定是否能够获取subject、score等相关信息，最终将结果输出。

#!/usr/bin/python3
import sys
### 脚本中指定输入和输出文件
# infile = file1 
# outfile = file2
### 在终端中指定输入和输出文件
infile = sys.argv[1] # 输入的 Blast 结果形式是 -outfmt "7 qacc sacc score evalue pident qcovs ppos"
outfile  = sys.argv[2]

outF = open(outfile,'w')
outF.write('%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\n' % ('Blast version','Query','Database','The number of \"hits found\"','Subject','Score','Evalue','Identity','Query coverage per subject','Postives'))#

RowNumber = 1
with open(infile) as inF:
    for line in inF:
        if line[0:7] == '# BLAST':
            if RowNumber > 1:
                templine = templine0
                RowNumber = 1
                templine0 = []
                templine0.append(line)
                if templine[3].strip()[0:9] == '# Fields:':
                    blastVersion = templine[0].strip()[2:]
                    query = templine[1].strip()[8:]
                    database = templine[2].strip()[11:]
                    numberHits = templine[4].strip().split()[1]
                    """
                    # 该注释部分是用于提取 Hit 上的所有行，如果需要提取所有结果，只需要把此处注释去掉，下面 subject ... outF.write 的7行注释掉即可。
                    for i in range(5,len(templine)):
                    	subject = templine[i].strip().split('\t')[1]
                        score = templine[i].strip().split('\t')[2]
                        evalue = templine[i].strip().split('\t')[3]
                        identity = templine[i].strip().split('\t')[4]
                        queryPerSubject = templine[i].strip().split('\t')[5]
                        positive = templine[i].strip().split('\t')[6]
                        outF.write('%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\n' % 
                                (blastVersion,query,database,numberHits,subject,score,evalue,identity,queryPerSubject,positive))
                     """
                    subject = templine[5].strip().split('\t')[1]
                    score = templine[5].strip().split('\t')[2]
                    evalue = templine[5].strip().split('\t')[3]
                    identity = templine[5].strip().split('\t')[4]
                    queryPerSubject = templine[5].strip().split('\t')[5]
                    positive = templine[5].strip().split('\t')[6]
                    outF.write('%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\n' % 
                                (blastVersion,query,database,numberHits,subject,score,evalue,identity,queryPerSubject,positive))
                else:
                    continue
            RowNumber = 1
            templine0 = []
            templine0.append(line)
        else:
            RowNumber += 1
            templine0.append(line)

运行方式：直接在脚本中指定infile和outfile或者在终端中通过指定输入文件或输出文件来运行。