# Path: kleborate/blastn.py
# def run_blastn(db, query, min_cov, min_ident):
#     build_blast_database_if_needed(db)
# 
#     cmd = 'blastn -task blastn -db {} -query {}'.format(db, query)
#     cmd += " -outfmt '6 sacc pident slen length bitscore qseq sstrand sstart send" \
#            " qacc qstart qend qframe'"
#     cmd += ' -dust no -evalue 1E-20 -word_size 32 -max_target_seqs 10000'
#     cmd += ' -perc_identity {}'.format(min_ident)
# 
#     # TODO: switch this over to subprocess
#     blast_hits = []
#     f = os.popen(cmd)
#     for line in f:
#         blast_hits.append(BlastHit(line))
#     f.close()
# 
#     # Toss out low identity and low coverage hits.
#     if min_ident is not None:
#         blast_hits = [h for h in blast_hits if h.pcid * 100 >= min_ident]
#     if min_cov is not None:
#         blast_hits = [h for h in blast_hits if h.ref_cov * 100 >= min_cov]
# 
#     # Clean up redundant hits so we only have one hit for each part of the genome.
#     blast_hits = cull_redundant_hits(blast_hits)
# 
#     # Get rid of any hits that start with 'delete_'. This is used in cases where there can be
#     # confounding homologous hits. E.g. searching for rmpA2 alleles but not wanting to be misled by
#     # homologous rmpA alleles.
#     blast_hits = [h for h in blast_hits if not h.gene_id.startswith('delete_')]
# 
#     return blast_hits
# 
# Path: kleborate/truncation.py
# def truncation_check(hit, cov_threshold=90.0):
#     """
#     Checks to see if the gene is truncated at the amino acid level.
#     """
#     nucl_seq, ref_start, _ = hit.get_seq_start_end_pos_strand()
# 
#     # The hit must start at the first base of the gene. If not, the gene is considered 0%.
#     if ref_start != 1:
#         return '-0%', 0.0, ''
# 
#     # If there are any ambiguous bases in the sequence, then they will break translation, probably
#     # resulting in truncation call.
#     ambiguous_bases = set(b for b in nucl_seq) - {'A', 'C', 'G', 'T'}
#     for b in ambiguous_bases:
#         nucl_seq = nucl_seq.split(b)[0]
# 
#     # BioPython doesn't like it if the sequence isn't a multiple of 3.
#     nucl_seq = nucl_seq[:len(nucl_seq) // 3 * 3]
# 
#     # The assumption is that the reference allele is a full CDS with a stop codon at the end. This
#     # isn't always true (the reference sequence is sometimes broken) but will serve to make our
#     # denominator for coverage.
#     ref_aa_length = (hit.ref_length - 3) // 3
# 
#     coding_dna = Seq(nucl_seq)
#     translation = str(coding_dna.translate(table='Bacterial', to_stop=True))
# 
#     coverage = 100.0 * len(translation) / ref_aa_length
#     if coverage >= cov_threshold:
#         return '', coverage, translation
#     else:
#         return '-{:.0f}%'.format(coverage), coverage, translation


import collections
import os
import re
from .blastn import run_blastn
from .truncation import truncation_check
        locus = gene_id_components[1]
        allele = gene_id_components[2]
    else:
        allele = hit.gene_id
        locus = hit.gene_id.split('_')[0]
    return allele, locus


def keep_only_one_hit_per_locus(hits):
    hits_per_locus = collections.defaultdict(list)
    for hit in hits:
        _, locus = get_allele_and_locus(hit)
        hits_per_locus[locus].append(hit)
    kept_hits = []
    for locus, locus_hits in hits_per_locus.items():
        locus_hits = sorted(locus_hits, key=lambda x: x.score, reverse=True)  # sort best to worst
        best_hit = locus_hits[0]
        kept_hits.append(best_hit)
    return kept_hits


def get_best_allele_per_locus(hits, check_for_truncation):
    best_scores = {}   # key = locus, value = BLAST score for best allele encountered so far
    best_alleles = {}  # key = locus, value = best allele (* if imprecise match)

    for hit in hits:
        allele, locus = get_allele_and_locus(hit)
        if hit.pcid < 100.00 or hit.alignment_length < hit.ref_length:
            allele += '*'  # inexact match
        if check_for_truncation:

