blast的结果中每条contig都会有很多hit,最好不用max_target_hit这个参数,因为它输出的是数据库中top的序列而不一定是比对最好的序列,即使所有版本都包含相同的最佳匹配结果,但是BLAST却返回不同的结果。而且以不同的方式对数据库进行排序,也会导致在将max_target_seqs参数设置为1时,BLAST返回不同的“top hit”。因此我们需要自己挑出所有hit中第一条打分最高的hit。
参考:https://www.jianshu.com/p/7eb530bc1a9c
nohup time blastn -db ref_prok_rep_genomes -query ../nextpolish/genome.nextpolish.fasta -perc_identity 90 -num_threads 16 -evalue 1e-5 -outfmt "6 std sscinames qcovus" -out nextpolish_blast_all.txt &
python blast_best.py nextpolish_blast_all.txt nextpolish_blast_best.txt
#!/usr/bin/python
# -*- coding: utf-8 -*-
# conda activate python3
"""
作者:徐诗芬
内容:读取并输出所有blast结果中最好的第一条hit结果
日期:2021.1.22
"""
import sys
def usage():
print('Usage: python blast_best.py [input_file] [outfile]')
def main():
# global name
inf = open(sys.argv[1], 'rt')
ouf = open(sys.argv[2], 'wt')
flag_list = []
for line in inf:
line = line.strip()
name = line.split("\t")[0]
#id = eval(line.split("\t")[2])
if name not in flag_list:
ouf.write(line + '\n')
else:
continue
flag_list.append(name)
inf.close()
ouf.close()
try:
main()
except IndexError:
usage()