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 等相关信息,最终将结果输出。
import sys
infile = sys.argv[1]
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]
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 或者在终端中通过指定输入文件或输出文件来运行。
|