import subprocess as sb
from sys import argv
from statistics import mean

genomename = argv[1]
leng = argv[2]

sb.call(['makeblastdb', 
         '-in', genomename, 
         '-dbtype', 'nucl',
         '-out', 'bacterium_genome_db',
         '-logfile', 'bacterium_genome_db_log'])
sb.call(['makenucseq', 'fasta::100randseqs.fasta', 
         '-amount', '100000', 
         '-length', leng,
         '-filter'])
sb.call(['blastn', 
         '-task', 'blastn', 
         '-query', '100randseqs.fasta', 
         '-db', 'bacterium_genome_db', 
         '-outfmt', '6 evalue',
         '-out', 'output.tsv'])

evaluecount = []
tsvfile = open('output.tsv', 'r')
for line in tsvfile:
    evalue = float(line.strip())
    if evalue < 0.1:
        evaluecount.append(1)
    else:
        evaluecount.append(0)
tsvfile.close()
    
print('В среднем blastn находит {:.2f} находок c E-value < 0.1 для случайной последовательности длины {} в Вашем геноме бактерии.'.format(mean(evaluecount), leng))