# Path: kleborate/blastn.py
# def run_blastn(db, query, min_cov, min_ident):
#     build_blast_database_if_needed(db)
# 
#     cmd = 'blastn -task blastn -db {} -query {}'.format(db, query)
#     cmd += " -outfmt '6 sacc pident slen length bitscore qseq sstrand sstart send" \
#            " qacc qstart qend qframe'"
#     cmd += ' -dust no -evalue 1E-20 -word_size 32 -max_target_seqs 10000'
#     cmd += ' -perc_identity {}'.format(min_ident)
# 
#     # TODO: switch this over to subprocess
#     blast_hits = []
#     f = os.popen(cmd)
#     for line in f:
#         blast_hits.append(BlastHit(line))
#     f.close()
# 
#     # Toss out low identity and low coverage hits.
#     if min_ident is not None:
#         blast_hits = [h for h in blast_hits if h.pcid * 100 >= min_ident]
#     if min_cov is not None:
#         blast_hits = [h for h in blast_hits if h.ref_cov * 100 >= min_cov]
# 
#     # Clean up redundant hits so we only have one hit for each part of the genome.
#     blast_hits = cull_redundant_hits(blast_hits)
# 
#     # Get rid of any hits that start with 'delete_'. This is used in cases where there can be
#     # confounding homologous hits. E.g. searching for rmpA2 alleles but not wanting to be misled by
#     # homologous rmpA alleles.
#     blast_hits = [h for h in blast_hits if not h.gene_id.startswith('delete_')]
# 
#     return blast_hits
# 
# Path: kleborate/truncation.py
# def truncation_check(hit, cov_threshold=90.0):
#     """
#     Checks to see if the gene is truncated at the amino acid level.
#     """
#     nucl_seq, ref_start, _ = hit.get_seq_start_end_pos_strand()
# 
#     # The hit must start at the first base of the gene. If not, the gene is considered 0%.
#     if ref_start != 1:
#         return '-0%', 0.0, ''
# 
#     # If there are any ambiguous bases in the sequence, then they will break translation, probably
#     # resulting in truncation call.
#     ambiguous_bases = set(b for b in nucl_seq) - {'A', 'C', 'G', 'T'}
#     for b in ambiguous_bases:
#         nucl_seq = nucl_seq.split(b)[0]
# 
#     # BioPython doesn't like it if the sequence isn't a multiple of 3.
#     nucl_seq = nucl_seq[:len(nucl_seq) // 3 * 3]
# 
#     # The assumption is that the reference allele is a full CDS with a stop codon at the end. This
#     # isn't always true (the reference sequence is sometimes broken) but will serve to make our
#     # denominator for coverage.
#     ref_aa_length = (hit.ref_length - 3) // 3
# 
#     coding_dna = Seq(nucl_seq)
#     translation = str(coding_dna.translate(table='Bacterial', to_stop=True))
# 
#     coverage = 100.0 * len(translation) / ref_aa_length
#     if coverage >= cov_threshold:
#         return '', coverage, translation
#     else:
#         return '-{:.0f}%'.format(coverage), coverage, translation


# Path: kleborate/mlstBLAST.py
import collections
import os
import re
from .blastn import run_blastn
from .truncation import truncation_check

        locus = gene_id_components[1]
        allele = gene_id_components[2]
    else:
        allele = hit.gene_id
        locus = hit.gene_id.split('_')[0]
    return allele, locus


def keep_only_one_hit_per_locus(hits):
    hits_per_locus = collections.defaultdict(list)
    for hit in hits:
        _, locus = get_allele_and_locus(hit)
        hits_per_locus[locus].append(hit)
    kept_hits = []
    for locus, locus_hits in hits_per_locus.items():
        locus_hits = sorted(locus_hits, key=lambda x: x.score, reverse=True)  # sort best to worst
        best_hit = locus_hits[0]
        kept_hits.append(best_hit)
    return kept_hits


def get_best_allele_per_locus(hits, check_for_truncation):
    best_scores = {}   # key = locus, value = BLAST score for best allele encountered so far
    best_alleles = {}  # key = locus, value = best allele (* if imprecise match)

    for hit in hits:
        allele, locus = get_allele_and_locus(hit)
        if hit.pcid < 100.00 or hit.alignment_length < hit.ref_length:
            allele += '*'  # inexact match
        if check_for_truncation:

            allele += truncation_check(hit)[0]