class SequenceSearcher:

    def __init__(self, search_hmm, aln_hmm=None):
        self.search_hmm = search_hmm
        self.aln_hmm = aln_hmm

    def _get_sequence_directions(self, search_result):
        sequence_directions = {}
        for result in search_result:
            for hit in result.each([SequenceSearchResult.QUERY_ID_FIELD, SequenceSearchResult.ALIGNMENT_DIRECTION]):
                sequence_directions[hit[0]]  ={"strand":[hit[1]], "entry": []}
        return sequence_directions


    def _hmmalign(self, input_path, directions, pipeline,
                  forward_reads_output_path, reverse_reads_output_path):
        '''
        Align reads to the aln_hmm. Receives unaligned sequences and
        aligns them.

        Parameters
        ----------
        input_path : str
            Filename of unaligned hits to be aligned
        directions : dict
            dictionary containing read names as keys, and complement
            as the entry (True=Forward, False=Reverse)
        pipeline: str
            either PIPELINE_AA = "P" or PIPELINE_NT = "D"
        forward_reads_output_fh: str
            Where to write aligned forward reads
        reverse_reads_output_fh: str
            Where to write aligned reverse reads
        Returns
        -------
        Nothing.
        '''
        if pipeline == PIPELINE_AA:
            reverse_direction_reads_present=False
        else:
            reverse_direction_reads_present=False in directions.values()

        with tempfile.NamedTemporaryFile(prefix='for_file', suffix='.fa') as for_file_fh:
            for_file = for_file_fh.name
            with tempfile.NamedTemporaryFile(prefix='rev_file', suffix='.fa') as rev_file_fh:
                rev_file = rev_file_fh.name
                # Align input reads to a specified hmm.
                if reverse_direction_reads_present:  # Any that are in the reverse direction would be True
                    reverse = []
                    forward = []
                    records = list(SeqIO.parse(open(input_path), 'fasta'))

                    # Split the reads into reverse and forward lists
                    for record in records:
                        read_id = record.id

                        if directions[read_id] == True:
                            forward.append(record)
                        elif directions[read_id] == False:
                            reverse.append(record)
                        else:
                            raise Exception(logging.error('Programming error: hmmalign'))
                            exit(1)

                    logging.debug("Found %i forward direction reads" % len(forward))
                    logging.debug("Found %i reverse direction reads" % len(reverse))

                    # Write reverse complement and forward reads to files
                    with open(for_file, 'w') as for_aln:
                        logging.debug("Writing forward direction reads to %s" % for_file)
                        for record in forward:
                            for_aln.write('>' + record.id + '\n')
                            for_aln.write(str(record.seq) + '\n')
                            # HMMalign and convert to fasta format
                    if any(forward):
                        self.hmmalign_sequences(self.aln_hmm, for_file, forward_reads_output_path)
                    else:
                        cmd = 'touch %s' % (forward_reads_output_path)
                        extern.run(cmd)
                    with open(rev_file, 'w') as rev_aln:
                        logging.debug("Writing reverse direction reads to %s" % rev_file)
                        for record in reverse:
                            if record.id and record.seq:
                                rev_aln.write('>' + record.id + '\n')
                                rev_aln.write(str(record.seq.reverse_complement()) + '\n')
                    self.hmmalign_sequences(self.aln_hmm, rev_file, reverse_reads_output_path)
                    conv_files = [forward_reads_output_path, reverse_reads_output_path]
                    return conv_files

                else:
                    # If there are only forward reads, just hmmalign and be done with it.
                    self.hmmalign_sequences(self.aln_hmm, input_path, forward_reads_output_path)
                    conv_files = [forward_reads_output_path]

                    return conv_files

    def hmmalign_sequences(self, hmm, sequences, output_file):
        '''Run hmmalign and convert output to aligned fasta format

        Parameters
        ----------
        hmm: str
            path to hmm file
        sequences: str
            path to file of sequences to be aligned
        output_file: str
            write sequences to this file

        Returns
        -------
        nothing
        '''

        cmd = 'hmmalign --trim %s %s' % (hmm, sequences)
        output = extern.run(cmd)
        with open(output_file, 'w') as f:
            SeqIO.write(SeqIO.parse(StringIO(output), 'stockholm'), f, 'fasta')

    def makeSequenceBinary(self, sequences, fm):
        cmd = 'makehmmerdb %s %s' % (sequences, fm)
        extern.run(cmd)

    def hmmsearch(self, output_path, input_path, unpack, seq_type, threads, cutoff, orfm):
        '''
        hmmsearch - Search raw reads for hits using search_hmm list

        Parameters
        ----------
        output_path : str
            path to output domtblout table
        input_path : str
            path to input sequences to search
        unpack : UnpackRawReads
            UnpackRawReads object, returns string command that will output
            sequences to stdout when called on command line
            (use: unpack.command_line())
        seq_type : str
            variable containing a string, either 'nucleotide' or 'aminoacid'.
            Tells the pipeline whether or not to call ORFs on the sequence.
            If sequence is 'nucleotide', ORFs are called. If not, no ORFs.
        threads : Integer
            Number of threads to use. Passed to HMMsearch command.
        cutoff : str
            cutoff for HMMsearch to use, either an evalue or --cut_tc, meaning
            use the TC score cutoff specified within the HMM. Passed to
            HMMsearch command.
        orfm : OrfM
            Object that builds the command chunk for calling ORFs on sequences
            coming through as stdin. Outputs to stdout. Calls command_line
            to construct final command line string.

        Returns
        -------
        output_table_list : array of HMMSearchResult
            Includes the name of the output domtblout table given by hmmer

        Raises
        ------
        hmmsearcher.NoInputSequencesException
            Raised if there are no sequences fed into the HMM.
        '''

        # Define the base hmmsearch command.
        logging.debug("Using %i HMMs to search" % (len(self.search_hmm)))
        output_table_list = []
        if len(self.search_hmm) > 1:
            for hmm in self.search_hmm:
                out = os.path.join(os.path.split(output_path)[0], os.path.basename(hmm).split('.')[0] + '_' + os.path.split(output_path)[1])
                output_table_list.append(out)
        elif len(self.search_hmm) == 1:
            output_table_list.append(output_path)
        else:
            raise Exception("Programming error: expected 1 or more HMMs")

        # Choose an input to this base command based off the file format found.
        if seq_type == 'nucleotide':  # If the input is nucleotide sequence
            input_cmd = orfm.command_line(input_path)
        elif seq_type == 'aminoacid':  # If the input is amino acid sequence
            input_cmd = unpack.command_line()
        else:
            raise Exception('Programming Error: error guessing input sequence type')

        # Run the HMMsearches
        if cutoff == "--cut_tc":
            searcher = HmmSearcher(threads, cutoff)
        else:
            searcher = HmmSearcher(threads, '--domE %s' % cutoff)
        searcher.hmmsearch(input_cmd, self.search_hmm, output_table_list)

        hmmtables = [HMMSearchResult.import_from_hmmsearch_table(x) for x in output_table_list]
        return hmmtables

    def merge_forev_aln(self, forward_aln_list, reverse_aln_list, outputs):
        '''
        merge_forev_aln - Merges forward and reverse alignments for a given run

        Parameters
        ----------
        forward_aln_list : array
            List of the forward or interleaved alignments for each of the runs
            given to graftM. 
        reverse_aln_list : array
            List of the reverse alignments for each of the runs
            given to graftM. 
            **MUST** be the the same length as forward_aln_list or empty (if
            interleaved)
        outputs : array
            List of paths to output file to which the merged alignments from the
            aln_list will go into. Must be exactly half the size of the aln_list
            (i.e. one output file for every forward and reverse file provided)

        Returns
        -------
        Nothing - output files are known.
        '''

        orfm_regex = OrfM.regular_expression()
        def remove_orfm_end(id):
            """ remove the orfm suffix from a read id """
            try:
                return orfm_regex.match(id).group(1)
            except AttributeError:
                return id

        dir_regex = re.compile(r'/[12]$')
        def trim_direction(id):
            """ remove fwd/rev suffix if present """
            return dir_regex.sub('',id)

        def get_read_base(id):
            return trim_direction(remove_orfm_end(id))

        def get_read_bases(record_iter):
            """ remove orfm suffix from keys in dict records """
            return {get_read_base(r.id):r for r in record_iter}

        def split_interleaved_reads(records):
            """ return list of fwd records and dict of reverse """
            fwd = []
            rev = {}
            for record in records:
                id = get_read_base(record.id)
                if remove_orfm_end(record.id).endswith('1'):
                    fwd.append(record)
                else:
                    rev[id] = record
            return fwd, rev

        if not reverse_aln_list:
            # create dummy array for interleaved
            reverse_aln_list = [False,]*len(forward_aln_list)

        for idx, (forward_path, reverse_path) in enumerate(zip(forward_aln_list, reverse_aln_list)):
            output_path = outputs[idx]
            logging.info('Merging pair %s, %s' %
                         (os.path.basename(forward_path),
                          os.path.basename(reverse_path) \
                          if reverse_path else None))
            if reverse_path:
                forward_reads = SeqIO.parse(forward_path, 'fasta')
                reverse_reads = get_read_bases(SeqIO.parse(reverse_path, 'fasta'))
            else:
                forward_reads, reverse_reads = \
                        split_interleaved_reads(SeqIO.parse(forward_path,
                                                            'fasta'))
                #logging.debug("Loaded %d fwd (EG %s) and %d rev (EG %s)",
                #              len(forward_reads),
                #              get_read_base(next(iter(forward_reads)).id),
                #              len(reverse_reads),
                #              next(iter(reverse_reads)))

            with open(output_path, 'w') as out:
                for forward_record in forward_reads:
                    id = get_read_base(forward_record.id)
                    forward_sequence = str(forward_record.seq)
                    logging.debug("Fwd id %s is %s in reverse dict",
                                  id, "" if id in reverse_reads else "not")
                    try:
                        reverse_sequence = str(reverse_reads[id].seq)
                        new_seq = ''
                        if len(forward_sequence) == len(reverse_sequence):
                            for f, r in zip(forward_sequence, reverse_sequence):
                                if f == r:
                                    new_seq += f
                                elif f == '-' and r != '-':
                                    new_seq += r
                                elif r == '-' and f != '-':
                                    new_seq += f
                                elif f != '-' and r != '-':
                                    if f != r:
                                        new_seq += f
                                else:
                                    new_seq += '-'
                        else:
                            logging.error('Alignments do not match')
                            raise Exception('Merging alignments failed: Alignments do not match')
                        out.write('>%s\n' % forward_record.id)
                        out.write('%s\n' % (new_seq))
                        del reverse_reads[id]
                    except:

                        out.write('>%s\n' % forward_record.id)
                        out.write('%s\n' % (forward_sequence))
                for record_id, record in reverse_reads.items():
                    out.write('>%s\n' % record.id)
                    out.write('%s\n' % (str(record.seq)))

    def nhmmer(self, output_path, unpack, threads, evalue):
        '''
        nhmmer - Search input path using nhmmer

        Parameters
        ----------
        output_path : str
            A string containing the path to the input sequences
        unpack : obj
            UnpackRawReads object, returns string command that will output
            sequences to stdout when called on command line
            (use: unpack.command_line())
        threads : str
            Number of threads to run. For compiling command line.
        evalue : str
            evalue to use. For compiling commmand line.

        Returns
        -------
        output_table_list : array
            Includes the name of the output domtblout table given by hmmer
        '''
        logging.debug("Using %i HMMs to search" % (len(self.search_hmm)))
        output_table_list = []
        if len(self.search_hmm) > 1:
            for hmm in self.search_hmm:
                out = os.path.join(os.path.split(output_path)[0], os.path.basename(hmm).split('.')[0] + '_' + os.path.split(output_path)[1])
                output_table_list.append(out)
        elif len(self.search_hmm) == 1:
            output_table_list.append(output_path)
        else:
            raise Exception("Programming error: Expected 1 or more HMMs")
        input_pipe = unpack.command_line()

        searcher = NhmmerSearcher(threads, extra_args='--incE %s -E %s' % (evalue, evalue))
        searcher.hmmsearch(input_pipe, self.search_hmm, output_table_list)

        hmmtables = [HMMSearchResult.import_from_nhmmer_table(x) for x in output_table_list]

        return hmmtables, output_table_list

    def _check_euk_contamination(self, hmm_hit_tables):
        '''
        check_euk_contamination - Check output HMM tables hits reads that hit
                                  the 18S HMM with a higher bit score.

        Parameters
        ----------
        hmm_hit_tables : array
            Array of paths to the output files produced by hmmsearch or
            nhmmer.
        run_stats : dict
            A dictionary to updatewith the number of unique 18S reads and reads
            detected by both 18S and non-18S HMMs

        Returns
        -------
        euk_reads : set
            Non-redundant set of all read names deemed to be eukaryotic
        '''

        euk_hit_table = HMMreader(hmm_hit_tables.pop(-1))
        other_hit_tables = [HMMreader(x) for x in hmm_hit_tables]
        reads_unique_to_eukaryotes = []
        reads_with_better_euk_hit = []

        for hit in euk_hit_table.names():
            bits = []
            for hit_table in other_hit_tables:
                if hit in hit_table.names():
                    bits.append(hit_table.bit(hit))
                else:
                    reads_unique_to_eukaryotes.append(hit)
            if bits:
                if any([x for x in bits if x > euk_hit_table.bit(hit)]):
                    continue
                else:
                    reads_with_better_euk_hit.append(hit)
            else:
                continue

        if len(reads_with_better_euk_hit) == 0:
            logging.info("No contaminating eukaryotic reads detected")
        else:
            logging.info("Found %s read(s) that may be eukaryotic" % len(reads_with_better_euk_hit + reads_unique_to_eukaryotes))

        euk_reads = set(reads_with_better_euk_hit + reads_unique_to_eukaryotes)

        return euk_reads


    def _extract_multiple_hits(self, hits, reads_path, output_path):
        '''
        splits out regions of a read that hit the HMM. For example when two of
        same gene are identified within the same contig, The regions mapping to
        the HMM will be split out and written out to a new file as a new record.

        Parameters
        ----------
        hits : dict
            A dictionary where the keys are the read names, the entry for each
            is a list of lists, each containing the range within the contig
            (or read) that mapped to the HMM. e.g.:

            {'read1': [[3, 126], [305, 413]],
             'read2': [[1, 125]],
             ...

        reads_path : str
            path to reads file containing each read or contig in FASTA format to
            be opened, and split.
        output_path : str
            path to file to which split reads will be written to in FASTA
            format.

        Returns
        -------
        Nothing, output path is known.
        '''

        complement_information = {}
        try:
            reads = SeqIO.to_dict(SeqIO.parse(reads_path, "fasta"))  # open up reads as dictionary
        except:
            logging.error("Multiple sequences found with the same ID. The input sequences are either ill formated or are interleaved. \
If you provided GraftM with an interleaved sequence file, please split them into forward and reverse reads, and provide to the the appropriate \
flags (--forward, --reverse). Otherwise, it appears that you have provided sequences with redundant IDs. GraftM doesn't know how to \
deal with these, so please remove/rename sequences with duplicate keys.")
            raise InterleavedFileError()
        with open(output_path, 'w') as out:
            for read_name, entry in hits.items():  # For each contig
                ranges = entry["entry"]
                complements = entry["strand"]
                index = 1
                if len(ranges) > 1:  # if there are multiple hits in that contig
                    for r, c in zip(ranges, complements):  # for each of those hits

                        new_record = reads[read_name][r[0] - 1:r[1]]  # subset the record by the span of that hit
                        new_record.id = new_record.id + '_split_%i' % index  # give that subset record a new header
                        SeqIO.write(new_record, out, "fasta")  # and write it to output
                        index += 1  # increment the split counter
                        complement_information[new_record.id]=c

                else:  # Otherwise, just write the read back to the file
                    complement_information[read_name] = entry["strand"][0]
                    SeqIO.write(reads[read_name], out, "fasta")

        return complement_information

    def _extract_from_raw_reads(self, output_path, input_reads, raw_sequences_path, input_file_format, hits):
        '''
        _extract_from_raw_reads - Extract hit sequences of the hmm/diamond
        search from a command which generates uncompressed FASTA sequences from
        the raw sequences file. Output into specified file

        Parameters
        ----------
        output_path : str
            Path of the desired output file
        input_reads : list
            list of read IDS to extract.
        raw_sequences_path : str
            Path to the raw sequences
        input_file_format : var
            Variable, either FORMAT_FASTA_GZ, FORMAT_FASTQ_GZ, FORMAT_FASTQ or
            FORMAT_FASTA, denoting the format of the input sequence
        hits : dict
            A hash with the readnames as the keys and the spans as the values

        Returns
        -------
        output_path: str
            Path to file containing extracted reads.

        direction_information: dict

            {read_1: False
             ...
             read n: True}

            where True = Forward direction
            and False = Reverse direction
        '''

        with tempfile.NamedTemporaryFile(prefix='_raw_extracted_reads.fa') as tmp:
            # Extract reads from original sequence file
            extract_cmd = "mfqe --output-uncompressed"
            if input_file_format in (FORMAT_FASTA, FORMAT_FASTA_GZ, FORMAT_FASTQ, FORMAT_FASTQ_GZ):
                extract_cmd += " --fasta-read-name-lists /dev/stdin --input-fasta {} --output-fasta-files '{}'".format(
                    raw_sequences_path, tmp.name)
            else:
                raise Exception("Programming error: Unexpected input file format {}".format(input_file_format))

            extern.run(extract_cmd, stdin='\n'.join(input_reads))
            complement_info = self._extract_multiple_hits(hits, tmp.name, output_path)  # split them into multiple reads

        return output_path, complement_info


    def alignment_correcter(self, alignment_file_list, output_file_name,
                            filter_minimum=None):
        '''
        Remove lower case insertions in alignment outputs from HMM align. Give
        a list of alignments, and an output file name, and each alignment will
        be corrected, and written to a single file, ready to be placed together
        using pplacer.

        Parameters
        ----------
        alignment_file_list : array
            List of strings, each the path to different alignments from the
            inputs provided to GraftM
        output_file_name : str
            The path and filename of the output file desired.
        filter_minimum : int
            minimum number of positions that must be aligned for each sequence
        Returns
        -------
        True or False, depending if reads were written to file

        '''

        corrected_sequences = {}
        for alignment_file in alignment_file_list:
            insert_list = []  # Define list containing inserted positions to be removed (lower case characters)
            with open(alignment_file) as f:
                sequence_list = list(SeqIO.parse(f, 'fasta'))
            for sequence in sequence_list:  # For each sequence in the alignment
                for idx, nt in enumerate(list(sequence.seq)):  # For each nucleotide in the sequence
                    if nt.islower():  # Check for lower case character
                        insert_list.append(idx)  # Add to the insert list if it is
            insert_list = list(OrderedDict.fromkeys(sorted(insert_list, reverse=True)))  # Reverse the list and remove duplicate positions
            for sequence in sequence_list:  # For each sequence in the alignment
                new_seq = list(sequence.seq)  # Define a list of sequences to be iterable list for writing
                for position in insert_list:  # For each position in the removal list
                    del new_seq[position]  # Delete that inserted position in every sequence
                corrected_sequences['>' + sequence.id + '\n'] = (''.join(new_seq) + '\n').replace('~', '-')

        pre_filter_count=len(corrected_sequences)

        if filter_minimum:
            # Use '>' not '>=' here because the sequence is on a single line,
            # but also includes a newline character at the end of the sequence
            corrected_sequences={key:item for key, item in iter(corrected_sequences.items()) if len(item.replace('-', '')) > filter_minimum}

        post_filter_count=len(corrected_sequences)
        logging.info("Filtered %i short sequences from the alignment" % \
                        (pre_filter_count-post_filter_count)
                    )
        logging.info("%i sequences remaining" % post_filter_count)

        if len(corrected_sequences) >= 1:
            with open(output_file_name, 'w') as output_file:  # Create an open file to write the new sequences to
                for fasta_id, fasta_seq in corrected_sequences.items():
                    output_file.write(fasta_id)
                    output_file.write(fasta_seq)
            return True
        else:
            return False


    def _extract_orfs(self, input_path, orfm, hit_readnames, output_path, search_method, sequence_frame_info_list=None):
        '''
        Call ORFs on a file with nucleotide sequences and extract the proteins
        whose name is in `hit_readnames`.

        Parameters
        ----------
        input_path : str
            Path to input nucleotide sequences in FASTA format.
        orfm: obj
            graftm.OrfM object with parameters already set
        hit_readnames : str
            path to a file containin the readnames of hits to the HMM, one per
            line.
        output_path : str
            Path to output orfs into, in FASTA format.
        search_method : str
            The method for searching, either 'hmmsearch' or 'diamond'
        sequence_frame_info : list
            A dataframe (list of lists) containing readname, alignment direction
            and alignment start point information
        '''

        if search_method == "hmmsearch":
            # Build and run command to extract ORF sequences:
            orfm_cmd = orfm.command_line()
            cmd = "mfqe --output-uncompressed --fasta-read-name-lists /dev/stdin --input-fasta <({} {}) --output-fasta-files {}".format(
                orfm_cmd, input_path, output_path)
            extern.run(cmd, stdin='\n'.join(hit_readnames))

        elif search_method == "diamond":
            sequence_frame_info_dict = {x[0]:[x[1], x[2], x[3]] for x in sequence_frame_info_list}
            records = SeqIO.parse(input_path, "fasta")
            with open(output_path, 'w') as open_output_path:
                for record in records:
                    entry=sequence_frame_info_dict[record.id]
                    indfrom=(min(entry[2], entry[1])-1)
                    indto=max(entry[2], entry[1])
                    if entry[0] == False:
                        record.seq = max([x for x in record.seq[indfrom:indto].reverse_complement().translate().split("*")], key=len)
                    else:
                        record.seq = max([x for x in record.seq[indfrom:indto].translate().split("*")], key=len)
                    SeqIO.write(record, open_output_path, "fasta")
                open_output_path.flush()


    def _get_read_names(self, search_result, max_range):
        '''
        _get_read_names - loops through hmm hits and their alignment spans to
        determine if they are potentially linked (for example, if one gene in a
        contig hits a hmm more than once, in two different conserved regions
        of that gene) and combines them into one 'hit' because they are
        technically the same. The total span of the hits deemed to be linked is
        returned.

        Parameters
        ----------
        search_result : obj
            SequenceSearchResult object with all paramaters defined. Used here
            to create rows containing information on alignment direction and
            alignment span.
        max_range : int
            Maximum range that a gene can extend within a contig. Any hits
            that extend beyond this length cannot be linked. max_range is
            set as 1.5 X the average length of all full length genes used
            in the search database. This is defined in the CONTENTS.json file
            within a gpkg.
        Returns
        -------
            Dictionary where keys are the contig/read name. The value for each
            entry is an array lists, one per hit in each contig, each with the
            span (min and max) of the alignment.
        '''

        splits = {}  # Define an output dictionary to be filled
        spans = []
        for result in search_result:  # Create a table (list of rows contain span, and complement information
            spans += list(
                         result.each(
                                    [SequenceSearchResult.QUERY_ID_FIELD,
                                        SequenceSearchResult.ALIGNMENT_DIRECTION,
                                        SequenceSearchResult.HIT_FROM_FIELD,
                                        SequenceSearchResult.HIT_TO_FIELD,
                                        SequenceSearchResult.QUERY_FROM_FIELD,
                                        SequenceSearchResult.QUERY_TO_FIELD]
                                       )
                         )

        for hit in spans:  # For each of these rows (i.e. hits)
            i = hit[0]  # set id to i
            c = hit[1]  # set complement to c
            ft = [min(hit[2:4]), max(hit[2:4])]  # set span as ft (i.e. from - to) - This is the amount covering the query
            qs = [min(hit[4:6]), max(hit[4:6])]  # seq the query span to qs - This is the amount covering the HMM

            if ft[0] == ft[1]: continue  # if the span covers none of the query, skip that entry (seen this before)

            if i not in splits:  # If the hit hasnt been seen yet
                splits[i] = {'span'       : [ft],
                             'strand'     : [c],
                             'query_span' : [qs]}  # add the span and complement as new entry

            else:  # otherwise (if it has been seen)
                for idx, entry in enumerate(splits[i]['span']):  # for each previously added entry
                    if splits[i]['strand'][idx] != c:  # If the hit is on the same complement strand
                        splits[i]['span'].append(ft)  # Add the new range to be split out in the future
                        splits[i]['strand'].append(c)  # Add the complement strand as well
                        splits[i]['query_span'].append(qs)
                        break

                    previous_qs      = splits[i]['query_span'][idx] # Get the query span of the previous hit

                    previous_q_range = set(range(previous_qs[0], previous_qs[1])) # Get the range of each
                    current_q_range  = set(range(qs[0], qs[1]))
                    query_overlap    = set(previous_q_range).intersection(current_q_range) # Find the intersection between the two ranges

                    previous_ft_span = set(range(entry[0], entry[1]))
                    current_ft_span  = set(range(ft[0], ft[1]))
                    if any(query_overlap): # If there is an overlap
                        ####################################################
                        # if the span over the actual read that hit the HMM
                        # for each hit overlap by > 25%, they are considered
                        # the same hit, and ignored
                        ####################################################
                        intersection_fraction = float(len(previous_ft_span.intersection(current_ft_span)))

                        if intersection_fraction / float(len(previous_ft_span)) >= PREVIOUS_SPAN_CUTOFF:
                            break
                        elif intersection_fraction / float(len(current_ft_span)) >= PREVIOUS_SPAN_CUTOFF:
                            break
                        else: # else (i.e. if the hit covers less that 25% of the sequence of the previous hit)
                            ####################################################
                            # If they made it this far, it means that the hits do not overlap.
                            # But one last check must be made to ensure they do not cover the same
                            # region in the HMM.
                            ####################################################
                            if len(query_overlap) > (len(current_q_range)*PREVIOUS_SPAN_CUTOFF): # if the overlap on the query HMM does not span over 25%
                                if (idx+1) == len(splits[i]['span']):
                                    splits[i]['span'].append(ft) # Add from-to as another entry, this is another hit.
                                    splits[i]['strand'].append(c) # Add strand info as well
                                    splits[i]['query_span'].append(qs)
                                    break # And break

                    if min(entry) < min(ft):  # if/else to determine which entry comes first (e.g. 1-5, 6-10 not 6-10, 1-5)
                        if max(ft) - min(entry) < max_range:  # Check if they lie within range of eachother
                            entry[1] = max(ft)  # ammend the entry if they are
                            break  # And break the loop
                    else:
                        if max(entry) - min(ft) < max_range:  # Check if they lie within range of eachother
                            entry[0] = min(ft)  # ammend the entry if they are
                            break  # And break the loop

                else:  # if no break occured (no overlap)
                    splits[i]['span'].append(ft)  # Add the new range to be split out in the future
                    splits[i]['strand'].append(c)  # Add the complement strand as well
                    splits[i]['query_span'].append(qs)
        return {key: {"entry":entry['span'], 'strand': entry['strand']} for key, entry in iter(splits.items())}  # return the dict, without strand information which isn't required.

    def _check_for_slash_endings(self, readnames):
        '''
        Provide a list of read names to be checked for the /1 or /2 endings
        which cause troubles downstream when working with paired reads. This
        function checks the list for ANY read ending in /1 or /2. The reads are
        first stripped of any comment lines, using the 'split' function

        Parameters
        ----------
        readnames : list
            list of strings, each a read name to be checked for the /1 or /2
            ending

        Returns
        -------
        Boolean. True if /1 or /2 endings are found, False if not.

        '''

        # Strip comment lines, if any
        readnames = [x.split()[0] for x in readnames]

        # Check for /1 and /2 endings:
        if any(x for x in readnames if x.endswith('/1') or x.endswith('/2')):
            result = True
        else:
            result = False

        return result

    @T.timeit
    def aa_db_search(self, files, base, unpack, search_method,
                     maximum_range, threads, evalue, min_orf_length,
                     restrict_read_length, diamond_database,
                     diamond_performance_parameters):
        '''
        Amino acid database search pipeline - pipeline where reads are searched
        as amino acids, and hits are identified using hmmsearch or diamond
        searches

        Parameters
        ----------
        files : obj
            graftm_output_paths object.
        base : str
            The name of the input file, stripped of all suffixes, and paths.
            Used for creating file names with 'files' object.
        unpack : obj
            UnpackRawReads object, returns string command that will output
            sequences to stdout when called on command line
            (use: unpack.command_line())
        search_method : str
            The method for searching, either 'hmmsearch' or 'diamond'
        maximum_range : int
            Maximum range that a gene can extend within a contig. Any hits
            that extend beyond this length cannot be linked. max_range is defined
            as 1.5 X the average length of all full length genes used in the
            search database. This is defined in the CONTENTS.json file within a
            gpkg.
        threads : int
            Number of threads for hmmer to use
        evalue : str
            evalue cutoff for hmmer to use
        min_orf_length : int
            minimum orf length for orfm to use
        restrict_read_length : int
            orf length to retrict orfm to.
        diamond_database : str
            Path to diamond database to use when searching. Set to 'None' if not
            using diamond pipeline
        diamond_performance_parameters : str
            Extra arguments to provide to DIAMOND
        Returns
        -------
        String path to amino acid fasta file of reads that hit
        '''
        # Define outputs
        if search_method == 'hmmsearch':
            output_search_file = files.hmmsearch_output_path(base)
        elif search_method == 'diamond':
            output_search_file = files.diamond_search_output_basename(base)
        hit_reads_fasta = files.fa_output_path(base)
        hit_reads_orfs_fasta = files.orf_fasta_output_path(base)

        return self.search_and_extract_orfs_matching_protein_database(\
                                                    unpack,
                                                    search_method,
                                                    maximum_range,
                                                    threads,
                                                    evalue,
                                                    min_orf_length,
                                                    restrict_read_length,
                                                    diamond_database,
                                                    output_search_file,
                                                    hit_reads_fasta,
                                                    hit_reads_orfs_fasta,
                                                    diamond_performance_parameters)

    def search_and_extract_orfs_matching_protein_database(self,
                                                      unpack,
                                                      search_method,
                                                      maximum_range,
                                                      threads,
                                                      evalue,
                                                      min_orf_length,
                                                      restrict_read_length,
                                                      diamond_database,
                                                      output_search_file,
                                                      hit_reads_fasta,
                                                      hit_reads_orfs_fasta,
                                                      diamond_performance_parameters):
        '''As per aa_db_search() except slightly lower level. Search an
        input read set (unpack) and then extract the proteins that hit together
        with their containing nucleotide sequences.

        Parameters
        ----------
        output_search_file: str
            path to hmmsearch output table or diamond basename
        hit_reads_fasta: str
            path to nucleotide sequences containing hit proteins
        hit_reads_orfs_fasta: str
            path to hit proteins, unaligned

        Returns
        -------
        direction_information: dict

            {read_1: False
             ...
             read n: True}

            where True = Forward direction
            and False = Reverse direction

        result: DBSearchResult object containing file locations and hit
        information

        '''

        orfm = OrfM(min_orf_length=min_orf_length,
                     restrict_read_length=restrict_read_length)
        extracting_orfm = OrfM(min_orf_length=min_orf_length,
                      restrict_read_length=restrict_read_length)

        if search_method == 'hmmsearch':
            # run hmmsearch
            search_result = self.hmmsearch(
                                           output_search_file,
                                           unpack.get_file_as_process(),
                                           unpack,
                                           unpack.sequence_type(),
                                           threads,
                                           evalue,
                                           orfm
                                           )

        elif search_method == 'diamond':
            # run diamond
            search_result = Diamond(
                                     database=diamond_database,
                                     threads=threads,
                                     evalue=evalue,
                                     ).run(
                                           unpack.get_file_as_process(),
                                           unpack.sequence_type(),
                                           daa_file_basename=output_search_file,
                                           extra_args=diamond_performance_parameters,
                                           )
            search_result = [search_result]

        else:  # if the search_method isn't recognised
            raise Exception("Programming error: unexpected search_method %s" % search_method)

        orfm_regex = OrfM.regular_expression()

        hits = self._get_sequence_directions(search_result)

        orf_hit_readnames = list(hits.keys()) # Orf read hit names
        if unpack.sequence_type() == 'nucleotide':
            hits={(orfm_regex.match(key).groups(0)[0] if orfm_regex.match(key) else key): item for key, item in iter(hits.items())}
            hit_readnames = list(hits.keys()) # Store raw read hit names
        else:
            hit_readnames=orf_hit_readnames

        hit_reads_fasta, direction_information = self._extract_from_raw_reads(
                                                       hit_reads_fasta,
                                                       hit_readnames,
                                                       unpack.get_file_as_process(),
                                                       unpack.format(),
                                                       hits
                                                       )


        if not hit_readnames:
            hit_read_counts = [0, len(hit_readnames)]
            result = DBSearchResult(None,
                                    search_result,
                                    hit_read_counts,
                                    None)

            return result, direction_information


        if unpack.sequence_type() == 'nucleotide':
            # Extract the orfs of these reads that hit the original search
            self._extract_orfs(
                               hit_reads_fasta,
                               extracting_orfm,
                               orf_hit_readnames,
                               hit_reads_orfs_fasta,
                               search_method,
                               list(search_result[0].each([SequenceSearchResult.QUERY_ID_FIELD,
                                                           SequenceSearchResult.ALIGNMENT_DIRECTION,
                                                           SequenceSearchResult.QUERY_FROM_FIELD,
                                                           SequenceSearchResult.QUERY_TO_FIELD])
                                    )
                               )

            hit_reads_fasta = hit_reads_orfs_fasta
        slash_endings=self._check_for_slash_endings(hit_readnames)
        result = DBSearchResult(hit_reads_fasta,
                                search_result,
                                [0, len([itertools.chain(*list(hits.values()))])],  # array of hits [euk hits, true hits]. Euk hits alway 0 unless searching from 16S
                                slash_endings)  # Any reads that end in /1 or /2

        if maximum_range:
            n_hits = sum([len(x["strand"]) for x in hits.values()])
        else:
            n_hits = len(hits.keys())
        logging.info("%s read(s) detected" % n_hits)

        return result, direction_information


    @T.timeit
    def nt_db_search(self, files, base, unpack, euk_check,
                     search_method, maximum_range, threads, evalue):
        '''
        Nucleotide database search pipeline - pipeline where reads are searched
        as nucleotides, and hits are identified using nhmmer searches

        Parameters
        ----------
        files : obj
            graftm_output_paths object.
        base : str
            The name of the input file, stripped of all suffixes, and paths.
            Used for creating file names with 'files' object.
        unpack : obj
            UnpackRawReads object, returns string command that will output
            sequences to stdout when called on command line
            (use: unpack.command_line())
        euk_check : bool
            True indicates the sample will be checked for eukaryotic reads,
            False indicates not.
        search_method : str
            The method for searching e.g. 'hmmsearch' or 'diamond'
        maximum_range : int
            Maximum range that a gene can extend within a contig. Any hits
            that extend beyond this length cannot be linked. max_range is defined
            as 1.5 X the average length of all full length genes used in the
            search database. This is defined in the CONTENTS.json file within a
            gpkg.
        threads : str
            Number of threads for hmmer to use
        evalue : str
            Evalue cutoff for hmmer to use
        Returns
        -------
        String path to amino acid fasta file of reads that hit
        '''

        # Define outputs
        hmmsearch_output_table = files.hmmsearch_output_path(base)
        hit_reads_fasta = files.fa_output_path(base)
        return \
            self.search_and_extract_nucleotides_matching_nucleotide_database(\
                                                      unpack,
                                                      euk_check,
                                                      search_method,
                                                      maximum_range,
                                                      threads,
                                                      evalue,
                                                      hmmsearch_output_table,
                                                      hit_reads_fasta)

    def search_and_extract_nucleotides_matching_nucleotide_database(self,
                                                      unpack,
                                                      euk_check,
                                                      search_method,
                                                      maximum_range,
                                                      threads,
                                                      evalue,
                                                      hmmsearch_output_table,
                                                      hit_reads_fasta):
        '''As per nt_db_search() except slightly lower level. Search an
        input read set (unpack) and then extract the sequences that hit.

        Parameters
        ----------
        hmmsearch_output_table: str
            path to hmmsearch output table
        hit_reads_fasta: str
            path to hit nucleotide sequences

        Returns
        -------
        direction_information: dict

            {read_1: False
             ...
             read n: True}

            where True = Forward direction
            and False = Reverse direction

        result: DBSearchResult object containing file locations and hit
        information
        '''

        if search_method == "hmmsearch":
            # First search the reads using the HMM
            search_result, table_list = self.nhmmer(
                                                    hmmsearch_output_table,
                                                    unpack,
                                                    threads,
                                                    evalue
                                                    )


        elif search_method == 'diamond':
            raise Exception("Diamond searches not supported for nucelotide databases yet")


        if maximum_range:

            hits = self._get_read_names(
                                        search_result,  # define the span of hits
                                        maximum_range
                                        )
        else:
            hits = self._get_sequence_directions(search_result)

        hit_readnames = list(hits.keys())

        if euk_check:
            euk_reads = self._check_euk_contamination(table_list)
            hit_readnames = set([read for read in hit_readnames if read not in euk_reads])
            hits = {key:item for key, item in  iter(hits.items()) if key in hit_readnames}
            hit_read_count = [len(euk_reads), len(hit_readnames)]
        else:
            hit_read_count = [0, len(hit_readnames)]

        hit_reads_fasta, direction_information = self._extract_from_raw_reads(
                                                       hit_reads_fasta,
                                                       hit_readnames,
                                                       unpack.get_file_as_process(),
                                                       unpack.format(),
                                                       hits
                                                       )

        if not hit_readnames:
            result = DBSearchResult(None,
                                  search_result,
                                  hit_read_count,
                                  None)
        else:
            slash_endings=self._check_for_slash_endings(hit_readnames)
            result = DBSearchResult(hit_reads_fasta,
                                    search_result,
                                    hit_read_count,
                                    slash_endings)

        if maximum_range:
            n_hits = sum([len(x["strand"]) for x in hits.values()])
        else:
            n_hits = len(hits)
        logging.info("%s read(s) detected" % n_hits)

        return result, direction_information

    @T.timeit
    def align(self, input_path, output_path, directions, pipeline,
              filter_minimum):
        '''align - Takes input path to fasta of unaligned reads, aligns them to
        a HMM, and returns the aligned reads in the output path

        Parameters
        ----------
        input_path : str
        output_path : str
        reverse_direction : dict
            A dictionary of read names, with the entries being the complement
            strand of the read (True = forward, False = reverse)
        pipeline : str
            Either "P" or "D" corresponding to the protein and nucleotide (DNA)
            pipelines, respectively.


        Returns
        -------
        N/A - output alignment path known.
        '''

        # HMMalign the forward reads, and reverse complement reads.
        with tempfile.NamedTemporaryFile(prefix='for_conv_file', suffix='.fa') as fwd_fh:
            fwd_conv_file = fwd_fh.name
            with tempfile.NamedTemporaryFile(prefix='rev_conv_file', suffix='.fa') as rev_fh:
                rev_conv_file = rev_fh.name
                alignments = self._hmmalign(
                    input_path,
                    directions,
                    pipeline,
                    fwd_conv_file,
                    rev_conv_file)
                alignment_result = self.alignment_correcter(alignments,
                                                            output_path,
                                                            filter_minimum)
                return alignment_result class DendropyTreeCleaner:
    def write_fasttree_newick(self, tree, output_file_io):
        '''Write the given dendropy tree to the open output file I/O, ensuring that is
        acceptable for FastTree by e.g. removing quoted taxon names.

        Parameters
        ----------
        tree: dendropy.Tree
            tree to be written
        output_file_io: stream
            stream to write the newick formatted tree to

        '''
        newick = tree.as_string(schema='newick',
                                suppress_internal_node_labels=True,
                                suppress_internal_taxon_labels=True,
                                preserve_spaces=False) # Labels must not be quoted.
        output_file_io.write(newick.replace("'",'').replace(' ','_').encode())

    def match_alignment_and_tree_sequence_ids(self, sequence_names, tree):
        '''Check to make sure that the sequences specified in the alignment
        and the tree are the same, otherwise raise an Exception detailing
        the problem for the user to fix

        Parameters
        ----------
        sequence_names: list of str
            names of sequences to ensure are in the tree
        tree: dendropy.Tree
            tree to find names in

        '''

        tip_names_count = {}
        for t in tree.leaf_node_iter():
            # replace spaces with underscores as this is how they are given to FastTree.
            name = t.taxon.label.replace(' ','_')
            if name in tip_names_count:
                raise Exception("Duplicate tip name found in tree: '%s'" % name)
            else:
                tip_names_count[name] = 1
        for name in sequence_names:
            if name not in tip_names_count:
                raise Exception("The alignment sequence '%s' was found in the alignment but not the tree" % name)
            elif tip_names_count[name] > 1:
                raise Exception("Found duplicate sequence name '%s'" % name)
            else:
                tip_names_count[name] += 1
        for name, count in tip_names_count.items():
            if count < 2:
                raise Exception("Sequence '%s' was found in the tree but not the alignment" % name)

    def remove_sequences(self, tree, sequence_names):
        '''Remove sequences with in the given sequence_names array from the tree in
        place. Assumes the sequences are found in the tree, and that they are
        all unique.

        Parameters
        ----------
        tree: dendropy.Tree
            tree to remove from
        sequence_names: list of str
            list of tip names to remove

        '''
        tree.prune_taxa_with_labels(sequence_names)
        tree.prune_taxa_with_labels([s.replace('_',' ') for s in sequence_names]) class TaxonomyExtractor:
    def taxonomy_from_annotated_tree(self, dendropy_tree):
        '''Given an annotated dendropy tree, return the taxonomy of each tip
        
        Parameters
        ---------
        tree: dendropy.Tree
            Decorated tree to extract taxonomy from       
        Returns
        -------
        dictionary of tip name to array of taxonomic affiliations'''
        tip_to_taxonomy = {}
        for tip in dendropy_tree.leaf_node_iter():
            tax = []
            n = tip.parent_node
            while n:
                node_taxonomy = self.taxonomy_from_node_name(n.label)
                if node_taxonomy and n.parent_node:
                    tax = [node_taxonomy]+tax
                n = n.parent_node
            tip_name = tip.taxon.label.replace(' ','_')
            if tip_name in tip_to_taxonomy:
                raise Exception("Found duplicate tip name '%s'" % tip_name)
            if len(tax)==0:
                tip_to_taxonomy[tip_name] = []
            else:
                tip_to_taxonomy[tip_name] = [t.strip() for t in '; '.join(tax).split(';')]
        return tip_to_taxonomy
    
    def taxonomy_from_node_name(self, node_name):
        '''return the taxonomy incorporated at a particular node, or None
        if it does not encode any taxonomy
        
        Parameters
        ----------
        node_name: str
            a node label.
            
        Returns
        -------
        Taxonomy as a string, or None if there is no taxonomy
        '''
        
        def is_float(s):
            try:
                float(s)
                return True
            except ValueError:
                return False
        
        if node_name is None:
            return None
        elif is_float(node_name):
            # no name, just a bootstrap
            return None
        else:
            bootstrap_regex = re.compile(r'[\d\.]+:(.*)')
            reg = bootstrap_regex.match(node_name)
            if reg:
                # bootstrap in name
                return reg.groups(0)[0]
            else:
                # bootstrap not in name
                return node_name class Getaxnseq:

    def _taxonomy_line(self, level_index, taxon_array, taxonomic_level_names):
        if level_index == 0:
            return '%s,Root,%s,%s,%s,%s%s' \
                % (taxon_array[level_index],
                    #parent
                   taxonomic_level_names[level_index],
                   taxon_array[level_index],
                   'Root',
                   taxon_array[level_index],
                   ''.join([','] * (len(taxonomic_level_names) - (level_index+1)) ))
        else:
            return '%s,%s,%s,%s,%s,%s%s' \
                % (taxon_array[level_index],
                   taxon_array[level_index-1],
                   taxonomic_level_names[level_index],
                   taxon_array[level_index],
                   'Root',
                   ','.join(taxon_array),
                   ''.join([','] * (len(taxonomic_level_names) - (level_index+1)) ))

    def read_taxtastic_taxonomy_and_seqinfo(self, taxonomy_io, seqinfo_io):
        '''Read the taxonomy and seqinfo files into a dictionary of
        sequence_name => taxonomy, where the taxonomy is an array of lineages
        given to that sequence.

        Possibly this method is unable to handle the full definition of these
        files? It doesn't return what each of the ranks are, for starters.
        Doesn't deal with duplicate taxon names either.
        '''

        # read in taxonomy file
        lineages = [] #array of hashes where each taxon points to its parent taxon's name
        taxon_to_lineage_index = {}
        expected_number_of_fields = None
        for line in taxonomy_io:
            splits = line.strip().split(',')
            if expected_number_of_fields is None:
                expected_number_of_fields = len(splits)
                lineages = [{}]* (expected_number_of_fields-4)
                continue #this is the header line
            elif len(splits) != expected_number_of_fields:
                raise Exception("Encountered error parsing taxonomy file, expected %i fields but found %i on line: %s" %
                                (expected_number_of_fields, len(splits), line))
            # e.g. 'tax_id,parent_id,rank,tax_name,root,kingdom,phylum,class,order,family,genus,species
            tax_id = splits[0]
            parent_id = splits[1]

            try:
                lineage_index = splits.index('')-5
            except ValueError:
                lineage_index = len(splits)-5
            taxon_to_lineage_index[tax_id] = lineage_index
            lineages[lineage_index][tax_id] = parent_id

        taxonomy_dictionary = {}
        for i, line in enumerate(seqinfo_io):
            if i==0: continue #skip header line

            splits = line.strip().split(',')
            if len(splits) != 2:
                raise Exception("Bad formatting of seqinfo file on this line: %s" % line)

            seq_name = splits[0]
            taxon = splits[1]
            lineage_index = taxon_to_lineage_index[taxon]
            if lineage_index==0:
                # Don't include Root in the taxonomy
                taxonomy_dictionary[seq_name] = []
            else:
                full_taxonomy_rev = []
                while lineage_index > 0:
                    full_taxonomy_rev.append(taxon)
                    taxon = lineages[lineage_index][taxon]
                    lineage_index = lineage_index-1
                taxonomy_dictionary[seq_name] = list(reversed(full_taxonomy_rev))

        return taxonomy_dictionary

    def write_taxonomy_and_seqinfo_files(self, taxonomies, output_taxonomy_file, output_seqinfo_file):
        '''Write out taxonomy and seqinfo files as required by taxtastic
        from known taxonomies

        Parameters
        ----------
        taxonomies:
            hash of taxon_id to array of taxonomic information
        output_taxonomy_file:
            write taxtastic-compatible 'taxonomy' file here
        output_seqinfo_file:
            write taxtastic-compatible 'seqinfo' file here'''

        first_pass_id_and_taxonomies = []
        tc=TaxonomyCleaner()
        max_number_of_ranks = 0

        for taxon_id, tax_split in taxonomies.items():
            # Replace spaces with underscores e.g. 'Candidatus my_genus'
            for idx, item in enumerate(tax_split):
                tax_split[idx] = re.sub('\s+', '_', item.strip())

            # Remove 'empty' taxononomies e.g. 's__'
            tax_split = tc.remove_empty_ranks(tax_split)

            # Add this fixed up list to the list
            first_pass_id_and_taxonomies.append([taxon_id]+tax_split)

            if len(tax_split) > max_number_of_ranks:
                max_number_of_ranks = len(tax_split)

        # Find taxons that have multiple parents, building a hash of parents as we go (i.e. a tree of taxonomies embedded in a hash)
        #
        # Assumes that no two taxonomic labels are the same when they are from different
        # taxonomic levels. When there are children with multiple parents at the
        # same taxonomic label then these are warned about and worked around.
        parents = {} #hash of taxon to immediate parent
        known_duplicates = set()
        for j, array in enumerate(first_pass_id_and_taxonomies):
            taxonomy = array[1:]
            for i, tax in enumerate(taxonomy):
                if i==0: continue #top levels don't have parents
                ancestry = taxonomy[i-1]
                if tax in parents:
                    if parents[tax] != ancestry:

                        dup = "%s%s" %(parents[tax], ancestry)
                        # don't report the same problem several times
                        if dup not in known_duplicates:
                            print(" %s '%s' with multiple parents %s and %s" % (array[0], tax, parents[tax], ancestry))
                            known_duplicates.add(dup)

                        # fix the current one
                        new_name_id = 1
                        new_name = "%s_graftm_%s" % (tax, new_name_id)

                        while new_name in parents and parents[new_name] != ancestry:
                            new_name_id += 1
                            new_name = "%s_graftm_%s" % (tax, new_name_id)

                        first_pass_id_and_taxonomies[j][i+1] = new_name
                        taxonomy[i] = new_name
                        parents[new_name] = ancestry
                else:
                    # normal case, seeing a new taxon and parent for the first time

                    parents[tax] = ancestry

        # Write the sequence file
        with open(output_seqinfo_file, 'w') as seqout:
            # write header
            seqout.write('seqname,tax_id\n')
            # write each taxonomic association
            for array in first_pass_id_and_taxonomies:
                if len(array)==1:
                    most_specific_taxonomic_affiliation = 'Root'
                else:
                    most_specific_taxonomic_affiliation = array[-1]

                seqout.write("%s,%s\n" % (array[0],
                                          most_specific_taxonomic_affiliation))

        # Write the taxonomy file
        noted_taxonomies = set()
        taxonomic_level_names = ["rank_%i" % rank for rank in range(max_number_of_ranks)]
        with open(output_taxonomy_file, 'w') as seqout:
            # write header and root line
            seqout.write(','.join(['tax_id','parent_id','rank','tax_name','root'] + taxonomic_level_names) +'\n')
            seqout.write(','.join(['Root','Root','root','Root','Root'])+ ''.join([',']*max_number_of_ranks) + '\n')
            # write all the taxonomies
            for array in first_pass_id_and_taxonomies:
                taxons = array[1:]
                for i, tax in enumerate(taxons):
                    line = self._taxonomy_line(i, taxons[:(i+1)], taxonomic_level_names)
                    if line not in noted_taxonomies:
                        seqout.write(line+"\n")
                        noted_taxonomies.add(line) class Deduplicator:
    r"""Deduplicates sequences"""
    
    def deduplicate(self, aligned_sequence_objects):
        r"""Sort the given aligned_sequence objects into an array of arrays,
        where input sequences are grouped iff they have the same sequence
        
        Parameters
        ----------
        aligned_sequence_objects: array of Sequence objects
            input sequences
            
        Returns
        -------
        Array of arrays of Sequence objects"""
        
        sequence_to_groups = {}
        for s in aligned_sequence_objects:
            try:
                sequence_to_groups[s.seq].append(s)
            except KeyError:
                sequence_to_groups[s.seq] = [s]
        return list(sequence_to_groups.values())
    
    def lca_taxonomy(self, deduplicated_sequences, taxonomy_hash):
        r'''Given a set of deduplicated sequences and a taxonomy hash,
        return the respective LCAs of taxonomy
        
        Parameters
        ----------
        deduplicated_sequences: Array of arrays of Sequence objects
            as output from deduplicate()
        taxonomy_hash: dictionary 
            of sequence names to taxonomy array (i.e. array of str)
        
        Returns
        -------
        Array of taxonomy LCAs'''
        
        to_return = []
        for dup_group in deduplicated_sequences:
            lca = taxonomy_hash[dup_group[0].name]
            for s in dup_group[1:]:
                for i, tax in enumerate(taxonomy_hash[s.name]):
                    if i >= len(lca) or tax != lca[i]:
                        lca = lca[:i]
                        break
                if len(lca) > len(taxonomy_hash[s.name]):
                    lca = lca[:len(taxonomy_hash[s.name])]
            to_return.append(lca)
        return to_return class SequenceIO:
    # Stolen from https://github.com/lh3/readfq/blob/master/readfq.py
    def each(self, fp): # this is a generator function
        last = None # this is a buffer keeping the last unprocessed line
        while True: # mimic closure; is it a bad idea?
            if not last: # the first record or a record following a fastq
                for l in fp: # search for the start of the next record
                    if l[0] in '>@': # fasta/q header line
                        last = l[:-1] # save this line
                        break
            if not last: break
            name, seqs, last = last[1:].partition(" ")[0], [], None
            for l in fp: # read the sequence
                if l[0] in '@+>':
                    last = l[:-1]
                    break
                seqs.append(l[:-1])
            if not last or last[0] != '+': # this is a fasta record
                yield name, ''.join(seqs), None # yield a fasta record
                if not last: break
            else: # this is a fastq record
                seq, leng, seqs = ''.join(seqs), 0, []
                for l in fp: # read the quality
                    seqs.append(l[:-1])
                    leng += len(l) - 1
                    if leng >= len(seq): # have read enough quality
                        last = None
                        yield name, seq, ''.join(seqs); # yield a fastq record
                        break
                if last: # reach EOF before reading enough quality
                    yield name, seq, None # yield a fasta record instead
                    break

    def each_sequence(self, fp):
        '''Like each except iterate over Sequence objects'''
        for name, seq, _ in self.each(fp):
            yield Sequence(name, seq)

    def read_fasta_file(self, path_to_fasta_file):
        seqs = []
        with open(path_to_fasta_file) as f:
            for name, seq, _ in self.each(f):
                seqs.append(Sequence(name, seq))
        return seqs

    def write_fasta_file(self, sequence_objects, path_to_fasta_file):
        with open(path_to_fasta_file,'w') as f:
            self.write_fasta(sequence_objects, f)

    def write_fasta(self, sequence_objects, io):
        for s in sequence_objects:
            io.write(">")
            io.write(s.name)
            io.write("\n")
            io.write(s.seq)
            io.write("\n") class GraftMPackageVersion3(GraftMPackageVersion2):

    version = 3

    def unaligned_sequence_database_path(self):
        if self._contents_hash[GraftMPackage.UNALIGNED_SEQUENCE_DATABASE_KEY]:
            return os.path.join(self._base_directory,
                                self._contents_hash[GraftMPackage.UNALIGNED_SEQUENCE_DATABASE_KEY])
        else:
            return None

    def create_diamond_db(self):
        '''Create a diamond database from the unaligned sequences in this package.

        Returns
        -------
        path to the created diamond db e.g. 'my_sequences.dmnd'
        '''
        base = self.unaligned_sequence_database_path()
        cmd = "diamond makedb --in '%s' -d '%s'" % (self.unaligned_sequence_database_path(), base)
        extern.run(cmd)
        diamondb = '%s.dmnd' % base
        # Mostly this moves a file to it's current location because Create
        # follows this same logic, but there's a specially crafted
        # test/data/mcrA.gpkg which is slightly different.
        os.rename(diamondb, self.diamond_database_path())
        return diamondb

    def taxonomy_hash(self):
        '''Read in the taxonomy and return as a hash of name: taxonomy,
        where taxonomy is an array of strings.'''
        gtns = Getaxnseq()
        with open(self.taxtastic_taxonomy_path()) as tax:
            with open(self.taxtastic_seqinfo_path()) as seqinfo:
                return gtns.read_taxtastic_taxonomy_and_seqinfo(
                    tax, seqinfo)

    @staticmethod
    def graftm_package_is_protein(graftm_package):
        '''Return true if this package is an Amino Acid alignment package, otherwise
        False i.e. it is a nucleotide package. In general it is best to use
        'is_protein_package' instead.

        '''
        found = None
        with open(graftm_package.alignment_hmm_path()) as f:
            r = f.read().split("\n")
        for line in r:
            if line=='ALPH  DNA':
                found = False
                break
            elif line=='ALPH  amino':
                found = True
                break
        if found is None:
            raise Exception("Unable to determine whether the HMM was amino acid or dna")
        return found

    def is_protein_package(self):
        '''Return true if this package is an Amino Acid alignment package, otherwise
        False i.e. it is a nucleotide package. Cache the result for speed.

        '''
        if not hasattr(self, '_is_protein_package'):
            self._is_protein_package = GraftMPackageVersion3.graftm_package_is_protein(self)
        return self._is_protein_package


    @staticmethod
    def compile(output_package_path, refpkg_path, hmm_path, diamond_database_file,
                max_range, unaligned_sequence_database=None,
                trusted_cutoff=False, search_hmm_files=None):
        '''Create a new GraftM package with the given inputs. Any files
        specified as parameters are copied into the final package so can
        be removed after calling this function.

        Parameters
        ----------
        output_package_path: str
            path to the package being created (must not exist)
        refpkg_path: str
            path to pplacer reference package
        hmm_path: str
            path to the search and align HMM
        diamond_database_file: str
            path to diamond DB file, or None for nucleotide packages
        max_rage: str
            as per maximum_range()
        unaligned_sequence_database: str
            path to unaligned sequence database
        trusted_cutoff: boolean
            set TC in search HMM
        search_hmm_files: list of str or None
            use these HMMs for search instead of the hmm_path. All
            basenames of these paths must be unique, and not the same as
            hmm_path.

        Returns
        -------
        Nothing
        '''

        if os.path.exists(output_package_path):
            raise Exception("Not writing new GraftM package to already existing file/directory with name %s" % output_package_path)
        os.mkdir(output_package_path)

        hmm_file_in_gpkg = os.path.basename(hmm_path)
        shutil.copyfile(hmm_path, os.path.join(output_package_path, hmm_file_in_gpkg))
        # Copy unaligned sequence database into graftm package
        if unaligned_sequence_database:
            unaligned_sequence_database_in_gpkg = os.path.join(output_package_path, os.path.basename(unaligned_sequence_database))
            shutil.copyfile(unaligned_sequence_database, unaligned_sequence_database_in_gpkg)

        if diamond_database_file:
            diamond_database_file_in_gpkg = os.path.basename(diamond_database_file)
            shutil.copyfile(diamond_database_file, os.path.join(output_package_path, diamond_database_file_in_gpkg))
        else:
            diamond_database_file_in_gpkg = diamond_database_file
        refpkg_in_gpkg = os.path.basename(refpkg_path)
        shutil.copytree(refpkg_path, os.path.join(output_package_path, refpkg_in_gpkg))

        if search_hmm_files:
            search_hmm_files_in_gpkg_names =\
                [os.path.basename(path) for path in search_hmm_files]
            if hmm_file_in_gpkg in search_hmm_files_in_gpkg_names:
                raise Exception("Search and align HMMs must all have different basenames to create a new gpkg")
            if len(set(search_hmm_files_in_gpkg_names)) != len(search_hmm_files_in_gpkg_names):
                raise Exception("Search HMMs must have different basenames to create a new gpkg")
            for i, search_hmm in enumerate(search_hmm_files):
                shutil.copyfile(search_hmm, os.path.join(output_package_path, search_hmm_files_in_gpkg_names[i]))
        else:
            search_hmm_files_in_gpkg_names = [hmm_file_in_gpkg]

        contents = {GraftMPackage.VERSION_KEY: GraftMPackageVersion3.version,
                    GraftMPackage.ALIGNMENT_HMM_KEY: hmm_file_in_gpkg,
                    GraftMPackage.SEARCH_HMM_KEY: search_hmm_files_in_gpkg_names,
                    GraftMPackage.REFERENCE_PACKAGE_KEY: refpkg_in_gpkg,
                    GraftMPackage.HMM_TRUSTED_CUTOFF_KEY: trusted_cutoff,
                    GraftMPackage.RANGE_KEY: max_range,
                    GraftMPackage.UNALIGNED_SEQUENCE_DATABASE_KEY: (os.path.basename(unaligned_sequence_database)
                                                                    if unaligned_sequence_database else None),
                    GraftMPackage.DIAMOND_DATABASE_KEY: diamond_database_file_in_gpkg}

        with open(os.path.join(
                output_package_path, GraftMPackage._CONTENTS_FILE_NAME),
                  'w') as f:
            json.dump(contents, f) class GraftMPackage:
    '''Package to represent a GraftM package. To start using a package, run

    pkg = GraftMPackage.acquire('/path/to/a.gpkg')

    and then retrieve values with e.g.

    pkg.alignment_hmm_path()
      #=> '/path/to/a.gpkg/the.hmm'
    '''

    _CONTENTS_FILE_NAME = 'CONTENTS.json'

    # The key names are unlikely to change across package format versions,
    # so store them here in the superclass
    DIAMOND_DATABASE_KEY = 'diamond_database'
    VERSION_KEY = 'graftm_package_version'
    ALIGNMENT_HMM_KEY = 'align_hmm'
    SEARCH_HMM_KEY = "search_hmms"
    REFERENCE_PACKAGE_KEY = "refpkg"
    HMM_TRUSTED_CUTOFF_KEY = "trusted_cutoff"
    RANGE_KEY = "range"
    UNALIGNED_SEQUENCE_DATABASE_KEY = "unaligned_sequence_database"
    _CONTENTS_FILE_NAME = 'CONTENTS.json'

    _CURRENT_VERSION = 3

    _REQUIRED_KEYS = {'2': [
                             VERSION_KEY,
                             ALIGNMENT_HMM_KEY,
                             SEARCH_HMM_KEY,
                             REFERENCE_PACKAGE_KEY,
                             HMM_TRUSTED_CUTOFF_KEY,
                             RANGE_KEY
                             ],
                      '3': [
                             VERSION_KEY,
                             ALIGNMENT_HMM_KEY,
                             SEARCH_HMM_KEY,
                             REFERENCE_PACKAGE_KEY,
                             HMM_TRUSTED_CUTOFF_KEY,
                             RANGE_KEY,
                             UNALIGNED_SEQUENCE_DATABASE_KEY
                             ]
                      }


    @staticmethod
    def acquire(graftm_package_path):
        '''Acquire a new graftm Package

        Parameters
        ----------
        graftm_output_path: str
            path to base directory of graftm
        '''

        with open(os.path.join(
                graftm_package_path,
                GraftMPackage._CONTENTS_FILE_NAME
                )) as f:
            contents_hash = json.load(f)

        v=contents_hash[GraftMPackage.VERSION_KEY]
        logging.debug("Loading version %i GraftM package: %s" % (v, graftm_package_path))
        if v == 2:
            pkg = GraftMPackageVersion2()
        elif v == 3:
            pkg = GraftMPackageVersion3()
        else:
            raise InsufficientGraftMPackageException("Bad version: %s" % v)

        pkg._contents_hash = contents_hash
        pkg._base_directory = graftm_package_path
        # check we are at current version otherwise choke
        pkg.check_universal_keys(v)
        pkg.check_required_keys(GraftMPackage._REQUIRED_KEYS[str(v)])
        return pkg

    def check_universal_keys(self, version):
        h = self._contents_hash
        try:
            v = h[GraftMPackage.VERSION_KEY]
        except KeyError:
            raise InsufficientGraftMPackageException("No version information in graftm package")
        if v != version:
            raise InsufficientGraftMPackageException("Bad version: %s" % v)

    def check_required_keys(self, required_keys):
        '''raise InsufficientGraftMPackageException if this package does not
        conform to the standard of the given package'''
        h = self._contents_hash
        for key in required_keys:
            if key not in h:
                raise InsufficientGraftMPackageException("Package missing key %s" % key)

    def __getitem__(self, key):
        '''Return the value of the given key from the contents file'''
        return self.contents_hash[key]

    def contents_file_path(self):
        return os.path.join(self._base_directory, GraftMPackage._CONTENTS_FILE_NAME) class Decorator:
    '''re-root a tree and decorate it using old taxonomy, for the graftM 
    decorate pipeline.
    
    ALL trees provided to this class must be in newick format. No other formats 
    are supported. If any taxonomic decoration is already in the tree to be 
    decorated by this class is WILL be overwritten. Of course, bootstrap values 
    will remain untouched.'''

    def __init__(self, **kwargs):
        '''
        Parameters
        ----------
        reference_tree_path: str
            Path to the file containing the reference tree, which is used to
            retroot the tree tree provided to tree
        tree_path: str
            Path to the file containing the tree to be re-rooted. This tree will
            be rerooted at the same position as the tree porovided to the 
            reference_tree
        '''
        reference_tree_path = kwargs.pop('reference_tree_path', None)
        tree_path = kwargs.pop('tree_path')
        
        logging.debug("Importing old tree from file: %s" 
                        % tree_path)
        self.tree = Tree.get(path=tree_path,
                            schema='newick')
        if reference_tree_path:
            logging.debug("Importing reference tree from file: %s" 
                            % reference_tree_path)
            self.reference_tree = Tree.get(path=reference_tree_path,
                                           schema='newick')
        else:
            self.reference_tree = reference_tree_path
        if len(kwargs) > 0:
            raise Exception("Unexpected arguments provided to Decorator class: %s" % kwargs)
        
    def _reroot(self):
        '''Run the re-rooting algorithm in the Rerooter class.'''
        rerooter = Rerooter()
        self.tree = rerooter.reroot_by_tree(self.reference_tree, 
                                                self.tree)
        
    def main(self, taxonomy, output_tree, output_tax, no_unique_tax,
             decorate, seqinfo):
        '''Decorate and if necessary, re-root the tree. If an old reference tree
        is provided it is assumed that re-rooting is desired
        
        Parameters
        ----------
        taxonomy: str
            Path to a GreenGenes formatted, or taxtastic formatted taxonomy 
            file, containing the taxonomy of all or some if the sequences used 
            to construct the tree. This taxonomy will be used by the 
            tree_decorator.TreeDecorator class to decorate self.tree
        output_tree: str
            Path to file to which the decorated tree will be written to.
        output_tax: str
            Path to file to which the decorated taxonomy will be written to.
        seqinfo: str
            path to taxtastic seqinfo file to be used with taxonomy to annotate
            the tree.'''
        # Reroot
        if self.reference_tree:
            self._reroot()
        # Decorate
        td =  TreeDecorator(self.tree,
                            taxonomy,
                            seqinfo)
        if decorate:
            td.decorate(output_tree, output_tax, no_unique_tax)
        else:
            logging.info("Writing tree to file: %s" % (output_tree))
            self.tree.write(
                            output_tree,
                            format = "newick"   
                            ) class GreenGenesTaxonomy:
    # A hash of name to array of taxonomy strings
    taxonomy = {}

    def __init__(self, taxonomy={}):
        self.taxonomy = taxonomy

    @staticmethod
    def read(input_taxonomy_io):
        '''Parse in a taxonomy from the given I/O.

        Parameters
        ----------
        input_taxonomy_io: io
            an open io object that is iterable. This stream is neither opened
            nor closed by this method

        Returns
        -------
        A GreenGenesTaxonomy instance with the taxonomy parsed in

        Raises
        ------
        DuplicateTaxonomyException: if there is duplicate taxonomy entries
        MalformedGreenGenesTaxonomyException: if there is something else amiss'''
        tax = {}
        for line in input_taxonomy_io:
            if len(line.strip()) == 0: continue #ignore empty lines
            splits = line.split("\t")
            if len(splits) != 2:
                raise MalformedGreenGenesTaxonomyException("Unexpected number of tab-separated fields found in taxonomy file, on line %s" % line)
            name = splits[0].strip()
            taxonomy = [t.strip() for t in splits[1].split(';')]
            while len(taxonomy) > 0 and taxonomy[-1] == '':
                taxonomy = taxonomy[:-1]
            for lineage in taxonomy:
                if lineage == '':
                    raise MalformedGreenGenesTaxonomyException("Encountered a taxonomy string with the middle of the taxonomy string missing: %s" % line)

            if name in tax:
                raise DuplicateTaxonomyException("Duplicate definition of taxonomy for %s" % name)
            tax[name] = taxonomy
        return GreenGenesTaxonomy(tax)

    @staticmethod
    def read_file(input_filename):
        '''Like read() except uses a file rather than an IO stream, for convenience'''
        with open(input_filename) as f:
            g = GreenGenesTaxonomy.read(f)
        return g

    def write(self, output_io):
        '''Write a taxonomy to an open stream out in GG format. Code calling this
        function must open and close the io object.'''
        for name, tax in self.taxonomy.items():
            output_io.write("%s\t%s\n" % (name, '; '.join(tax))) class SequenceSearcher:

    def __init__(self, search_hmm, aln_hmm=None):
        self.search_hmm = search_hmm
        self.aln_hmm = aln_hmm

    def _get_sequence_directions(self, search_result):
        sequence_directions = {}
        for result in search_result:
            for hit in result.each([SequenceSearchResult.QUERY_ID_FIELD, SequenceSearchResult.ALIGNMENT_DIRECTION]):
                sequence_directions[hit[0]]  ={"strand":[hit[1]], "entry": []}
        return sequence_directions


    def _hmmalign(self, input_path, directions, pipeline,
                  forward_reads_output_path, reverse_reads_output_path):
        '''
        Align reads to the aln_hmm. Receives unaligned sequences and
        aligns them.

        Parameters
        ----------
        input_path : str
            Filename of unaligned hits to be aligned
        directions : dict
            dictionary containing read names as keys, and complement
            as the entry (True=Forward, False=Reverse)
        pipeline: str
            either PIPELINE_AA = "P" or PIPELINE_NT = "D"
        forward_reads_output_fh: str
            Where to write aligned forward reads
        reverse_reads_output_fh: str
            Where to write aligned reverse reads
        Returns
        -------
        Nothing.
        '''
        if pipeline == PIPELINE_AA:
            reverse_direction_reads_present=False
        else:
            reverse_direction_reads_present=False in directions.values()

        with tempfile.NamedTemporaryFile(prefix='for_file', suffix='.fa') as for_file_fh:
            for_file = for_file_fh.name
            with tempfile.NamedTemporaryFile(prefix='rev_file', suffix='.fa') as rev_file_fh:
                rev_file = rev_file_fh.name
                # Align input reads to a specified hmm.
                if reverse_direction_reads_present:  # Any that are in the reverse direction would be True
                    reverse = []
                    forward = []
                    records = list(SeqIO.parse(open(input_path), 'fasta'))

                    # Split the reads into reverse and forward lists
                    for record in records:
                        read_id = record.id

                        if directions[read_id] == True:
                            forward.append(record)
                        elif directions[read_id] == False:
                            reverse.append(record)
                        else:
                            raise Exception(logging.error('Programming error: hmmalign'))
                            exit(1)

                    logging.debug("Found %i forward direction reads" % len(forward))
                    logging.debug("Found %i reverse direction reads" % len(reverse))

                    # Write reverse complement and forward reads to files
                    with open(for_file, 'w') as for_aln:
                        logging.debug("Writing forward direction reads to %s" % for_file)
                        for record in forward:
                            for_aln.write('>' + record.id + '\n')
                            for_aln.write(str(record.seq) + '\n')
                            # HMMalign and convert to fasta format
                    if any(forward):
                        self.hmmalign_sequences(self.aln_hmm, for_file, forward_reads_output_path)
                    else:
                        cmd = 'touch %s' % (forward_reads_output_path)
                        extern.run(cmd)
                    with open(rev_file, 'w') as rev_aln:
                        logging.debug("Writing reverse direction reads to %s" % rev_file)
                        for record in reverse:
                            if record.id and record.seq:
                                rev_aln.write('>' + record.id + '\n')
                                rev_aln.write(str(record.seq.reverse_complement()) + '\n')
                    self.hmmalign_sequences(self.aln_hmm, rev_file, reverse_reads_output_path)
                    conv_files = [forward_reads_output_path, reverse_reads_output_path]
                    return conv_files

                else:
                    # If there are only forward reads, just hmmalign and be done with it.
                    self.hmmalign_sequences(self.aln_hmm, input_path, forward_reads_output_path)
                    conv_files = [forward_reads_output_path]

                    return conv_files

    def hmmalign_sequences(self, hmm, sequences, output_file):
        '''Run hmmalign and convert output to aligned fasta format

        Parameters
        ----------
        hmm: str
            path to hmm file
        sequences: str
            path to file of sequences to be aligned
        output_file: str
            write sequences to this file

        Returns
        -------
        nothing
        '''

        cmd = 'hmmalign --trim %s %s' % (hmm, sequences)
        output = extern.run(cmd)
        with open(output_file, 'w') as f:
            SeqIO.write(SeqIO.parse(StringIO(output), 'stockholm'), f, 'fasta')

    def makeSequenceBinary(self, sequences, fm):
        cmd = 'makehmmerdb %s %s' % (sequences, fm)
        extern.run(cmd)

    def hmmsearch(self, output_path, input_path, unpack, seq_type, threads, cutoff, orfm):
        '''
        hmmsearch - Search raw reads for hits using search_hmm list

        Parameters
        ----------
        output_path : str
            path to output domtblout table
        input_path : str
            path to input sequences to search
        unpack : UnpackRawReads
            UnpackRawReads object, returns string command that will output
            sequences to stdout when called on command line
            (use: unpack.command_line())
        seq_type : str
            variable containing a string, either 'nucleotide' or 'aminoacid'.
            Tells the pipeline whether or not to call ORFs on the sequence.
            If sequence is 'nucleotide', ORFs are called. If not, no ORFs.
        threads : Integer
            Number of threads to use. Passed to HMMsearch command.
        cutoff : str
            cutoff for HMMsearch to use, either an evalue or --cut_tc, meaning
            use the TC score cutoff specified within the HMM. Passed to
            HMMsearch command.
        orfm : OrfM
            Object that builds the command chunk for calling ORFs on sequences
            coming through as stdin. Outputs to stdout. Calls command_line
            to construct final command line string.

        Returns
        -------
        output_table_list : array of HMMSearchResult
            Includes the name of the output domtblout table given by hmmer

        Raises
        ------
        hmmsearcher.NoInputSequencesException
            Raised if there are no sequences fed into the HMM.
        '''

        # Define the base hmmsearch command.
        logging.debug("Using %i HMMs to search" % (len(self.search_hmm)))
        output_table_list = []
        if len(self.search_hmm) > 1:
            for hmm in self.search_hmm:
                out = os.path.join(os.path.split(output_path)[0], os.path.basename(hmm).split('.')[0] + '_' + os.path.split(output_path)[1])
                output_table_list.append(out)
        elif len(self.search_hmm) == 1:
            output_table_list.append(output_path)
        else:
            raise Exception("Programming error: expected 1 or more HMMs")

        # Choose an input to this base command based off the file format found.
        if seq_type == 'nucleotide':  # If the input is nucleotide sequence
            input_cmd = orfm.command_line(input_path)
        elif seq_type == 'aminoacid':  # If the input is amino acid sequence
            input_cmd = unpack.command_line()
        else:
            raise Exception('Programming Error: error guessing input sequence type')

        # Run the HMMsearches
        if cutoff == "--cut_tc":
            searcher = HmmSearcher(threads, cutoff)
        else:
            searcher = HmmSearcher(threads, '--domE %s' % cutoff)
        searcher.hmmsearch(input_cmd, self.search_hmm, output_table_list)

        hmmtables = [HMMSearchResult.import_from_hmmsearch_table(x) for x in output_table_list]
        return hmmtables

    def merge_forev_aln(self, forward_aln_list, reverse_aln_list, outputs):
        '''
        merge_forev_aln - Merges forward and reverse alignments for a given run

        Parameters
        ----------
        forward_aln_list : array
            List of the forward or interleaved alignments for each of the runs
            given to graftM. 
        reverse_aln_list : array
            List of the reverse alignments for each of the runs
            given to graftM. 
            **MUST** be the the same length as forward_aln_list or empty (if
            interleaved)
        outputs : array
            List of paths to output file to which the merged alignments from the
            aln_list will go into. Must be exactly half the size of the aln_list
            (i.e. one output file for every forward and reverse file provided)

        Returns
        -------
        Nothing - output files are known.
        '''

        orfm_regex = OrfM.regular_expression()
        def remove_orfm_end(id):
            """ remove the orfm suffix from a read id """
            try:
                return orfm_regex.match(id).group(1)
            except AttributeError:
                return id

        dir_regex = re.compile(r'/[12]$')
        def trim_direction(id):
            """ remove fwd/rev suffix if present """
            return dir_regex.sub('',id)

        def get_read_base(id):
            return trim_direction(remove_orfm_end(id))

        def get_read_bases(record_iter):
            """ remove orfm suffix from keys in dict records """
            return {get_read_base(r.id):r for r in record_iter}

        def split_interleaved_reads(records):
            """ return list of fwd records and dict of reverse """
            fwd = []
            rev = {}
            for record in records:
                id = get_read_base(record.id)
                if remove_orfm_end(record.id).endswith('1'):
                    fwd.append(record)
                else:
                    rev[id] = record
            return fwd, rev

        if not reverse_aln_list:
            # create dummy array for interleaved
            reverse_aln_list = [False,]*len(forward_aln_list)

        for idx, (forward_path, reverse_path) in enumerate(zip(forward_aln_list, reverse_aln_list)):
            output_path = outputs[idx]
            logging.info('Merging pair %s, %s' %
                         (os.path.basename(forward_path),
                          os.path.basename(reverse_path) \
                          if reverse_path else None))
            if reverse_path:
                forward_reads = SeqIO.parse(forward_path, 'fasta')
                reverse_reads = get_read_bases(SeqIO.parse(reverse_path, 'fasta'))
            else:
                forward_reads, reverse_reads = \
                        split_interleaved_reads(SeqIO.parse(forward_path,
                                                            'fasta'))
                #logging.debug("Loaded %d fwd (EG %s) and %d rev (EG %s)",
                #              len(forward_reads),
                #              get_read_base(next(iter(forward_reads)).id),
                #              len(reverse_reads),
                #              next(iter(reverse_reads)))

            with open(output_path, 'w') as out:
                for forward_record in forward_reads:
                    id = get_read_base(forward_record.id)
                    forward_sequence = str(forward_record.seq)
                    logging.debug("Fwd id %s is %s in reverse dict",
                                  id, "" if id in reverse_reads else "not")
                    try:
                        reverse_sequence = str(reverse_reads[id].seq)
                        new_seq = ''
                        if len(forward_sequence) == len(reverse_sequence):
                            for f, r in zip(forward_sequence, reverse_sequence):
                                if f == r:
                                    new_seq += f
                                elif f == '-' and r != '-':
                                    new_seq += r
                                elif r == '-' and f != '-':
                                    new_seq += f
                                elif f != '-' and r != '-':
                                    if f != r:
                                        new_seq += f
                                else:
                                    new_seq += '-'
                        else:
                            logging.error('Alignments do not match')
                            raise Exception('Merging alignments failed: Alignments do not match')
                        out.write('>%s\n' % forward_record.id)
                        out.write('%s\n' % (new_seq))
                        del reverse_reads[id]
                    except:

                        out.write('>%s\n' % forward_record.id)
                        out.write('%s\n' % (forward_sequence))
                for record_id, record in reverse_reads.items():
                    out.write('>%s\n' % record.id)
                    out.write('%s\n' % (str(record.seq)))

    def nhmmer(self, output_path, unpack, threads, evalue):
        '''
        nhmmer - Search input path using nhmmer

        Parameters
        ----------
        output_path : str
            A string containing the path to the input sequences
        unpack : obj
            UnpackRawReads object, returns string command that will output
            sequences to stdout when called on command line
            (use: unpack.command_line())
        threads : str
            Number of threads to run. For compiling command line.
        evalue : str
            evalue to use. For compiling commmand line.

        Returns
        -------
        output_table_list : array
            Includes the name of the output domtblout table given by hmmer
        '''
        logging.debug("Using %i HMMs to search" % (len(self.search_hmm)))
        output_table_list = []
        if len(self.search_hmm) > 1:
            for hmm in self.search_hmm:
                out = os.path.join(os.path.split(output_path)[0], os.path.basename(hmm).split('.')[0] + '_' + os.path.split(output_path)[1])
                output_table_list.append(out)
        elif len(self.search_hmm) == 1:
            output_table_list.append(output_path)
        else:
            raise Exception("Programming error: Expected 1 or more HMMs")
        input_pipe = unpack.command_line()

        searcher = NhmmerSearcher(threads, extra_args='--incE %s -E %s' % (evalue, evalue))
        searcher.hmmsearch(input_pipe, self.search_hmm, output_table_list)

        hmmtables = [HMMSearchResult.import_from_nhmmer_table(x) for x in output_table_list]

        return hmmtables, output_table_list

    def _check_euk_contamination(self, hmm_hit_tables):
        '''
        check_euk_contamination - Check output HMM tables hits reads that hit
                                  the 18S HMM with a higher bit score.

        Parameters
        ----------
        hmm_hit_tables : array
            Array of paths to the output files produced by hmmsearch or
            nhmmer.
        run_stats : dict
            A dictionary to updatewith the number of unique 18S reads and reads
            detected by both 18S and non-18S HMMs

        Returns
        -------
        euk_reads : set
            Non-redundant set of all read names deemed to be eukaryotic
        '''

        euk_hit_table = HMMreader(hmm_hit_tables.pop(-1))
        other_hit_tables = [HMMreader(x) for x in hmm_hit_tables]
        reads_unique_to_eukaryotes = []
        reads_with_better_euk_hit = []

        for hit in euk_hit_table.names():
            bits = []
            for hit_table in other_hit_tables:
                if hit in hit_table.names():
                    bits.append(hit_table.bit(hit))
                else:
                    reads_unique_to_eukaryotes.append(hit)
            if bits:
                if any([x for x in bits if x > euk_hit_table.bit(hit)]):
                    continue
                else:
                    reads_with_better_euk_hit.append(hit)
            else:
                continue

        if len(reads_with_better_euk_hit) == 0:
            logging.info("No contaminating eukaryotic reads detected")
        else:
            logging.info("Found %s read(s) that may be eukaryotic" % len(reads_with_better_euk_hit + reads_unique_to_eukaryotes))

        euk_reads = set(reads_with_better_euk_hit + reads_unique_to_eukaryotes)

        return euk_reads


    def _extract_multiple_hits(self, hits, reads_path, output_path):
        '''
        splits out regions of a read that hit the HMM. For example when two of
        same gene are identified within the same contig, The regions mapping to
        the HMM will be split out and written out to a new file as a new record.

        Parameters
        ----------
        hits : dict
            A dictionary where the keys are the read names, the entry for each
            is a list of lists, each containing the range within the contig
            (or read) that mapped to the HMM. e.g.:

            {'read1': [[3, 126], [305, 413]],
             'read2': [[1, 125]],
             ...

        reads_path : str
            path to reads file containing each read or contig in FASTA format to
            be opened, and split.
        output_path : str
            path to file to which split reads will be written to in FASTA
            format.

        Returns
        -------
        Nothing, output path is known.
        '''

        complement_information = {}
        try:
            reads = SeqIO.to_dict(SeqIO.parse(reads_path, "fasta"))  # open up reads as dictionary
        except:
            logging.error("Multiple sequences found with the same ID. The input sequences are either ill formated or are interleaved. \
If you provided GraftM with an interleaved sequence file, please split them into forward and reverse reads, and provide to the the appropriate \
flags (--forward, --reverse). Otherwise, it appears that you have provided sequences with redundant IDs. GraftM doesn't know how to \
deal with these, so please remove/rename sequences with duplicate keys.")
            raise InterleavedFileError()
        with open(output_path, 'w') as out:
            for read_name, entry in hits.items():  # For each contig
                ranges = entry["entry"]
                complements = entry["strand"]
                index = 1
                if len(ranges) > 1:  # if there are multiple hits in that contig
                    for r, c in zip(ranges, complements):  # for each of those hits

                        new_record = reads[read_name][r[0] - 1:r[1]]  # subset the record by the span of that hit
                        new_record.id = new_record.id + '_split_%i' % index  # give that subset record a new header
                        SeqIO.write(new_record, out, "fasta")  # and write it to output
                        index += 1  # increment the split counter
                        complement_information[new_record.id]=c

                else:  # Otherwise, just write the read back to the file
                    complement_information[read_name] = entry["strand"][0]
                    SeqIO.write(reads[read_name], out, "fasta")

        return complement_information

    def _extract_from_raw_reads(self, output_path, input_reads, raw_sequences_path, input_file_format, hits):
        '''
        _extract_from_raw_reads - Extract hit sequences of the hmm/diamond
        search from a command which generates uncompressed FASTA sequences from
        the raw sequences file. Output into specified file

        Parameters
        ----------
        output_path : str
            Path of the desired output file
        input_reads : list
            list of read IDS to extract.
        raw_sequences_path : str
            Path to the raw sequences
        input_file_format : var
            Variable, either FORMAT_FASTA_GZ, FORMAT_FASTQ_GZ, FORMAT_FASTQ or
            FORMAT_FASTA, denoting the format of the input sequence
        hits : dict
            A hash with the readnames as the keys and the spans as the values

        Returns
        -------
        output_path: str
            Path to file containing extracted reads.

        direction_information: dict

            {read_1: False
             ...
             read n: True}

            where True = Forward direction
            and False = Reverse direction
        '''

        with tempfile.NamedTemporaryFile(prefix='_raw_extracted_reads.fa') as tmp:
            # Extract reads from original sequence file
            extract_cmd = "mfqe --output-uncompressed"
            if input_file_format in (FORMAT_FASTA, FORMAT_FASTA_GZ, FORMAT_FASTQ, FORMAT_FASTQ_GZ):
                extract_cmd += " --fasta-read-name-lists /dev/stdin --input-fasta {} --output-fasta-files '{}'".format(
                    raw_sequences_path, tmp.name)
            else:
                raise Exception("Programming error: Unexpected input file format {}".format(input_file_format))

            extern.run(extract_cmd, stdin='\n'.join(input_reads))
            complement_info = self._extract_multiple_hits(hits, tmp.name, output_path)  # split them into multiple reads

        return output_path, complement_info


    def alignment_correcter(self, alignment_file_list, output_file_name,
                            filter_minimum=None):
        '''
        Remove lower case insertions in alignment outputs from HMM align. Give
        a list of alignments, and an output file name, and each alignment will
        be corrected, and written to a single file, ready to be placed together
        using pplacer.

        Parameters
        ----------
        alignment_file_list : array
            List of strings, each the path to different alignments from the
            inputs provided to GraftM
        output_file_name : str
            The path and filename of the output file desired.
        filter_minimum : int
            minimum number of positions that must be aligned for each sequence
        Returns
        -------
        True or False, depending if reads were written to file

        '''

        corrected_sequences = {}
        for alignment_file in alignment_file_list:
            insert_list = []  # Define list containing inserted positions to be removed (lower case characters)
            with open(alignment_file) as f:
                sequence_list = list(SeqIO.parse(f, 'fasta'))
            for sequence in sequence_list:  # For each sequence in the alignment
                for idx, nt in enumerate(list(sequence.seq)):  # For each nucleotide in the sequence
                    if nt.islower():  # Check for lower case character
                        insert_list.append(idx)  # Add to the insert list if it is
            insert_list = list(OrderedDict.fromkeys(sorted(insert_list, reverse=True)))  # Reverse the list and remove duplicate positions
            for sequence in sequence_list:  # For each sequence in the alignment
                new_seq = list(sequence.seq)  # Define a list of sequences to be iterable list for writing
                for position in insert_list:  # For each position in the removal list
                    del new_seq[position]  # Delete that inserted position in every sequence
                corrected_sequences['>' + sequence.id + '\n'] = (''.join(new_seq) + '\n').replace('~', '-')

        pre_filter_count=len(corrected_sequences)

        if filter_minimum:
            # Use '>' not '>=' here because the sequence is on a single line,
            # but also includes a newline character at the end of the sequence
            corrected_sequences={key:item for key, item in iter(corrected_sequences.items()) if len(item.replace('-', '')) > filter_minimum}

        post_filter_count=len(corrected_sequences)
        logging.info("Filtered %i short sequences from the alignment" % \
                        (pre_filter_count-post_filter_count)
                    )
        logging.info("%i sequences remaining" % post_filter_count)

        if len(corrected_sequences) >= 1:
            with open(output_file_name, 'w') as output_file:  # Create an open file to write the new sequences to
                for fasta_id, fasta_seq in corrected_sequences.items():
                    output_file.write(fasta_id)
                    output_file.write(fasta_seq)
            return True
        else:
            return False


    def _extract_orfs(self, input_path, orfm, hit_readnames, output_path, search_method, sequence_frame_info_list=None):
        '''
        Call ORFs on a file with nucleotide sequences and extract the proteins
        whose name is in `hit_readnames`.

        Parameters
        ----------
        input_path : str
            Path to input nucleotide sequences in FASTA format.
        orfm: obj
            graftm.OrfM object with parameters already set
        hit_readnames : str
            path to a file containin the readnames of hits to the HMM, one per
            line.
        output_path : str
            Path to output orfs into, in FASTA format.
        search_method : str
            The method for searching, either 'hmmsearch' or 'diamond'
        sequence_frame_info : list
            A dataframe (list of lists) containing readname, alignment direction
            and alignment start point information
        '''

        if search_method == "hmmsearch":
            # Build and run command to extract ORF sequences:
            orfm_cmd = orfm.command_line()
            cmd = "mfqe --output-uncompressed --fasta-read-name-lists /dev/stdin --input-fasta <({} {}) --output-fasta-files {}".format(
                orfm_cmd, input_path, output_path)
            extern.run(cmd, stdin='\n'.join(hit_readnames))

        elif search_method == "diamond":
            sequence_frame_info_dict = {x[0]:[x[1], x[2], x[3]] for x in sequence_frame_info_list}
            records = SeqIO.parse(input_path, "fasta")
            with open(output_path, 'w') as open_output_path:
                for record in records:
                    entry=sequence_frame_info_dict[record.id]
                    indfrom=(min(entry[2], entry[1])-1)
                    indto=max(entry[2], entry[1])
                    if entry[0] == False:
                        record.seq = max([x for x in record.seq[indfrom:indto].reverse_complement().translate().split("*")], key=len)
                    else:
                        record.seq = max([x for x in record.seq[indfrom:indto].translate().split("*")], key=len)
                    SeqIO.write(record, open_output_path, "fasta")
                open_output_path.flush()


    def _get_read_names(self, search_result, max_range):
        '''
        _get_read_names - loops through hmm hits and their alignment spans to
        determine if they are potentially linked (for example, if one gene in a
        contig hits a hmm more than once, in two different conserved regions
        of that gene) and combines them into one 'hit' because they are
        technically the same. The total span of the hits deemed to be linked is
        returned.

        Parameters
        ----------
        search_result : obj
            SequenceSearchResult object with all paramaters defined. Used here
            to create rows containing information on alignment direction and
            alignment span.
        max_range : int
            Maximum range that a gene can extend within a contig. Any hits
            that extend beyond this length cannot be linked. max_range is
            set as 1.5 X the average length of all full length genes used
            in the search database. This is defined in the CONTENTS.json file
            within a gpkg.
        Returns
        -------
            Dictionary where keys are the contig/read name. The value for each
            entry is an array lists, one per hit in each contig, each with the
            span (min and max) of the alignment.
        '''

        splits = {}  # Define an output dictionary to be filled
        spans = []
        for result in search_result:  # Create a table (list of rows contain span, and complement information
            spans += list(
                         result.each(
                                    [SequenceSearchResult.QUERY_ID_FIELD,
                                        SequenceSearchResult.ALIGNMENT_DIRECTION,
                                        SequenceSearchResult.HIT_FROM_FIELD,
                                        SequenceSearchResult.HIT_TO_FIELD,
                                        SequenceSearchResult.QUERY_FROM_FIELD,
                                        SequenceSearchResult.QUERY_TO_FIELD]
                                       )
                         )

        for hit in spans:  # For each of these rows (i.e. hits)
            i = hit[0]  # set id to i
            c = hit[1]  # set complement to c
            ft = [min(hit[2:4]), max(hit[2:4])]  # set span as ft (i.e. from - to) - This is the amount covering the query
            qs = [min(hit[4:6]), max(hit[4:6])]  # seq the query span to qs - This is the amount covering the HMM

            if ft[0] == ft[1]: continue  # if the span covers none of the query, skip that entry (seen this before)

            if i not in splits:  # If the hit hasnt been seen yet
                splits[i] = {'span'       : [ft],
                             'strand'     : [c],
                             'query_span' : [qs]}  # add the span and complement as new entry

            else:  # otherwise (if it has been seen)
                for idx, entry in enumerate(splits[i]['span']):  # for each previously added entry
                    if splits[i]['strand'][idx] != c:  # If the hit is on the same complement strand
                        splits[i]['span'].append(ft)  # Add the new range to be split out in the future
                        splits[i]['strand'].append(c)  # Add the complement strand as well
                        splits[i]['query_span'].append(qs)
                        break

                    previous_qs      = splits[i]['query_span'][idx] # Get the query span of the previous hit

                    previous_q_range = set(range(previous_qs[0], previous_qs[1])) # Get the range of each
                    current_q_range  = set(range(qs[0], qs[1]))
                    query_overlap    = set(previous_q_range).intersection(current_q_range) # Find the intersection between the two ranges

                    previous_ft_span = set(range(entry[0], entry[1]))
                    current_ft_span  = set(range(ft[0], ft[1]))
                    if any(query_overlap): # If there is an overlap
                        ####################################################
                        # if the span over the actual read that hit the HMM
                        # for each hit overlap by > 25%, they are considered
                        # the same hit, and ignored
                        ####################################################
                        intersection_fraction = float(len(previous_ft_span.intersection(current_ft_span)))

                        if intersection_fraction / float(len(previous_ft_span)) >= PREVIOUS_SPAN_CUTOFF:
                            break
                        elif intersection_fraction / float(len(current_ft_span)) >= PREVIOUS_SPAN_CUTOFF:
                            break
                        else: # else (i.e. if the hit covers less that 25% of the sequence of the previous hit)
                            ####################################################
                            # If they made it this far, it means that the hits do not overlap.
                            # But one last check must be made to ensure they do not cover the same
                            # region in the HMM.
                            ####################################################
                            if len(query_overlap) > (len(current_q_range)*PREVIOUS_SPAN_CUTOFF): # if the overlap on the query HMM does not span over 25%
                                if (idx+1) == len(splits[i]['span']):
                                    splits[i]['span'].append(ft) # Add from-to as another entry, this is another hit.
                                    splits[i]['strand'].append(c) # Add strand info as well
                                    splits[i]['query_span'].append(qs)
                                    break # And break

                    if min(entry) < min(ft):  # if/else to determine which entry comes first (e.g. 1-5, 6-10 not 6-10, 1-5)
                        if max(ft) - min(entry) < max_range:  # Check if they lie within range of eachother
                            entry[1] = max(ft)  # ammend the entry if they are
                            break  # And break the loop
                    else:
                        if max(entry) - min(ft) < max_range:  # Check if they lie within range of eachother
                            entry[0] = min(ft)  # ammend the entry if they are
                            break  # And break the loop

                else:  # if no break occured (no overlap)
                    splits[i]['span'].append(ft)  # Add the new range to be split out in the future
                    splits[i]['strand'].append(c)  # Add the complement strand as well
                    splits[i]['query_span'].append(qs)
        return {key: {"entry":entry['span'], 'strand': entry['strand']} for key, entry in iter(splits.items())}  # return the dict, without strand information which isn't required.

    def _check_for_slash_endings(self, readnames):
        '''
        Provide a list of read names to be checked for the /1 or /2 endings
        which cause troubles downstream when working with paired reads. This
        function checks the list for ANY read ending in /1 or /2. The reads are
        first stripped of any comment lines, using the 'split' function

        Parameters
        ----------
        readnames : list
            list of strings, each a read name to be checked for the /1 or /2
            ending

        Returns
        -------
        Boolean. True if /1 or /2 endings are found, False if not.

        '''

        # Strip comment lines, if any
        readnames = [x.split()[0] for x in readnames]

        # Check for /1 and /2 endings:
        if any(x for x in readnames if x.endswith('/1') or x.endswith('/2')):
            result = True
        else:
            result = False

        return result

    @T.timeit
    def aa_db_search(self, files, base, unpack, search_method,
                     maximum_range, threads, evalue, min_orf_length,
                     restrict_read_length, diamond_database,
                     diamond_performance_parameters):
        '''
        Amino acid database search pipeline - pipeline where reads are searched
        as amino acids, and hits are identified using hmmsearch or diamond
        searches

        Parameters
        ----------
        files : obj
            graftm_output_paths object.
        base : str
            The name of the input file, stripped of all suffixes, and paths.
            Used for creating file names with 'files' object.
        unpack : obj
            UnpackRawReads object, returns string command that will output
            sequences to stdout when called on command line
            (use: unpack.command_line())
        search_method : str
            The method for searching, either 'hmmsearch' or 'diamond'
        maximum_range : int
            Maximum range that a gene can extend within a contig. Any hits
            that extend beyond this length cannot be linked. max_range is defined
            as 1.5 X the average length of all full length genes used in the
            search database. This is defined in the CONTENTS.json file within a
            gpkg.
        threads : int
            Number of threads for hmmer to use
        evalue : str
            evalue cutoff for hmmer to use
        min_orf_length : int
            minimum orf length for orfm to use
        restrict_read_length : int
            orf length to retrict orfm to.
        diamond_database : str
            Path to diamond database to use when searching. Set to 'None' if not
            using diamond pipeline
        diamond_performance_parameters : str
            Extra arguments to provide to DIAMOND
        Returns
        -------
        String path to amino acid fasta file of reads that hit
        '''
        # Define outputs
        if search_method == 'hmmsearch':
            output_search_file = files.hmmsearch_output_path(base)
        elif search_method == 'diamond':
            output_search_file = files.diamond_search_output_basename(base)
        hit_reads_fasta = files.fa_output_path(base)
        hit_reads_orfs_fasta = files.orf_fasta_output_path(base)

        return self.search_and_extract_orfs_matching_protein_database(\
                                                    unpack,
                                                    search_method,
                                                    maximum_range,
                                                    threads,
                                                    evalue,
                                                    min_orf_length,
                                                    restrict_read_length,
                                                    diamond_database,
                                                    output_search_file,
                                                    hit_reads_fasta,
                                                    hit_reads_orfs_fasta,
                                                    diamond_performance_parameters)

    def search_and_extract_orfs_matching_protein_database(self,
                                                      unpack,
                                                      search_method,
                                                      maximum_range,
                                                      threads,
                                                      evalue,
                                                      min_orf_length,
                                                      restrict_read_length,
                                                      diamond_database,
                                                      output_search_file,
                                                      hit_reads_fasta,
                                                      hit_reads_orfs_fasta,
                                                      diamond_performance_parameters):
        '''As per aa_db_search() except slightly lower level. Search an
        input read set (unpack) and then extract the proteins that hit together
        with their containing nucleotide sequences.

        Parameters
        ----------
        output_search_file: str
            path to hmmsearch output table or diamond basename
        hit_reads_fasta: str
            path to nucleotide sequences containing hit proteins
        hit_reads_orfs_fasta: str
            path to hit proteins, unaligned

        Returns
        -------
        direction_information: dict

            {read_1: False
             ...
             read n: True}

            where True = Forward direction
            and False = Reverse direction

        result: DBSearchResult object containing file locations and hit
        information

        '''

        orfm = OrfM(min_orf_length=min_orf_length,
                     restrict_read_length=restrict_read_length)
        extracting_orfm = OrfM(min_orf_length=min_orf_length,
                      restrict_read_length=restrict_read_length)

        if search_method == 'hmmsearch':
            # run hmmsearch
            search_result = self.hmmsearch(
                                           output_search_file,
                                           unpack.get_file_as_process(),
                                           unpack,
                                           unpack.sequence_type(),
                                           threads,
                                           evalue,
                                           orfm
                                           )

        elif search_method == 'diamond':
            # run diamond
            search_result = Diamond(
                                     database=diamond_database,
                                     threads=threads,
                                     evalue=evalue,
                                     ).run(
                                           unpack.get_file_as_process(),
                                           unpack.sequence_type(),
                                           daa_file_basename=output_search_file,
                                           extra_args=diamond_performance_parameters,
                                           )
            search_result = [search_result]

        else:  # if the search_method isn't recognised
            raise Exception("Programming error: unexpected search_method %s" % search_method)

        orfm_regex = OrfM.regular_expression()

        hits = self._get_sequence_directions(search_result)

        orf_hit_readnames = list(hits.keys()) # Orf read hit names
        if unpack.sequence_type() == 'nucleotide':
            hits={(orfm_regex.match(key).groups(0)[0] if orfm_regex.match(key) else key): item for key, item in iter(hits.items())}
            hit_readnames = list(hits.keys()) # Store raw read hit names
        else:
            hit_readnames=orf_hit_readnames

        hit_reads_fasta, direction_information = self._extract_from_raw_reads(
                                                       hit_reads_fasta,
                                                       hit_readnames,
                                                       unpack.get_file_as_process(),
                                                       unpack.format(),
                                                       hits
                                                       )


        if not hit_readnames:
            hit_read_counts = [0, len(hit_readnames)]
            result = DBSearchResult(None,
                                    search_result,
                                    hit_read_counts,
                                    None)

            return result, direction_information


        if unpack.sequence_type() == 'nucleotide':
            # Extract the orfs of these reads that hit the original search
            self._extract_orfs(
                               hit_reads_fasta,
                               extracting_orfm,
                               orf_hit_readnames,
                               hit_reads_orfs_fasta,
                               search_method,
                               list(search_result[0].each([SequenceSearchResult.QUERY_ID_FIELD,
                                                           SequenceSearchResult.ALIGNMENT_DIRECTION,
                                                           SequenceSearchResult.QUERY_FROM_FIELD,
                                                           SequenceSearchResult.QUERY_TO_FIELD])
                                    )
                               )

            hit_reads_fasta = hit_reads_orfs_fasta
        slash_endings=self._check_for_slash_endings(hit_readnames)
        result = DBSearchResult(hit_reads_fasta,
                                search_result,
                                [0, len([itertools.chain(*list(hits.values()))])],  # array of hits [euk hits, true hits]. Euk hits alway 0 unless searching from 16S
                                slash_endings)  # Any reads that end in /1 or /2

        if maximum_range:
            n_hits = sum([len(x["strand"]) for x in hits.values()])
        else:
            n_hits = len(hits.keys())
        logging.info("%s read(s) detected" % n_hits)

        return result, direction_information


    @T.timeit
    def nt_db_search(self, files, base, unpack, euk_check,
                     search_method, maximum_range, threads, evalue):
        '''
        Nucleotide database search pipeline - pipeline where reads are searched
        as nucleotides, and hits are identified using nhmmer searches

        Parameters
        ----------
        files : obj
            graftm_output_paths object.
        base : str
            The name of the input file, stripped of all suffixes, and paths.
            Used for creating file names with 'files' object.
        unpack : obj
            UnpackRawReads object, returns string command that will output
            sequences to stdout when called on command line
            (use: unpack.command_line())
        euk_check : bool
            True indicates the sample will be checked for eukaryotic reads,
            False indicates not.
        search_method : str
            The method for searching e.g. 'hmmsearch' or 'diamond'
        maximum_range : int
            Maximum range that a gene can extend within a contig. Any hits
            that extend beyond this length cannot be linked. max_range is defined
            as 1.5 X the average length of all full length genes used in the
            search database. This is defined in the CONTENTS.json file within a
            gpkg.
        threads : str
            Number of threads for hmmer to use
        evalue : str
            Evalue cutoff for hmmer to use
        Returns
        -------
        String path to amino acid fasta file of reads that hit
        '''

        # Define outputs
        hmmsearch_output_table = files.hmmsearch_output_path(base)
        hit_reads_fasta = files.fa_output_path(base)
        return \
            self.search_and_extract_nucleotides_matching_nucleotide_database(\
                                                      unpack,
                                                      euk_check,
                                                      search_method,
                                                      maximum_range,
                                                      threads,
                                                      evalue,
                                                      hmmsearch_output_table,
                                                      hit_reads_fasta)

    def search_and_extract_nucleotides_matching_nucleotide_database(self,
                                                      unpack,
                                                      euk_check,
                                                      search_method,
                                                      maximum_range,
                                                      threads,
                                                      evalue,
                                                      hmmsearch_output_table,
                                                      hit_reads_fasta):
        '''As per nt_db_search() except slightly lower level. Search an
        input read set (unpack) and then extract the sequences that hit.

        Parameters
        ----------
        hmmsearch_output_table: str
            path to hmmsearch output table
        hit_reads_fasta: str
            path to hit nucleotide sequences

        Returns
        -------
        direction_information: dict

            {read_1: False
             ...
             read n: True}

            where True = Forward direction
            and False = Reverse direction

        result: DBSearchResult object containing file locations and hit
        information
        '''

        if search_method == "hmmsearch":
            # First search the reads using the HMM
            search_result, table_list = self.nhmmer(
                                                    hmmsearch_output_table,
                                                    unpack,
                                                    threads,
                                                    evalue
                                                    )


        elif search_method == 'diamond':
            raise Exception("Diamond searches not supported for nucelotide databases yet")


        if maximum_range:

            hits = self._get_read_names(
                                        search_result,  # define the span of hits
                                        maximum_range
                                        )
        else:
            hits = self._get_sequence_directions(search_result)

        hit_readnames = list(hits.keys())

        if euk_check:
            euk_reads = self._check_euk_contamination(table_list)
            hit_readnames = set([read for read in hit_readnames if read not in euk_reads])
            hits = {key:item for key, item in  iter(hits.items()) if key in hit_readnames}
            hit_read_count = [len(euk_reads), len(hit_readnames)]
        else:
            hit_read_count = [0, len(hit_readnames)]

        hit_reads_fasta, direction_information = self._extract_from_raw_reads(
                                                       hit_reads_fasta,
                                                       hit_readnames,
                                                       unpack.get_file_as_process(),
                                                       unpack.format(),
                                                       hits
                                                       )

        if not hit_readnames:
            result = DBSearchResult(None,
                                  search_result,
                                  hit_read_count,
                                  None)
        else:
            slash_endings=self._check_for_slash_endings(hit_readnames)
            result = DBSearchResult(hit_reads_fasta,
                                    search_result,
                                    hit_read_count,
                                    slash_endings)

        if maximum_range:
            n_hits = sum([len(x["strand"]) for x in hits.values()])
        else:
            n_hits = len(hits)
        logging.info("%s read(s) detected" % n_hits)

        return result, direction_information

    @T.timeit
    def align(self, input_path, output_path, directions, pipeline,
              filter_minimum):
        '''align - Takes input path to fasta of unaligned reads, aligns them to
        a HMM, and returns the aligned reads in the output path

        Parameters
        ----------
        input_path : str
        output_path : str
        reverse_direction : dict
            A dictionary of read names, with the entries being the complement
            strand of the read (True = forward, False = reverse)
        pipeline : str
            Either "P" or "D" corresponding to the protein and nucleotide (DNA)
            pipelines, respectively.


        Returns
        -------
        N/A - output alignment path known.
        '''

        # HMMalign the forward reads, and reverse complement reads.
        with tempfile.NamedTemporaryFile(prefix='for_conv_file', suffix='.fa') as fwd_fh:
            fwd_conv_file = fwd_fh.name
            with tempfile.NamedTemporaryFile(prefix='rev_conv_file', suffix='.fa') as rev_fh:
                rev_conv_file = rev_fh.name
                alignments = self._hmmalign(
                    input_path,
                    directions,
                    pipeline,
                    fwd_conv_file,
                    rev_conv_file)
                alignment_result = self.alignment_correcter(alignments,
                                                            output_path,
                                                            filter_minimum)
                return alignment_result
