选择blast比对第一条best hit的结果

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()

推荐阅读更多精彩内容