from Bio import SeqIO
from Bio.Blast import NCBIWWW
from Bio.Blast import NCBIXML

E_VALUE_THRESH = 0.04

def print_seqs():
    for seq_record in SeqIO.parse("ls_orchid.gbk", "genbank"):
        print(seq_record.id)
        print(repr(seq_record.seq))
        print(len(seq_record))
    
def convert():
    records = (rec for rec in SeqIO.parse("ls_orchid.gbk", "genbank"))
    SeqIO.write(records, "ls_orchid.fasta", "fasta")

def blast():   
    record = SeqIO.read("m_cold.fasta", format="fasta")
    result_handle = NCBIWWW.qblast("blastn", "nt", record.format("fasta"))
    with open("blast.xml","w") as blastlog:
        blastlog.write(result_handle.read())

def parse_blast():
    result_handle = open("blast.xml")
    blast_record = NCBIXML.read(result_handle)
    for alignment in blast_record.alignments:
       #print(alignment.title)
       for hsp in alignment.hsps:
            if hsp.expect < E_VALUE_THRESH:
                print('****Alignment****')
                print('sequence:', alignment.title)
                print('length:', alignment.length)
                print('e value:', hsp.expect)
                print(hsp.query[0:75] + '...')
                print(hsp.match[0:75] + '...')
                print(hsp.sbjct[0:75] + '...')
        

print_seqs()
convert()
blast()
parse_blast()