annotate - method for annotating a BioPython sequence

        :param sequence: The input consensus sequence.
        :type sequence: Seq
        :param locus: The gene locus associated with the sequence.
        :type locus: ``str``
        :param nseqs: The number of blast sequences to use.
        :type nseqs: ``int``
        :param alignseqs: The number of sequences to use for targeted alignments.
        :type alignseqs: ``int``
        :param skip: A list of alleles to skip for using as a reference. This is used for validation and testing.
        :type skip: ``List``
        :rtype: :ref:`ann`

        Returns:
            The annotate function return an ``Annotation`` object that
            contains the sequence features and names associated with them.

            Example output::

                {
                     'complete_annotation': True,
                     'annotation': {'exon_1': SeqRecord(seq=Seq('AGAGACTCTCCCG', SingleLetterAlphabet()), id='HLA:HLA00630', name='HLA:HLA00630', description='HLA:HLA00630 DQB1*03:04:01 597 bp', dbxrefs=[]),
                                    'exon_2': SeqRecord(seq=Seq('AGGATTTCGTGTACCAGTTTAAGGCCATGTGCTACTTCACCAACGGGACGGAGC...GAG', SingleLetterAlphabet()), id='HLA:HLA00630', name='HLA:HLA00630', description='HLA:HLA00630 DQB1*03:04:01 597 bp', dbxrefs=[]),
                                    'exon_3': SeqRecord(seq=Seq('TGGAGCCCACAGTGACCATCTCCCCATCCAGGACAGAGGCCCTCAACCACCACA...ATG', SingleLetterAlphabet()), id='HLA:HLA00630', name='<unknown name>', description='HLA:HLA00630', dbxrefs=[])},
                     'features': {'exon_1': SeqFeature(FeatureLocation(ExactPosition(0), ExactPosition(13), strand=1), type='exon_1'),
                                  'exon_2': SeqFeature(FeatureLocation(ExactPosition(13), ExactPosition(283), strand=1), type='exon_2')
                                  'exon_3': SeqFeature(FeatureLocation(ExactPosition(283), ExactPosition(503), strand=1), type='exon_3')},
                     'method': 'nt_search and clustalo',
                     'gfe': 'HLA-Aw2-1-1-1-1-1-1-1-1-1-1-1-1-1-1-1-4',
                     'seq': SeqRecord(seq=Seq('AGAGACTCTCCCGAGGATTTCGTGTACCAGTTTAAGGCCATGTGCTACTTCACC...ATG', SingleLetterAlphabet()), id='HLA:HLA00630', name='HLA:HLA00630', description='HLA:HLA00630 DQB1*03:04:01 597 bp', dbxrefs=[])
                }


        Example usage:

            >>> from Bio.Seq import Seq
            >>> from seqann import BioSeqAnn
            >>> sequence = Seq('AGAGACTCTCCCGAGGATTTCGTGTACCAGTTTAAGGCCATGTGCTACTTCACC')
            >>> seqann = BioSeqAnn()
            >>> ann = seqann.annotate(sequence)
            >>> for f in ann.annotation:
            ...    print(f, ann.method, str(ann.annotation[f].seq), sep="\t")
            exon_2  nt_search and clustalo  AGGATTTCGTGTACCAGTTTAAGGCCATGTGCTACTTCACCAACGGGACGGAGCGCGTGCGTTATGTGACCAGATACATCTATAACCGAGAGGAGTACGCACGCTTCGACAGCGACGTGGAGGTGTACCGGGCGGTGACGCCGCTGGGGCCGCCTGCCGCCGAGTACTGGAACAGCCAGAAGGAAGTCCTGGAGAGGACCCGGGCGGAGTTGGACACGGTGTGCAGACACAACTACCAGTTGGAGCTCCGCACGACCTTGCAGCGGCGAG
            exon_3  nt_search and clustalo  TGGAGCCCACAGTGACCATCTCCCCATCCAGGACAGAGGCCCTCAACCACCACAACCTGCTGGTCTGCTCAGTGACAGATTTCTATCCAGCCCAGATCAAAGTCCGGTGGTTTCGGAATGACCAGGAGGAGACAACCGGCGTTGTGTCCACCCCCCTTATTAGGAACGGTGACTGGACCTTCCAGATCCTGGTGATGCTGGAAATGACTCCCCAGCATGGAGACGTCTACACCTGCCACGTGGAGCACCCCAGCCTCCAGAACCCCATCACCGTGGAGTGGC
            exon_1  nt_search and clustalo  AGAGACTCTCCCG
            exon_4  nt_search and clustalo  GGGCTCAGTCTGAATCTGCCCAGAGCAAGATG