import os
import shutil
import tempfile
import logging
import subprocess
import extern
import json
import signal
import re
import itertools
from dendropy import Tree
from Bio import SeqIO
from io import StringIO
from graftm.sequence_searcher import SequenceSearcher
from graftm.dendropy_tree_cleaner import DendropyTreeCleaner
from graftm.taxonomy_extractor import TaxonomyExtractor
from graftm.getaxnseq import Getaxnseq
from graftm.deduplicator import Deduplicator
from graftm.sequence_io import SequenceIO
from graftm.graftm_package import GraftMPackageVersion3, GraftMPackage
from graftm.decorator import Decorator
from graftm.greengenes_taxonomy import GreenGenesTaxonomy
from graftm.sequence_searcher import SequenceSearcher
            search_re = re.compile(r'([^ACDEFGHIKLMNPQRSTVWY\-X])',re.IGNORECASE)
            replace_char = 'X'
        elif package_type == Create._NUCLEOTIDE_PACKAGE_TYPE:
            search_re = re.compile(r'([^ATGC\-N])',re.IGNORECASE)
            replace_char = 'N'

        for s in sequences:
            if package_type == Create._NUCLEOTIDE_PACKAGE_TYPE:
                s.seq = s.seq.replace('U','T')
                s.seq = s.seq.replace('u','t')
            newseq = s.seq
            m = search_re.search(newseq)
            while m:
                logging.warning(
                    "Found a non-standard character in the "
                    "sequence of %s: e.g. '%s'" % (
                        s.name, m.group()))
                newseq = newseq[:m.start()] + replace_char + newseq[(m.start()+1):]
                m = search_re.search(newseq)
            s.seq = newseq

    def _check_for_duplicate_sequence_names(self, fasta_file_path):
        """Test if the given fasta file contains sequences with duplicate
        sequence names.

        Parameters
        ----------
        fasta_file_path: string
            path to file that is to be checked

        Returns
        -------
        The name of the first duplicate sequence found, else False.

        """
        found_sequence_names = set()
        for record in SeqIO.parse(fasta_file_path, 'fasta'):
            name = record.name
            if name in found_sequence_names:
                return name
            found_sequence_names.add(name)
        return False

    def _test_package(self, package_path):
        '''Give a GraftM package a spin, and see if it works in reality with default
        parameters (i.e. pplacer). If it does not work, then raise an error.

        Parameters
        ----------
        package_path: str
            path to graftm_package to be tested
        '''
        pkg = GraftMPackage.acquire(package_path)
        with tempfile.TemporaryDirectory() as graftM_graft_test_dir_name:
            # Take a subset of sequences for testing
            with tempfile.NamedTemporaryFile(suffix=".fa",mode='w') as tf:
                seqio = SequenceIO()
                with open(pkg.unaligned_sequence_database_path()) as f:
                    seqio.write_fasta(
                        itertools.islice(seqio.each_sequence(f), 10),
                        tf)
                tf.flush()
                cmd = "graftM graft --forward %s --graftm_package %s --output_directory %s --force" %(
                    tf.name, package_path, graftM_graft_test_dir_name)
                extern.run(cmd)


    def main(self, **kwargs):
        alignment = kwargs.pop('alignment',None)
        sequences = kwargs.pop('sequences',None)
        taxonomy = kwargs.pop('taxonomy',None)
        rerooted_tree = kwargs.pop('rerooted_tree',None)
        unrooted_tree = kwargs.pop('unrooted_tree',None)
        tree_log = kwargs.pop('tree_log', None)
        no_tree = kwargs.pop('no_tree', False)
        prefix = kwargs.pop('prefix', None)
        rerooted_annotated_tree = kwargs.pop('rerooted_annotated_tree', None)
        user_hmm = kwargs.pop('hmm', None)
        search_hmm_files = kwargs.pop('search_hmm_files',None)
        min_aligned_percent = kwargs.pop('min_aligned_percent',0.01)
        taxtastic_taxonomy = kwargs.pop('taxtastic_taxonomy', None)
        taxtastic_seqinfo = kwargs.pop('taxtastic_seqinfo', None)
        force_overwrite = kwargs.pop('force',False)
        graftm_package = kwargs.pop('graftm_package',False)
        dereplication_level = kwargs.pop('dereplication_level',False)
        threads = kwargs.pop('threads',5)

        if len(kwargs) > 0:
            raise Exception("Unexpected arguments detected: %s" % kwargs)
        seqio = SequenceIO()
        locus_name = (os.path.basename(sequences).split('.')[0]
                      if sequences
                      else os.path.basename(alignment).split('.')[0])
        tmp = tempfile.TemporaryDirectory()
        base = os.path.join(tmp.name, locus_name)
        insufficiently_aligned_sequences = [None]
        removed_sequence_names = []
        tempfiles_to_close = []

        if prefix:
            output_gpkg_path = prefix
        else:
            output_gpkg_path = "%s.gpkg" % locus_name

        if os.path.exists(output_gpkg_path):
            if force_overwrite:
                logging.warning("Deleting previous directory %s" % output_gpkg_path)
                shutil.rmtree(output_gpkg_path)
            else:
                raise Exception("Cowardly refusing to overwrite gpkg to already existing %s" % output_gpkg_path)
        logging.info("Building gpkg for %s" % output_gpkg_path)

        # Read in taxonomy somehow
        gtns = Getaxnseq()
        if rerooted_annotated_tree:
            logging.info("Building seqinfo and taxonomy file from input annotated tree")
            taxonomy_definition = TaxonomyExtractor().taxonomy_from_annotated_tree(\
                    Tree.get(path=rerooted_annotated_tree, schema='newick'))
        elif taxonomy:
            logging.info("Building seqinfo and taxonomy file from input taxonomy")

            taxonomy_definition = GreenGenesTaxonomy.read_file(taxonomy).taxonomy