Gets the a list of alleles that are the most similar to the input sequence

    :param sequences: The input sequence record.
    :type sequences: SeqRecord
    :param locus: The gene locus associated with the sequence.
    :type locus: ``str``
    :param nseqs: The incomplete annotation from a previous iteration.
    :type nseqs: ``int``
    :param evalue: The evalue to use (default = 10)
    :type evalue: ``int``
    :param kir: Run with KIR or not
    :type kir: ``bool``
    :param verbose: Run in versboe
    :type verbose: ``bool``
    :param refdata: An object with reference data
    :type refdata: :ref:`ref`
    :rtype: :ref:`bl`

    Example usage:

        >>> from Bio.Seq import Seq
        >>> from seqann.blast_cmd import blastn
        >>> sequence = Seq('AGAGACTCTCCCGAGGATTTCGTGTACCAGTTTAAGGCCATGTGCTACTTCACC')
        >>> blast = blastn(sequence, locus, nseqs)

def blastn(sequences, locus, nseqs, kir=False,
           verbose=False, refdata=None, evalue=10):
    """
    Gets the a list of alleles that are the most similar to the input sequence

    :param sequences: The input sequence record.
    :type sequences: SeqRecord
    :param locus: The gene locus associated with the sequence.
    :type locus: ``str``
    :param nseqs: The incomplete annotation from a previous iteration.
    :type nseqs: ``int``
    :param evalue: The evalue to use (default = 10)
    :type evalue: ``int``
    :param kir: Run with KIR or not
    :type kir: ``bool``
    :param verbose: Run in versboe
    :type verbose: ``bool``
    :param refdata: An object with reference data
    :type refdata: :ref:`ref`
    :rtype: :ref:`bl`

    Example usage:

        >>> from Bio.Seq import Seq
        >>> from seqann.blast_cmd import blastn
        >>> sequence = Seq('AGAGACTCTCCCGAGGATTTCGTGTACCAGTTTAAGGCCATGTGCTACTTCACC')
        >>> blast = blastn(sequence, locus, nseqs)

    """
    logger = logging.getLogger("Logger." + __name__)

    if not refdata:
        refdata = ReferenceData()

    file_id = str(randomid())
    input_fasta = file_id + ".fasta"
    output_xml = file_id + ".xml"
    SeqIO.write(sequences, input_fasta, "fasta")
    blastn_cline = NcbiblastnCommandline(query=input_fasta,
                                         db=refdata.blastdb,
                                         evalue=evalue,
                                         outfmt=5,
                                         reward=1,
                                         penalty=-3,
                                         gapopen=5,
                                         gapextend=2,
                                         dust='yes',
                                         out=output_xml)

    stdout, stderr = blastn_cline()
    loc = locus
    if not kir:
        loc = locus.split("-")[1]

    blast_qresult = SearchIO.read(output_xml, 'blast-xml')

    #   Delete files
    cleanup(file_id)

    # TODO: Use logging
    if len(blast_qresult.hits) == 0:
        if verbose:
            logger.error("Failed blast! No hits!")
            logger.error(stderr)
        return Blast(failed=True)

    alleles = []
    full_sequences = []
    load_blast = 70 if nseqs < 70 else nseqs
    l = len(blast_qresult.hits) if load_blast > len(blast_qresult.hits) else load_blast

    # TODO: update all blast files to have HLA-
    if locus in refdata.hla_loci and not kir:
        alleles = [blast_qresult[i].id.split("_")[0] for i in range(0, l)
                   if blast_qresult[i].id.split("*")[0] == locus or "HLA-" + blast_qresult[i].id.split("*")[0] == locus]
        alleles = ["HLA-" + a if not has_hla(a) else a for a in alleles]
    if kir:
        alleles = [blast_qresult[i].id.split("_")[0] for i in range(0, l)
                   if blast_qresult[i].id.split("*")[0] == locus]

    if verbose:
        logger.info("Blast alleles: " + ",".join(alleles))

    # TODO: sort alleles by number of features they contain and evalue
    # Use biosql db if provided
    # otherwise use IMGT dat file
    final_seqs = []
    rmax = refdata.structure_max[locus]
    if refdata.server_avail:
        db = refdata.server[refdata.dbversion + "_" + loc]
        full_sequences = []
        for n in alleles:
            if n in refdata.hla_names:
                try:
                    seq = db.lookup(name=n)
                    full_sequences.append(seq)
                except:
                    logger.error("Allele doesnt exist in IMGT BioSQL DB!! "
                                 + n)
    else:
        if verbose:
            logger.info("Getting sequences from HLA.dat file")

        full_sequences = [refdata.hlaref[a] for a in alleles
                          if a in refdata.hlaref]

        for s in full_sequences:
            s.name = s.description.split(",")[0]

    i = 1
    last_seq = []
    max_f = 0
    added_max = False
    full_feats = False
    for s in full_sequences:

        fs = len([f.type for f in s.features
                  if not f.type in ['source', 'CDS']])
        if i <= nseqs:
            final_seqs.append(s)
            max_f = fs if fs > max_f else max_f

        if i <= nseqs and max_f < rmax:
            full_feats = True

        if(i >= nseqs and fs == max_f and not added_max):
            if len(last_seq) >= 10:
                last_seq.insert(3, s)
            else:
                last_seq.append(s)
            added_max = True

        if(fs > max_f and len(last_seq) < 10
           and i >= nseqs and len(last_seq) < 10):
            last_seq.append(s)

        i += 1

    if full_feats:
        for s in last_seq:
            final_seqs.append(s)

    #   Build Blast object
    blast_o = Blast(match_seqs=final_seqs, alleles=alleles)
    return blast_o