def annotate(self, sequence: Seq=None, locus: str=None,
                 nseqs: int=20, alignseqs: int=10,
                 skip: List=[],
                 rerun: bool=True,
                 full: bool=True) -> Annotation:
        """
        annotate - method for annotating a BioPython sequence

        :param sequence: The input consensus sequence.
        :type sequence: Seq
        :param locus: The gene locus associated with the sequence.
        :type locus: ``str``
        :param nseqs: The number of blast sequences to use.
        :type nseqs: ``int``
        :param alignseqs: The number of sequences to use for targeted alignments.
        :type alignseqs: ``int``
        :param skip: A list of alleles to skip for using as a reference. This is used for validation and testing.
        :type skip: ``List``
        :rtype: :ref:`ann`

        Returns:
            The annotate function return an ``Annotation`` object that
            contains the sequence features and names associated with them.

            Example output::

                {
                     'complete_annotation': True,
                     'annotation': {'exon_1': SeqRecord(seq=Seq('AGAGACTCTCCCG', SingleLetterAlphabet()), id='HLA:HLA00630', name='HLA:HLA00630', description='HLA:HLA00630 DQB1*03:04:01 597 bp', dbxrefs=[]),
                                    'exon_2': SeqRecord(seq=Seq('AGGATTTCGTGTACCAGTTTAAGGCCATGTGCTACTTCACCAACGGGACGGAGC...GAG', SingleLetterAlphabet()), id='HLA:HLA00630', name='HLA:HLA00630', description='HLA:HLA00630 DQB1*03:04:01 597 bp', dbxrefs=[]),
                                    'exon_3': SeqRecord(seq=Seq('TGGAGCCCACAGTGACCATCTCCCCATCCAGGACAGAGGCCCTCAACCACCACA...ATG', SingleLetterAlphabet()), id='HLA:HLA00630', name='<unknown name>', description='HLA:HLA00630', dbxrefs=[])},
                     'features': {'exon_1': SeqFeature(FeatureLocation(ExactPosition(0), ExactPosition(13), strand=1), type='exon_1'),
                                  'exon_2': SeqFeature(FeatureLocation(ExactPosition(13), ExactPosition(283), strand=1), type='exon_2')
                                  'exon_3': SeqFeature(FeatureLocation(ExactPosition(283), ExactPosition(503), strand=1), type='exon_3')},
                     'method': 'nt_search and clustalo',
                     'gfe': 'HLA-Aw2-1-1-1-1-1-1-1-1-1-1-1-1-1-1-1-4',
                     'seq': SeqRecord(seq=Seq('AGAGACTCTCCCGAGGATTTCGTGTACCAGTTTAAGGCCATGTGCTACTTCACC...ATG', SingleLetterAlphabet()), id='HLA:HLA00630', name='HLA:HLA00630', description='HLA:HLA00630 DQB1*03:04:01 597 bp', dbxrefs=[])
                }


        Example usage:

            >>> from Bio.Seq import Seq
            >>> from seqann import BioSeqAnn
            >>> sequence = Seq('AGAGACTCTCCCGAGGATTTCGTGTACCAGTTTAAGGCCATGTGCTACTTCACC')
            >>> seqann = BioSeqAnn()
            >>> ann = seqann.annotate(sequence)
            >>> for f in ann.annotation:
            ...    print(f, ann.method, str(ann.annotation[f].seq), sep="\t")
            exon_2  nt_search and clustalo  AGGATTTCGTGTACCAGTTTAAGGCCATGTGCTACTTCACCAACGGGACGGAGCGCGTGCGTTATGTGACCAGATACATCTATAACCGAGAGGAGTACGCACGCTTCGACAGCGACGTGGAGGTGTACCGGGCGGTGACGCCGCTGGGGCCGCCTGCCGCCGAGTACTGGAACAGCCAGAAGGAAGTCCTGGAGAGGACCCGGGCGGAGTTGGACACGGTGTGCAGACACAACTACCAGTTGGAGCTCCGCACGACCTTGCAGCGGCGAG
            exon_3  nt_search and clustalo  TGGAGCCCACAGTGACCATCTCCCCATCCAGGACAGAGGCCCTCAACCACCACAACCTGCTGGTCTGCTCAGTGACAGATTTCTATCCAGCCCAGATCAAAGTCCGGTGGTTTCGGAATGACCAGGAGGAGACAACCGGCGTTGTGTCCACCCCCCTTATTAGGAACGGTGACTGGACCTTCCAGATCCTGGTGATGCTGGAAATGACTCCCCAGCATGGAGACGTCTACACCTGCCACGTGGAGCACCCCAGCCTCCAGAACCCCATCACCGTGGAGTGGC
            exon_1  nt_search and clustalo  AGAGACTCTCCCG
            exon_4  nt_search and clustalo  GGGCTCAGTCTGAATCTGCCCAGAGCAAGATG

        """

        # If sequence is now a biopython
        # sequence record convert it to one
        if isinstance(sequence, Seq):
            sequence = SeqRecord(seq=sequence,
                                 id="NO_ID")

        if isinstance(sequence, str):
            sequence = SeqRecord(seq=Seq(sequence),
                                 id="NO_ID")

        # If sequence contains any characters
        # other than ATCG then the GFE notation
        # can not be created
        create_gfe = checkseq(sequence.seq)

        if self.verbose and not create_gfe:
            self.logger.warning(self.logname + " Sequence alphabet "
                                + "contains non DNA")
            self.logger.warning(self.logname
                                + " No GFE string will be generated")

        # Check it the locus exists
        if not locus:
            if self.verbose:
                self.logger.info(self.logname + " No locus provided! ")

            # Guessing locus with blastn
            locus = get_locus(sequence, kir=self.kir, refdata=self.refdata)

            if locus and self.verbose:
                self.logger.info(self.logname + " Locus prediction = " + locus)

            if not locus:
                if self.verbose:
                    self.logger.error(self.logname
                                      + " Locus could not be determined!")
                # TODO: Raise exception
                #raise NoLocusException("")
                return

        # Exact match found
        matched_annotation = self.refdata.search_refdata(sequence, locus)
        if matched_annotation and not skip:

            matched_annotation.exact = True

            # TODO: return name of allele
            if self.verbose:
                self.logger.info(self.logname + " exact match found")

            # Create GFE
            if create_gfe:
                feats, gfe = self.gfe.get_gfe(matched_annotation, locus)
                matched_annotation.gfe = gfe
                matched_annotation.structure = feats
            return matched_annotation

        # Run blast to get ref sequences
        blast = blastn(sequence, locus, nseqs,
                       kir=self.kir, verbose=self.verbose,
                       refdata=self.refdata)

        # If the blastn fails..
        if blast.failed:

            if self.verbose:
                self.logger.info(self.logname + " Blastn failed!!")

            # Try and determine the locus and rerun. This is
            # useful for cases when the sequences is associated
            # with the wrong locus.
            locus = get_locus(sequence, kir=self.kir, refdata=self.refdata)

            if locus and self.verbose:
                self.logger.info(self.logname + " Locus prediction = " + locus)

            # Check if the locus could not be found
            if not locus:
                self.logger.error(self.logname
                                  + " Locus could not be determined!")
                # TODO: Raise exception
                return
            return self.annotate(sequence, locus,
                                 nseqs=nseqs,
                                 alignseqs=alignseqs)
            return

        # Do seq_search first on all blast sequences
        # retain what the largest ref seq is
        leastmissing = 100
        partial_ann = None
        leastmissing_feat = None
        found = blast.match_seqs
        for i in range(0, len(found)):

            run = 0
            if i == len(found)-1:
                run = 1

            # Skip a reference
            # * For validation *
            if found[i].name in skip:
                continue

            if self.verbose:
                self.logger.info(self.logname
                                 + " Running seq_search with "
                                 + found[i].name)

            # * Running sequence search *
            # This does a simple string search for the
            # reference features within the provided sequence
            ann = self.seqsearch.search_seqs(found[i],
                                             sequence, locus,
                                             partial_ann=partial_ann,
                                             run=run)

            if ann.complete_annotation:
                if self.verbose:
                    self.logger.info(self.logname
                                     + " Finished annotation using "
                                     + found[i].name)

                # Add alignment flag is specified
                if self.align:
                    if self.verbose:
                        self.logger.info(self.logname + " Adding alignment")
                    ann = self.add_alignment(found[i], ann)

                if self.verbose and self.verbosity > 0:
                    self.logger.info(self.logname
                                     + " Features annotated = "
                                     + ",".join(list(ann
                                                     .annotation.keys())))
                    if self.verbosity > 3:
                        for f in ann.features:
                            self.logger.info(self.logname
                                             + " " + f + " = "
                                             + str(ann
                                                   .annotation[f].seq))

                # Create GFE
                if create_gfe:
                    feats, gfe = self.gfe.get_gfe(ann, locus)
                    ann.gfe = gfe
                    ann.structure = feats
                ann.clean()
                return ann
            else:
                partial_ann = ann

                if hasattr(partial_ann, 'refmissing'):
                    if len(partial_ann.refmissing) < leastmissing:
                        leastmissing_feat = found[i]
                        leastmissing = len(partial_ann.refmissing)
                else:
                    leastmissing_feat = found[i]
                    leastmissing = 0

                if self.verbose and self.verbosity > 0:
                    self.logger.info(self.logname
                                     + " Using partial annotation * run "
                                     + str(i) + " *")
                    self.logger.info(self.logname
                                     + " Features found ("
                                     + str(len(ann.features.keys())) + ") = "
                                     + ",".join(list(ann
                                                     .features.keys())))
                    self.logger.info(self.logname
                                     + " Features missing ("
                                     + str(len(ann.missing.keys())) + ") = "
                                     + ",".join(list(ann
                                                     .missing.keys())))
                    self.logger.info(self.logname
                                     + " Sequence unmapped = "
                                     + str(ann.covered))
                    self.logger.info(self.logname + " ############" +
                                     "##################")

        # The number of sequences being used for alignment
        # can't be greater than the number of sequences
        # to be returned from the blast results
        if alignseqs > len(found):
            alignseqs = len(found)-1

        # * HARD CODED LOGIC * #
        # > After testing with multiple thresholds
        #   this value seemed to work best.
        #
        # Aligned % cutoff
        align_cutoff = .90
        if((not hasattr(partial_ann, 'features') or
           len(partial_ann.features) == 0)
           and len(sequence) > 700 and self.safemode):
            self.logger.error("No feature matches!")
            self.logger.error("Running in safe mode. " +
                              "No alignments will be done!")

            if rerun:
                # Check to see if reverse comp
                # TODO: Add note for reverse complement
                self.logger.info("Running with reverse complement.")
                sequence = sequence.reverse_complement()
                return self.annotate(sequence=sequence,
                                     locus=locus,
                                     rerun=False)
            return

        # Now loop through doing alignment
        for i in range(0, alignseqs):

            # Skip a reference
            # * For validation *
            if found[i].name in skip:
                continue

            if self.verbose:
                self.logger.info(self.logname
                                 + " running ref_align with "
                                 + found[i].name)

            aligned_ann = self.ref_align(found[i], sequence, locus,
                                         annotation=partial_ann,
                                         run=i,
                                         cutoff=align_cutoff)

            # * HARD CODED LOGIC * #
            # > If sequences are very novel, then the alignment
            #   cutoff may be to stringent. Incrementally decreasing
            #   the cutoff improves the likelihood of these sequences
            #   being annotated.
            align_cutoff -= .01
            if aligned_ann and aligned_ann.complete_annotation:
                if self.align:
                    if self.verbose:
                        self.logger.info(self.logname + " Adding alignment")
                    aligned_ann = self.add_alignment(found[i], aligned_ann)

                if self.verbose:
                    self.logger.info(self.logname
                                     + " Finished ref_align annotation using "
                                     + found[i].name)

                if self.verbose and self.verbosity > 0:
                    self.logger.info(self.logname
                                     + " Features annotated = "
                                     + ",".join(list(aligned_ann
                                                     .annotation.keys())))
                    if self.verbosity > 2:
                        for f in aligned_ann.annotation:
                            self.logger.info(self.logname
                                             + " " + f + " = "
                                             + str(aligned_ann
                                                   .annotation[f].seq))

                # Create GFE
                if create_gfe:
                    feats, gfe = self.gfe.get_gfe(aligned_ann, locus)
                    aligned_ann.gfe = gfe
                    aligned_ann.structure = feats
                aligned_ann.clean()
                return aligned_ann
            elif(aligned_ann):

                if self.verbose and self.verbosity > 0:
                    self.logger.info(self.logname
                                     + " Using partial annotation "
                                     + "for alignment * run "
                                     + str(i) + " - cutoff = "
                                     + str(align_cutoff)
                                     + " *")
                    self.logger.info(self.logname
                                     + " Features found = "
                                     + ",".join(list(aligned_ann
                                                     .features.keys())))
                    self.logger.info(self.logname
                                     + " Features missing = "
                                     + ",".join(list(aligned_ann
                                                     .missing.keys())))
                    self.logger.info(self.logname + " ############" +
                                     "##################")
                partial_ann = aligned_ann

                if(hasattr(partial_ann, 'annotation')
                   and partial_ann.annotation):
                    exon_only = True
                    for f in partial_ann.annotation:
                        if re.search("intron", f) or re.search("UTR", f):
                            exon_only = False

                    if(is_classII(locus) and exon_only
                       and len(partial_ann.annotation.keys()) > 0
                       and align_cutoff < .9):
                        align_cutoff = .80

            if not is_classII(locus) and align_cutoff < .88:
                align_cutoff = .88

        # Don't run full
        # annotation if flag is passed
        if not full:
            return

        if self.verbose:
            self.logger.info(self.logname + " running full alignment")

        # Try doing full alignment
        full_align = self.ref_align(leastmissing_feat,
                                    sequence,
                                    locus,
                                    partial_ann=partial_ann,
                                    cutoff=.80)

        if self.verbose:
            self.logger.info(self.logname
                             + " Finished ref_align annotation using full "
                             + leastmissing_feat.name)

        # Check to see if an annotation was returned
        if(not isinstance(full_align, Annotation)
           or isinstance(full_align, str)):
            if(not rerun or len(sequence) > 4000):
                self.logger.info(self.logname + " Failed annotation!")
                return Annotation()
            else:
                if self.verbose and self.verbosity > 0:
                    self.logger.info(self.logname
                                     + " Reruning annotation!")
                return self.annotate(sequence=sequence,
                                     locus=locus,
                                     alignseqs=2,
                                     nseqs=nseqs+1,
                                     skip=[found[0].name],
                                     rerun=False)

        # Check if the annotation is complete
        if not full_align.complete_annotation and self.verbose:
            self.logger.info(self.logname + " Incomplete annotation!")

        # Add the alignment to the annotation
        if self.align and full_align.complete_annotation:
            if self.verbose:
                self.logger.info(self.logname + " Adding alignment")
            full_align = self.add_alignment(leastmissing_feat, full_align)

        if self.verbose and self.verbosity > 0:
            self.logger.info(self.logname
                             + " Features annotated = "
                             + ",".join(list(full_align
                                             .annotation.keys())))
            if self.verbosity > 2:
                for f in full_align.annotation:
                    self.logger.info(self.logname
                                     + " " + f + " = "
                                     + str(full_align
                                           .annotation[f].seq))
        # Create GFE
        if create_gfe and full_align.complete_annotation:
            feats, gfe = self.gfe.get_gfe(full_align, locus)
            full_align.gfe = gfe
            full_align.structure = feats

        full_align.clean()
        if(full_align.complete_annotation
           or not rerun or len(sequence) > 4000):
            return full_align
        else:
            if self.verbose and self.verbosity > 0:
                self.logger.info(self.logname
                                 + " Reruning annotation!")
            return self.annotate(sequence=sequence,
                                 locus=locus,
                                 nseqs=nseqs+1,
                                 alignseqs=2,
                                 skip=[found[0].name],
                                 rerun=False)