class SequenceSearchResult:
    QUERY_FROM_FIELD = 'query_from'
    QUERY_TO_FIELD = 'query_to'
    QUERY_LENGTH_FIELD = 'query_length'
    HIT_FROM_FIELD = 'hit_from'
    HIT_TO_FIELD = 'hit_to'
    ALIGNMENT_LENGTH_FIELD = 'alignment_length'
    ALIGNMENT_BIT_SCORE = 'alignment_bit_score'
    ALIGNMENT_DIRECTION = 'alignment_direction'
    HIT_ID_FIELD = 'hit_id'
    QUERY_ID_FIELD = 'query_id'
    HMM_NAME_FIELD = 'hmm_name'
    ACCESSION_ID_FIELD = 'accession_id'
    PERCENT_ID_FIELD = 'percent_id'
    MISMATCH_FIELD = "mismatch"
    EVALUE_FIELD = "evalue"



    def __init__(self):
        self.fields = []
        self.results = []

    def each(self, field_names):
        """Iterate over the results, yielding a list for each result, where
        each element corresponds to the field given in the field_name parameters

        Parameters
        ----------
        field_names: list of str
            The names of the fields to be returned during iteration

        Returns
        -------
        None

        Exceptions
        ----------
        raises something when a field name is not in self.fields
        """
        field_ids = []
        for f in field_names:
            # below raises error if the field name is not found, so
            # don't need to account for that.
            field_ids.append(self.fields.index(f))

        for r in self.results:
            yield([r[i] for i in field_ids]) class GraftMFiles:
    
    def __init__(self, old_title, outdir, direction):
        if direction in ['forward', 'reverse', 'interleaved']:
            self.basename = os.path.join(direction, old_title + '_' + direction)
        elif direction == False:
            self.basename = old_title
        else:
            raise Exception('Programming Error.')   
         
        self.outdir = outdir
    
    def search_otu_table(self):
        return os.path.join(self.outdir, "search_otu_table.txt")
    
    def hmmsearch_output_path(self, out_path):
        return os.path.join(self.outdir, out_path, "%s.hmmout.txt" % self.basename)
    
    def diamond_search_output_basename(self, out_path):
        '''Does not include the .daa part that diamond creates'''
        return os.path.join(self.outdir, out_path, "%s_diamond_search" % self.basename)
    
    def diamond_assignment_output_basename(self, out_path):
        '''Does not include the .daa part that diamond creates'''
        return os.path.join(self.outdir, out_path, "%s_diamond_assignment" % out_path)
    
    def read_tax_output_path(self, out_path):
        return os.path.join(self.outdir, out_path, "%s_read_tax.tsv" % out_path)

    def jplace_output_path(self):
        return "placements.jplace"
    
    def euk_free_path(self, out_path):
        return os.path.join(self.outdir, out_path, "%s_euk_free.fa" % self.basename)
    
    def euk_contam_path(self, out_path):
        return os.path.join(self.outdir, out_path, "%s_euk_contam.txt" % self.basename)
    
    def summary_table_output_path(self, out_path):
        return os.path.join(self.outdir, out_path, "%s_count_table.txt" % self.basename)
    
    def krona_output_path(self):
        return os.path.join(self.outdir, "krona.html")
    
    def aligned_fasta_output_path(self, out_path):
        return os.path.join(self.outdir, out_path, "%s_hits.aln.fa" % self.basename)

    def orf_output_path(self, out_path):
        return os.path.join(self.outdir, out_path, "%s_orf" % self.basename)

    def orf_titles_output_path(self, out_path):
        return os.path.join(self.outdir, out_path, "%s_orf.titles" % self.basename)

    def orf_fasta_output_path(self, out_path):
        return os.path.join(self.outdir, out_path, "%s_orf.fa" % self.basename)

    def conv_output_for_path(self, out_path):
        return os.path.join(self.outdir, out_path, "%s_conv_for.faa" % self.basename)
    
    def output_for_path(self, out_path):
        return os.path.join(self.outdir, out_path, "%s_for.faa" % self.basename)
    
    def output_rev_path(self, out_path):
        return os.path.join(self.outdir, out_path, "%s_rev.faa" % self.basename)
    
    def conv_output_rev_path(self, out_path):
        return os.path.join(self.outdir, out_path, "%s_conv_rev.faa" % self.basename)
    
    def comb_aln_fa(self):
        return os.path.join(self.outdir, "combined_alignment.aln.fa")

    def fa_output_path(self, out_path):
        return os.path.join(self.outdir, out_path, "%s_hits.fa" % self.basename)     
        
    def readnames_output_path(self, out_path):
        return os.path.join(self.outdir, out_path, "%s_readnames.txt" % self.basename)
    
    def sto_output_path(self, out_path):
        return os.path.join(self.outdir, out_path, "%s.sto" % self.basename)

    def basic_stats_path(self):
        return os.path.join(self.outdir, "basic_stats.txt")

    def for_aln_path(self, out_path):
        return os.path.join(self.outdir, out_path, "%s_for_aln.fa" % self.basename)
        
    def rev_aln_path(self, out_path):
        return os.path.join(self.outdir, out_path, "%s_rev_aln.fa" % self.basename)
    
    def combined_biom_output_path(self):
        return os.path.join(self.outdir, "graftm.biom")
    
    def combined_summary_table_output_path(self):
        return os.path.join(self.outdir, "combined_count_table.txt")
    
    def expand_search_hmm_path(self):
        return os.path.join(self.outdir, "expand_search.hmm")
      
    def base(self, out_path):
        return os.path.join(self.outdir, out_path, "%s" % self.basename) class SearchTableWriter:
    '''
    Class for writing the search output OTU table. Basically a summary
    of hits to the HMM/Diamond searched in the following format:

             #ID    Metagenome_1    Metagenome_2    ...
            HMM1    50              6
            HMM2    195             41
            HMM3    2               20120
            ...

    You just need to specify a series of SequenceSearchResult objects, and an
    output path.
    '''

    def _interpret_hits(self, results_list, base_list):
        '''Sort reads that hit multiple HMMs to the databases to which they had
        the highest bit score. Return a dictionary containing HMMs as keys, and
        number of hits as the values.

        This function is set up so that the read names could easily be returned
        instead of numbers, for future development of GraftM

        Parameters
        ----------
        results_list: list
            Iterable if SequenceSearchResult objects. e.g.
                [SequenceSearchResult_1, SequenceSearchResult_2, ...]

        base_list: list
            Iterable of the basenames for each sequence file provided to graftM
            e.g.
                [sample_1, sample_2, ...]

        Returns
        -------
        dictionary:
            Contains samples as entries. The value for each sample is another
            dictionary with HMM as the key, and number of hits as values:
                {"sample_1":{HMM_1: 12
                             HMM_2: 35
                             HMM_3: 1258
                             ...}
                 ...
                }

        '''
        logging.debug("Sorting reads into HMMs by bit score")

        run_results = {}

        ########################################################################
        ################## - Sort reads to best hit db - #######################
        for base, results in zip(base_list, results_list): # For each sample
            search_results = {}
            for search in results():
                search_list = list(
                                   search.each([SequenceSearchResult.QUERY_ID_FIELD,
                                                SequenceSearchResult.ALIGNMENT_BIT_SCORE,
                                                SequenceSearchResult.HMM_NAME_FIELD])
                                   )
                for hit in search_list:
                    if hit[0] in search_results:
                        if float(hit[1]) > search_results[hit[0]][0]:
                            search_results[hit[0]] = [float(hit[1]), hit[2]]
                    else:
                        search_results[hit[0]] = [float(hit[1]), hit[2]]
            run_results[base] = search_results

        ########################################################################
        ################## - Gather counts for each db - #######################
        db_count = {}
        for run in run_results.keys():
            run_count = {}
            for entry in list(run_results[run].values()):
                key = entry[1]
                if key in run_count:
                    run_count[key] += 1
                else:
                    run_count[key] = 1
            db_count[run] = run_count

        return db_count

    def _write_results(self, db_count, output_path):
        '''Write the table to the output_path directory

        db_count: dict
            Contains samples as entries. The value for each sample is another
            dictionary with HMM as the key, and number of hits as values:
                {"sample_1":{HMM_1: 12
                             HMM_2: 35
                             HMM_3: 1258
                             ...}
                 ...
                }

        output_path: str
            Path to output file to which the resultant output file will be
            written to.
        '''

        logging.debug("Writing search otu table to file: %s" % output_path)

        output_dict = {}

        for idx, value_dict in enumerate(db_count.values()):
            for database, count in value_dict.items():
                if database in output_dict:
                    output_dict[database].append(str(count))
                else:
                    output_dict[database] = ['0']*idx + [str(count)]

            for key, item in output_dict.items():
                if len(item) == idx:
                    output_dict[key].append('0')

        with open(output_path, 'w') as out:
            out.write('\t'.join(["#ID"] + list(db_count.keys())) + '\n')
            for key, item in output_dict.items():
                out.write("%s\t%s" % (key, '\t'.join(item)) + '\n' )

    def build_search_otu_table(self, search_results_list, base_list, output_path):
        '''
        Build an OTU from SequenceSearchResult objects

        Parameters
        ----------
        search_results_list: list
            Iterable if SequenceSearchResult objects. e.g.
                [SequenceSearchResult_1, SequenceSearchResult_2, ...]
        base_list: list
            Iterable of the basenames for each sequence file provided to graftM
            e.g.
                [sample_1, sample_2, ...]
        output_path: str
            Path to output file to which the resultant output file will be
            written to.
        '''

        db_count = self._interpret_hits(search_results_list,
                                        base_list)

        self._write_results(db_count, output_path) 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 NoInputSequencesException(Exception):
    def __init__(self, command):
        """Instantiate with the command used that went amiss"""
        self.command = command class HouseKeeping:
    ### Functions for setting up the graftM pipeline to run correctly, and
    ### general housekeeping things like adding and removing directories and
    ### files.

    HMMSEARCH_AND_DIAMOND_SEARCH_METHOD = 'hmmsearch+diamond'
    
    DIAMOND_SEARCH_METHOD   = 'diamond'
    HMMSEARCH_SEARCH_METHOD = 'hmmsearch'


    def file_basename(self, file):
        '''
        Strips the path and last extension from the file variable. If the
        extension is found to be valid, the basename will be returned. Otherwise
        an error will be raise and graftM will exit
        '''
        valid_extensions = set(".tree",
                               ".tre")
        split_file = os.path.basename(file).split('.')
        base, suffix = '.'.join(split_file[:-1]), split_file[-1]

        if suffix in valid_extensions:
            return base
        else:
            logging.error("Invalid file extension found on file: %s" % file)
            logging.error("For trees, please provide a file with one of the \
following extensions: %s" % ' '.join(valid_extensions.keys()))
            raise InvalidFileExtensionError


    def set_euk_hmm(self, args):
        'Set the hmm used by graftM to cross check for euks.'
        if hasattr(args, 'euk_hmm_file'):
            pass
        elif not hasattr(args, 'euk_hmm_file'):
            # set to path based on the location of bin/graftM, which has
            # a more stable relative path to the HMM when installed through
            # pip.
            setattr(args, 'euk_hmm_file', os.path.join(os.path.dirname(inspect.stack()[-1][1]),'..','share', '18S.hmm'))
        else:
            raise Exception('Programming Error: setting the euk HMM')

    def setpipe(self, hmm):
        ## Read in the hmm file and return the evalue trusted cutoff and the
        ## type of HMM
        hmm_tc = None
        hmm_read = [x.rstrip() for x in open(hmm, 'r').readlines() if x.startswith('ALPH') or x.startswith('TC')]
        if len(hmm_read) > 2:
            raise Exception('Possibly Programming Error: Misread HMM library.')
        for item in hmm_read:
            if item.startswith('ALPH'):
                if item.split()[1] == 'DNA' or item.split()[1] == 'RNA':
                    hmm_type = 'D'
                elif item.split()[1] == 'amino':
                    hmm_type = 'P'
            elif item.startswith('TC'):
                hmm_tc = item.split()[1]
            else:
                raise Exception('Possibly Programming Error: Misread HMM library.')
        return hmm_type, hmm_tc

    def delete(self, delete_list):
        for item in delete_list:
            try:
                os.remove(item)
            except:
                pass

    def make_working_directory(self, directory_path, force):
        if force:
            shutil.rmtree(directory_path, ignore_errors=True)
            os.mkdir(directory_path)
        else:
            try:
                os.mkdir(directory_path)
            except:
                logging.error('Directory %s already exists. Exiting to prevent over-writing' % directory_path)
                raise Exception('Directory %s already exists. Exiting to prevent over-writing'% directory_path)

    def parameter_checks(self, args):
        ## Check that the necessary files are in place
        if args.subparser_name == 'graft':
            # Check that the placement cutoff is between 0.5 and 1
            if float(args.placements_cutoff) < float(0.5) or float(args.placements_cutoff) > float(1.0):
                logging.info('Please specify a confidence level (-d) between 0.5 and 1.0! Found: %s' % args.placements_cutoff)
                exit(1)

            if args.interleaved and args.forward:
                logging.info('Please specify either reads with one of'
                             '--forward or --interleaved, not both')
                exit(1)
            if (not args.interleaved) and (not args.forward):
                logging.info('Please specify either reads with either'
                             '--forward or --interleaved')
                exit(1)
            if args.interleaved:
                self._check_file_existence(args.interleaved)
            if args.forward:
                self._check_file_existence(args.forward)
            if args.reverse:
                self._check_file_existence(args.reverse)


            # Determine the File format based on the suffix

            sequence_file_list = []
            if args.reverse:
                if len(args.forward) != len(args.reverse):
                    logging.error('Confusing input. There appears to be different numbers of forward and reverse files specified')
                for i, forward_file in enumerate(args.forward):
                    sequence_file_list.append([forward_file, args.reverse[i]])
            elif args.forward:
                sequence_file_list = [[f] for f in args.forward]
            if args.interleaved:
                sequence_file_list.extend([[f, None] for f in args.interleaved])

            return sequence_file_list

    def _check_file_existence(self, files):
        '''Iterate through files and exit(1) if any do not pass
        os.path.isfile'''
        for f in files:
            if not os.path.isfile(f):
                logging.error("The file '%s' does not appear to exist, stopping" % f)
                exit(1)



    def get_maximum_range(self, hmm):
        '''
        If no maximum range has been specified, and if using a hmm search, a
        maximum range can be determined by using the length of the HMM

        Parameters
        ----------
        hmm : str
            path to hmm profile

        Returns
        -------
        Length to search to when linking hits on a single contig
        '''
        length=int([x for x in open(hmm) if x.startswith("LENG")][0].split()[1])
        max_length=round(length*1.5, 0)

        return max_length

    def set_attributes(self, args):


        # Read graftM package and assign HMM and refpkg file
        if args.no_merge_reads:
            setattr(args, 'merge_reads', False)
        else:
            if args.reverse or args.interleaved:
                setattr(args, 'merge_reads', True)
            else:
                setattr(args, 'merge_reads', False)
                
        if args.graftm_package:
            if not os.path.isdir(args.graftm_package):
                raise Exception("%s does not exist. Are you sure you provided the correct path?" % args.graftm_package)
            else:
                gpkg = GraftMPackage.acquire(args.graftm_package)
                if hasattr(args, 'search_hmm_files'): # If a hmm is specified, overwrite the one graftM package
                    setattr(args, 'aln_hmm_file', gpkg.alignment_hmm_path())
                    setattr(args, 'reference_package', gpkg.reference_package_path())
                else:
                    setattr(args, 'search_hmm_files', [])
                    for hmm in gpkg.search_hmm_paths():
                        args.search_hmm_files.append(hmm)
                    setattr(args, 'aln_hmm_file', gpkg.alignment_hmm_path())
                    setattr(args, 'reference_package', gpkg.reference_package_path())

        elif hasattr(args, 'search_diamond_files'):
            if args.search_method == self.DIAMOND_SEARCH_METHOD:
                if hasattr(args, 'aln_hmm_file'):
                    pass
                else:
                    raise Exception("aln_hmm_file not specified")
            else:
                raise Exception("Specified DIAMOND databases when not using the diamond search pipeline. Using: %s" % (args.search_method))

        elif hasattr(args, 'search_hmm_files'):
            if args.search_method == self.HMMSEARCH_SEARCH_METHOD:
                if not hasattr(args, 'aln_hmm_file'):
                    if len(args.search_hmm_files) == 1:
                        if not args.search_only:
                            setattr(args, 'aln_hmm_file', args.search_hmm_files[0])
                    else:
                        raise Exception("Multiple search HMMs specified, but aln_hmm_file not specified")

            else:
                raise Exception("Specified HMM search_hmm_files when not using the hmmsearch pipeline. Using: %s" % (args.search_method))

        elif hasattr(args, 'search_hmm_list_file'):
            if args.search_method == self.HMMSEARCH_SEARCH_METHOD:
                setattr(args, 'search_hmm_files', [x.rstrip() for x in open(args.search_hmm_list_file).readlines()])
                if not hasattr(args, 'aln_hmm_file'):
                    if not args.search_only:
                        raise Exception("Multiple search HMMs specified, but aln_hmm_file not specified")
            else:
                raise Exception("Specified HMM search_hmm_files when not using the hmmsearch pipeline. Using: %s" % (args.search_method))

        else:
            if args.search_only:
                if args.search_diamond_file:
                    args.search_method = self.DIAMOND_SEARCH_METHOD
                    args.search_hmm_files = None
            else:
                raise Exception('No gpkg, HMM, or DIAMOND database was specified, so there is no reference database to search with.') class Stats_And_Summary:

    def __init__(self): pass

    def coverage_of_hmm(self, hmm, count_table, coverage_table, avg_read_length):

        for line in open(hmm):

            if line.startswith('LENG'):
                length = line.split()[1]

        with open(coverage_table, 'w') as ct:
            write = ['#ID     20110816_S1D    ConsensusLineage\n']
            for line in open(count_table):

                    if line.startswith('#'):
                        continue

                    splt = line.split()

                    cov = str(round((float(splt[1])*float(avg_read_length)) / float(length), 3))

                    write.append("%s\t%s\t%s\n" % (
                                    splt[0],
                                    cov,
                                    splt[2]))
            for entry in write:
                ct.write(entry)

    def readTax(self, placements, output_path):
        with open(output_path, 'w') as out:
            for read, tax in placements.items():
                out.write("%s\t%s\n" % (read, '; '.join(tax)))

    def build_basic_statistics(self, times, hit_read_count_list, placed_reads, base_list, output):

        output_lines = ["Basic run statistics (count):"]
        output_lines.append("Files:\t%s" % '\t'.join(base_list))
        if any([x[0] for x in hit_read_count_list if x[0] > 0]):
            output_lines.append("18S reads filtered:\t%s" % '\t'.join([str(x[0]) for x in hit_read_count_list]))
        output_lines.append("reads detected:\t%s" % '\t'.join([str(x[1]) for x in hit_read_count_list]))
        output_lines.append("reads placed in tree:\t%s" % '\t'.join([str(x) for x in placed_reads]))
        output_lines.append("Runtime (seconds):")
        output_lines.append("Search step:\t%s" % times[0])
        output_lines.append("Alignment step:\t%s" % times[1])
        output_lines.append("Tree insertion step:\t%s" % times[2])


        with open(output, 'w') as stats_file:
            for line in output_lines:
                stats_file.write(line + '\n')

    def _iterate_otu_table_rows(self, read_taxonomies):
        '''yield that which is required for an OTU table: taxonomy, and
        count of that taxonomy in each sample as an array

        Parameters
        ---------
        read_taxonomies:
            a list of hashes, where the position in the list corresponds to the
            sample list, the key is the read name, and the value is an array
            of taxonomic info

        Yield
        -----
        OTU_ID: int, made up by this function
        taxonomy: array where each element is a taxonomic level
        count: number of observations of that taxonomy in each sample respectively

        Return
        ------
        Nothing, use this as an iterator'''

        taxonomy_string_to_taxonomy_array = {}
        taxonomy_string_to_counts = {}
        sample_index = 0
        num_samples = len(read_taxonomies)
        for read_to_taxonomy in read_taxonomies: # For each sample
            for taxonomy_array in read_to_taxonomy.values(): # For each read
                taxonomy_string = '; '.join(taxonomy_array)
                if taxonomy_string in taxonomy_string_to_taxonomy_array:
                    if taxonomy_string_to_taxonomy_array[taxonomy_string] != taxonomy_array:
                        raise Exception("Programming error: two different taxonomies had same taxonomy string")
                    try:
                        taxonomy_string_to_counts[taxonomy_string][sample_index] += 1
                    except KeyError:
                        taxonomy_string_to_counts[taxonomy_string][sample_index] = 1
                else:
                    taxonomy_string_to_taxonomy_array[taxonomy_string] = taxonomy_array
                    taxonomy_string_to_counts[taxonomy_string] = [0]*num_samples
                    taxonomy_string_to_counts[taxonomy_string][sample_index] = 1
            sample_index += 1

        otu_id = 1
        for tax_string, counts_array in taxonomy_string_to_counts.items():
            array = []
            for i in range(num_samples):
                try:
                    array.append(counts_array[i])
                except IndexError:
                    array.append(0)
            yield otu_id,\
                taxonomy_string_to_taxonomy_array[tax_string],\
                array
            otu_id += 1

    def write_biom(self, sample_names, read_taxonomies, biom_file_io):
        '''Write the OTU info to a biom IO output stream

        Parameters
        ----------
        sample_names: String
            names of each sample (sample_ids for biom)
        read_taxonomies: Array of hashes as per _iterate_otu_table_rows()
        biom_file_io: io
            open writeable stream to write biom contents to

        Returns True if successful, else False'''
        counts = []
        observ_metadata = []
        otu_ids = []
        for otu_id, tax, count in self._iterate_otu_table_rows(read_taxonomies):
            if len(count) != len(sample_names):
                raise Exception("Programming error: mismatched sample names and counts")
            counts.append(count)
            observ_metadata.append({'taxonomy': tax})
            otu_ids.append(str(otu_id))
        if len(counts) == 0:
            logging.info("Not writing BIOM file since no sequences were assigned taxonomy")
            return True
        table = Table(np.array(counts),
                      otu_ids, sample_names, observ_metadata,
                      [{}]*len(sample_names), table_id='GraftM Taxonomy Count Table')
        try:
            table.to_hdf5(biom_file_io, 'GraftM graft')
            return True
        except RuntimeError as e:
            logging.warn("Error writing BIOM output, file not written. The specific error was: %s" % e)
            return False

    def write_tabular_otu_table(self, sample_names, read_taxonomies, combined_output_otu_table_io):
        '''A function that takes a hash of trusted placements, and compiles them
        into an OTU-esque table.'''
        delim = '\t'
        combined_output_otu_table_io.write(delim.join(['#ID',
                                                       delim.join(sample_names),
                                                       'ConsensusLineage']))
        combined_output_otu_table_io.write("\n")
        for otu_id, tax, counts in self._iterate_otu_table_rows(read_taxonomies):
            combined_output_otu_table_io.write(delim.join(\
                (str(otu_id),
                 delim.join([str(c) for c in counts]),
                 '; '.join(tax)))+"\n")

    def write_krona_plot(self, sample_names, read_taxonomies, output_krona_filename):
        '''Creates krona plot at the given location. Assumes the krona executable
        ktImportText is available on the shell PATH'''
        tempfiles = []
        for n in sample_names:
            tempfiles.append(
                tempfile.NamedTemporaryFile(
                    prefix='GraftMkronaInput',
                    suffix=n,
                    mode='w'))

        delim='\t'
        for _, tax, counts in self._iterate_otu_table_rows(read_taxonomies):
            for i, c in enumerate(counts):
                if c != 0:
                    tempfiles[i].write(delim.join((str(c),
                                                  delim.join(tax)
                                                  ))+"\n")

        for t in tempfiles:
            t.flush()

        cmd = ["ktImportText",'-o',output_krona_filename]
        for i, tmp in enumerate(tempfiles):
            cmd.append(','.join([tmp.name,sample_names[i]]))

        # run the actual krona
        cmd = ' '.join(cmd)
        extern.run(cmd)

        # close tempfiles
        for t in tempfiles:
            t.close() class Pplacer:
    ### Contains function related to processing alignment files to jplace files
    ### and running comparisons between forward and revere reads if reverse
    ### reads are provided.

    def __init__(self, refpkg):
        self.refpkg = refpkg
        self.hk = HouseKeeping()

    # Run pplacer
    def pplacer(self, output_file, output_path, input_path, threads):
        ## Runs pplacer on concatenated alignment file
        cmd = "pplacer -j %s --verbosity 0 --out-dir %s -c %s %s" % (str(threads), output_path, self.refpkg, input_path) # Set command
        extern.run(cmd)
        output_path = '.'.join(input_path.split('.')[:-1]) + '.jplace'
        return output_path

    def alignment_merger(self, alignment_files, output_alignment_path):
        ## Concatenate aligned read_files into one file. Each read with it's
        ## own unique identifier assigning it to a particular origin file
        alias_hash = {} # Set up a hash with file names and their unique identifier
        file_number = 0 # file counter (unique identifier)
        with open(output_alignment_path, 'w') as output:
            for alignment_file in alignment_files: # For each alignment
                if alignment_file is not None:
                    alignments = list(SeqIO.parse(open(alignment_file, 'r'), 'fasta')) # read list
                    for record in alignments: # For each record in the read list
                        record.id = record.id + '_' + str(file_number) # append the unique identifier to the record id
                    SeqIO.write(alignments, output, "fasta") # And write the reads to the file
                    alias_hash[str(file_number)] = {'output_path': os.path.join(os.path.dirname(alignment_file), 'placements.jplace')}
                file_number += 1
        return alias_hash

    def convert_cluster_dict_keys_to_aliases(self, cluster_dict, alias_hash):
        '''
        Parameters
        ----------
        cluster_dict : dict
            dictionary stores information on pre-placement clustering
        alias_hash : dict
            Stores information on each input read file given to GraftM, the
            corresponding reads found within each file, and their taxonomy

        Returns
        --------
        updated cluster_dict dict containing alias indexes for keys.
        '''
        output_dict = {}
        directory_to_index_dict = {os.path.split(item["output_path"])[0] : key
                            for key, item in iter(alias_hash.items())}
        for key, item in cluster_dict.items():
            cluster_file_directory = os.path.split(key)[0]
            cluster_idx = directory_to_index_dict[cluster_file_directory]
            output_dict[cluster_idx] = item
        return output_dict



    def jplace_split(self, original_jplace, cluster_dict):
        '''
        To make GraftM more efficient, reads are dereplicated and merged into
        one file prior to placement using pplacer. This function separates the
        single jplace file produced by this process into the separate jplace
        files, one per input file (if multiple were provided) and backfills
        abundance (re-replicates?) into the placement file so analyses can be
        done using the placement files.

        Parameters
        ----------
        original_jplace : dict (json)
            json .jplace file from the pplacer step.
        cluster_dict : dict
            dictionary stores information on pre-placement clustering

        Returns
        -------
        A dict containing placement hashes to write to
        new jplace file. Each key represents a file alias
        '''
        output_hash = {}

        for placement in original_jplace['placements']: # for each placement
            alias_placements_list = []
            nm_dict = {}

            p = placement['p']
            if 'nm' in placement.keys():
                nm = placement['nm']
            elif 'n' in placement.keys():
                nm = placement['n']
            else:
                raise Exception("Unexpected jplace format: Either 'nm' or 'n' are expected as keys in placement jplace .JSON file")

            for nm_entry in nm:
                nm_list = []
                placement_read_name, plval = nm_entry
                read_alias_idx = placement_read_name.split('_')[-1] # Split the alias
                                    # index out of the read name, which
                                    # corresponds to the input file from
                                    # which the read originated.
                read_name = '_'.join(placement_read_name.split('_')[:-1])
                read_cluster = cluster_dict[read_alias_idx][read_name]
                for read in read_cluster:
                    nm_list.append([read.name, plval])
                if read_alias_idx not in nm_dict:
                    nm_dict[read_alias_idx] = nm_list
                else:
                    nm_dict[read_alias_idx] += nm_entry

            for alias_idx, nm_list in nm_dict.items():
                placement_hash = {'p': p,
                                  'nm': nm_list}
                if alias_idx not in output_hash:
                    output_hash[alias_idx] = [placement_hash]
                else:
                    output_hash[alias_idx].append(placement_hash)
        return output_hash

    def write_jplace(self, original_jplace, alias_hash):
        # Write the jplace file to their respective file paths.
        for alias_idx in alias_hash.keys():
            output = {'fields'     : original_jplace['fields'],
                      'version'    : original_jplace['version'],
                      'tree'       : original_jplace['tree'],
                      'placements' : alias_hash[alias_idx]['place'],
                      'metadata'   : original_jplace['metadata']}
            with open(alias_hash[alias_idx]['output_path'], 'w') as output_io:
                json.dump(output, output_io, ensure_ascii=False, indent=3, separators=(',', ': '))

    @T.timeit
    def place(self, reverse_pipe, seqs_list, resolve_placements, files, args,
              slash_endings, tax_descr, clusterer):
        '''
        placement - This is the placement pipeline in GraftM, in aligned reads
                    are placed into phylogenetic trees, and the results interpreted.
                    If reverse reads are used, this is where the comparisons are made
                    between placements, for the summary tables to be build in the
                    next stage.

        Parameters
        ----------
        reverse_pipe : bool
            True: reverse reads are placed separately
            False: no reverse reads to place.
        seqs_list : list
            list of paths to alignment fastas to be placed into the tree
        resolve_placements : bool
            True:resolve placements to their most trusted taxonomy
            False: classify reads to their most trusted taxonomy, until the
                   confidence cutoff is reached.
        files : list
            graftM output file name object
        args : obj
            argparse object
        Returns
        -------
        trusted_placements : dict
            dictionary of reads and their trusted placements
        '''
        trusted_placements = {}
        files_to_delete = []
        # Merge the alignments so they can all be placed at once.
        alias_hash = self.alignment_merger(seqs_list, files.comb_aln_fa())
        files_to_delete += seqs_list
        files_to_delete.append(files.comb_aln_fa())
        if os.path.getsize(files.comb_aln_fa()) == 0:
            logging.debug("Combined alignment file has 0 size, not running pplacer")
            to_return = {}
            for idx, file in enumerate(seqs_list):
                base_file=os.path.basename(file).replace('_forward_hits.aln.fa', '')
                to_return[base_file] = {}
            return to_return

        # Run pplacer on merged file
        jplace = self.pplacer(files.jplace_output_path(), args.output_directory, files.comb_aln_fa(), args.threads)
        files_to_delete.append(jplace)
        logging.info("Placements finished")

        #Read the json of refpkg
        logging.info("Reading classifications")
        classifications=Classify(tax_descr).assignPlacement(
                                                           jplace,
                                                           args.placements_cutoff,
                                                           resolve_placements
                                                           )
        logging.info("Reads classified")
        # If the reverse pipe has been specified, run the comparisons between the two pipelines. If not then just return.

        for idx, file in enumerate(seqs_list):

            if reverse_pipe:
                base_file=os.path.basename(file).replace('_forward_hits.aln.fa', '')
                forward_gup=classifications.pop(sorted(classifications.keys())[0])
                reverse_gup=classifications.pop(sorted(classifications.keys())[0])
                seqs_list.pop(idx+1)
                placements_hash = Compare().compare_placements(
                                                               forward_gup,
                                                               reverse_gup,
                                                               args.placements_cutoff,
                                                               slash_endings,
                                                               base_file
                                                               )
                trusted_placements[base_file]=placements_hash['trusted_placements']

            else: # Set the trusted placements as
                base_file=os.path.basename(file).replace('_hits.aln.fa', '')
                trusted_placements[base_file]={}
                if str(idx) in classifications:
                    for read, entry in classifications[str(idx)].items():
                        trusted_placements[base_file][read] = entry['placement']
        # Split the original jplace file
        # and write split jplaces to separate file directories
        with open(jplace) as f: jplace_json = json.load(f)
        cluster_dict = self.convert_cluster_dict_keys_to_aliases(clusterer.seq_library,
                                                                 alias_hash)
        hash_with_placements = self.jplace_split(jplace_json,
                                                 cluster_dict)

        for file_alias, placement_entries_list in hash_with_placements.items():
            alias_hash[file_alias]['place'] = placement_entries_list

        for k, v in alias_hash.items():
            if 'place' not in v:
                alias_hash[k]['place'] = []
        self.write_jplace(jplace_json,
                          alias_hash)

        self.hk.delete(files_to_delete)# Remove combined split, not really useful

        return trusted_placements class Create:
    _PROTEIN_PACKAGE_TYPE = 'protein_package_type'
    _NUCLEOTIDE_PACKAGE_TYPE = 'nucleotide_package_type'

    def __init__(self, commands):
        '''
        Class for managing the creation of GraftM packages.

        Parameters
        ----------
        commands: ExternalProgramSuite object
            Contains each of the command line commands as attributes. e.g.
            fasttree:
                commands.fasttree = "FastTreeMP" or "fasttree" depending on
                what is installed.

        '''
        self.the_trash=[]
        self.fasttree = commands.fasttree

    def _parse_contents(self, contents_file_path):
        '''
        Parse the contents .json file and return the dictionary

        Parameters
        ----------
        contents_file_path: str
            Path to file containing .json file to parse.

        Returns
        -------
        contents_dict: dict
            parsed .json file
        '''
        logging.debug("Parsing %s" % (contents_file_path))
        contents_dict = json.load(open(contents_file_path))
        return contents_dict

    def _check_reads_hit(self, alignment_io, min_aligned_fraction):
        '''Given an alignment return a list of sequence names that are less
        than the min_aligned_fraction'''
        to_return = []
        alignment_length = None
        for s in SeqIO.parse(alignment_io, "fasta"):
            if not alignment_length:
                alignment_length = len(s.seq)
                min_length = int(min_aligned_fraction * alignment_length)
                logging.debug("Determined min number of aligned bases to be %s" % min_length)
            elif len(s.seq) != alignment_length:
                raise Exception("Alignment file appears to not be of uniform length")

            num_unaligned = s.seq.count('-')
            num_aligned = alignment_length-num_unaligned
            logging.debug("Sequence %s has %d aligned positions" % (s.name, alignment_length-num_unaligned))
            if num_aligned <= min_length:
                to_return.append(s.name)
        return to_return

    def _align_sequences(self, input_sequences_path, output_alignment_path,
                         threads):
        '''Align sequences into alignment_file

        Parameters
        ----------
        input_sequences_path: str
            path to input sequences in fasta format
        output_alignment_path: str
            path to output alignment path
        threads: str
            number of threads to use
        Returns
        -------
        Nothing
        '''
        logging.debug("Aligning sequences using mafft")
        cmd = "mafft --anysymbol --thread %s --auto /dev/stdin > %s" % (
            threads,
            output_alignment_path)
        inputs = []
        with open(input_sequences_path) as f:
            for name,seq,_ in SequenceIO().each(f):
                inputs.append('>%s' % name)
                # Do not include * characters in the HMM, as this means tree
                # insertion fails.
                inputs.append(seq.replace('*',''))
        extern.run(cmd, stdin="\n".join(inputs))

    def _get_hmm_from_alignment(self, alignment, hmm_filename, output_alignment_filename):
        '''Return a HMM file and alignment of sequences to that HMM

        Parameters
        ----------
        alignment: str
            path to aligned proteins
        hmm_filename: str
            write the hmm to this file path
        output_alignment_filename: str
            write the output alignment to this file path

        Returns
        -------
        Return the pipeline type of the HMM.
        '''
        logging.info("Building HMM from alignment")

        with tempfile.NamedTemporaryFile(
                suffix='.fasta',prefix='graftm',mode='w') as tempaln:

            cmd = "hmmbuild -O /dev/stdout -o /dev/stderr '%s' '%s'" % (
                hmm_filename, alignment)
            output = extern.run(cmd)

            SeqIO.write(SeqIO.parse(StringIO(output), 'stockholm'), tempaln, 'fasta')
            tempaln.flush()

            ptype, _ = self._pipe_type(hmm_filename)
            SequenceSearcher(hmm_filename).alignment_correcter([tempaln.name],
                                                    output_alignment_filename)
            return ptype

    def _align_sequences_to_hmm(self, hmm_file, sequences_file, output_alignment_file):
        '''Align sequences to an HMM, and write an alignment of
        these proteins after cleanup so that they can be used for tree-making

        Parameters
        ----------
        sequences_file: str
            path to file of unaligned protein sequences
        hmm_file: str
            path to hmm file
        output_alignment_file: str
            write alignment to this file

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

        ss = SequenceSearcher(hmm_file)
        with tempfile.NamedTemporaryFile(prefix='graftm', suffix='.aln.fasta') as tempalign:
            ss.hmmalign_sequences(hmm_file, sequences_file, tempalign.name)
            ss.alignment_correcter([tempalign.name], output_alignment_file)

    def _pipe_type(self, hmm):
        logging.debug("Setting pipeline type.")
        with open(hmm) as f:
            hmm_type=[x.split() for x in f.readlines() if x.startswith('ALPH') or x.startswith('LENG')]
        for item in hmm_type:
            if item[0]=='ALPH':
                if item[1]=='amino':
                    ptype=Create._PROTEIN_PACKAGE_TYPE
                elif item[1]=='DNA' or 'RNA':
                    ptype=Create._NUCLEOTIDE_PACKAGE_TYPE
                else:
                    raise Exception("Unfamiliar HMM type: %s" % (item[1]))
            elif item[0]=='LENG':
                leng=item[1]
            else:
                raise Exception("Programming Error: Misread HMM file")
        logging.debug("Set pipeline type as: %s " % ptype)
        logging.debug("Found alignment length as: %s" % leng)
        return ptype, leng

    def _check_aln_length(self, alignment):
        return len(list(SeqIO.parse(open(alignment, 'r'), 'fasta'))[0].seq)

    def _build_tree(self, alignment, base, ptype, fasttree):
        log_file = base + ".tre.log"
        tre_file = base + ".tre"
        if ptype == Create._NUCLEOTIDE_PACKAGE_TYPE: # If it's a nucleotide sequence
            cmd = "%s -quiet -gtr -nt -log %s -out %s %s" % (fasttree,
                                                             log_file,
                                                             tre_file,
                                                             alignment)
            extern.run(cmd)
        else: # Or if its an amino acid sequence
            cmd = "%s -quiet -log %s -out %s %s" % (fasttree,
                                                    log_file,
                                                    tre_file,
                                                    alignment)
            extern.run(cmd)

        self.the_trash += [log_file, tre_file]
        return log_file, tre_file

    def _taxit_create(self, base, aln_file, tre, tre_log, tax, seq,
                      refpkg, no_reroot):
        cmd = "taxit create -f %s -P %s -t %s -s %s -c -l  %s -T %s -i %s"\
            % (aln_file, refpkg, tre, tre_log, base, tax, seq)

        if no_reroot:
            cmd += ' --no-reroot'
            logging.debug("Calling command assuming pre-rerooting: %s" % cmd)

            extern.run(cmd)

        else:
            logging.debug("Calling command: %s" % cmd)
            logging.info("Attempting to run taxit create with rerooting capabilities")
            def exit_gracefully(base):
                local_base = os.path.basename(base)
                nontemp_tax_file = 'graftm_create_taxonomy.%s.csv' % local_base
                shutil.copy(tax, nontemp_tax_file)
                nontemp_seqinfo_file = 'graftm_create_seqinfo.%s.csv' % local_base
                shutil.copy(seq, nontemp_seqinfo_file)
                nontemp_aln_file = 'graftm_create_alignment.%s.faa' % local_base
                shutil.copy(aln_file, nontemp_aln_file)
                nontemp_tre_file = 'graftm_create_tree.%s.tree' % local_base
                shutil.copy(tre, nontemp_tre_file)
                logging.error('''
taxit create failed to run in a small amount of time suggesting that
rerooting was unsuccessful. Unfortunately this tree will need to be rerooted
manually yourself using a tree editor such as ARB or FigTree.
Once you have a rerooted newick format tree, rerun graftm create
specifying the new tree with --rerooted_tree. The tree file to be rerooted is \'%s\'''' % nontemp_tre_file)

                logging.error('''
When rerunning, please use the following flags for the command line to account
for the fact that some sequences may have been removed during the deduplication
process.

graftM create --taxtastic_taxonomy %s --taxtastic_seqinfo %s --alignment %s  --rerooted_tree <REROOTED_TREE>

(plus other relevant arguments).
''' % (
       #nontemp_align_hmm_file,
       #nontemp_search_hmm_file,
       nontemp_tax_file,
       nontemp_seqinfo_file,
       nontemp_aln_file,
       #nontemp_seq_file
       )
                              )

                exit(2)


            try:
                process = subprocess.Popen(['bash','-c',cmd],
                                           stdout=subprocess.PIPE,
                                           stderr=subprocess.PIPE,
                                           preexec_fn=os.setsid)
                stdout, stderr = process.communicate(timeout=60)
                logging.debug("taxit create STDOUT was: %s" % stdout)
                logging.debug("taxit create STDERR was: %s" % stderr)

                if process.returncode != 0 or re.match(r'^Error: ',stdout.decode()):
                    raise extern.ExternCalledProcessError(
                        process, cmd)
            except subprocess.TimeoutExpired:
                os.killpg(process.pid, signal.SIGTERM)
                exit_gracefully(base)
            except extern.ExternCalledProcessError:
                exit_gracefully(base)

        return refpkg
    
    def _taxit_create_no_tree(self, base, aln_file, tax, seq, refpkg):
        cmd = "taxit create -f %s -P %s -c -l  %s -T %s -i %s"\
            % (aln_file, refpkg, base, tax, seq)

        logging.debug("Calling command without tree: %s" % cmd)
        extern.run(cmd)

        return refpkg

    def _cleanup(self, the_trashcan):
        for every_piece_of_junk in the_trashcan:
            if os.path.isdir(every_piece_of_junk):
                shutil.rmtree(every_piece_of_junk)
            else:
                os.remove(every_piece_of_junk)

    def _define_range(self, sequences):
        '''
        define_range - define the maximum range within which two hits in a db
        search can be linked. This is defined as 1.5X the average length of all
        reads in the database.

        Parameters
        ----------
        sequences : str
            A path to the sequences in FASTA format. This fasta file is assumed
            to be in the correct format. i.e. headers start with '>'

        Returns
        -------
        max_range : int
            As described above, 1.5X the size of the average length of genes
            within the database
        '''
        sequence_count = 0
        total_sequence = 0

        with open(sequences) as f:
            for record in SeqIO.parse(f, 'fasta'): # For every sequence in the output
                total_sequence += 1 # increment the total sequence count
                sequence_count += len(record.seq) # add the length of that sequence to the total string count

        # Get the average, multiply by 1.5, and return
        max_range = (sequence_count/total_sequence)*1.5
        return max_range

    def _generate_tree_log_file(self, tree, alignment, output_tree_file_path,
                               output_log_file_path, residue_type, fasttree):
        '''Generate the FastTree log file given a tree and the alignment that
        made that tree

        Returns
        -------
        Nothing. The log file as parameter is written as the log file.
        '''
        if residue_type==Create._NUCLEOTIDE_PACKAGE_TYPE:
            cmd = "%s -quiet -gtr -nt -nome -mllen -intree '%s' -log %s -out %s %s" %\
                                       (fasttree, tree, output_log_file_path,
                                        output_tree_file_path, alignment)
        elif residue_type==Create._PROTEIN_PACKAGE_TYPE:
            cmd = "%s -quiet -nome -mllen -intree '%s' -log %s -out %s %s" %\
                                       (fasttree, tree, output_log_file_path,
                                        output_tree_file_path, alignment)
        extern.run(cmd)

    def _remove_sequences_from_alignment(self, sequence_names, input_alignment_file, output_alignment_file):
        '''Remove sequences from the alignment file that have names in
        sequence_names

        Parameters
        ----------
        sequence_names: list of str
            names of sequences to remove
        input_alignment_file: str
            path to alignment file to remove from
        output_alignment_file: str
            path to alignment file to write to

        Returns
        -------
        int: number of sequences written to file'''
        nameset = set(sequence_names)
        num_written = 0
        with open(output_alignment_file, 'w') as output:
            with open(input_alignment_file) as input_f:
                for s in SeqIO.parse(input_f, "fasta"):
                    if s.name not in nameset:
                        SeqIO.write(s, output, "fasta")
                        num_written += 1
        logging.debug("After removing sequences from alignment, %i remain" % num_written)
        return num_written

    def _create_search_hmm(self, sequences, taxonomy_definition,
                           search_hmm, dereplication_level, threads):
        logging.debug("Constructing search HMM")
        with tempfile.NamedTemporaryFile(prefix='graftm', suffix='.aln.faa') as temporary_alignment_fh:
            temporary_alignment = temporary_alignment_fh.name
            with tempfile.NamedTemporaryFile(prefix='graftm', suffix='.faa') as temporary_dereplicated_fh:
                temporary_dereplicated = temporary_dereplicated_fh.name
                if dereplication_level > 0:
                    dereplication_index = dereplication_level - 1

                    #####################
                    ### Dereplication ###
                    logging.debug("Dereplicating sequences at rank %i in the taxonomy file provided" % (dereplication_level))
                    dereplicated_sequence_ids = []
                    seen_taxons = set()
                    prefix_re = re.compile(r'^[a-zA-Z]__$')
                    for sequence_id, taxonomy in taxonomy_definition.items():
                        try:
                            taxon = taxonomy[dereplication_index]
                        except IndexError:
                            taxon = None
                        if taxon is None or \
                           taxon == "" or \
                           prefix_re.match(taxon):
                            pass

                        elif taxon in seen_taxons:
                            logging.debug("Sequence %s redundant at %i rank in the taxonomy file level: %s" % (sequence_id, dereplication_level, taxon) )
                            dereplicated_sequence_ids.append(sequence_id)
                        else:
                            seen_taxons.add(taxon)
                    if len(taxonomy_definition) - len(dereplicated_sequence_ids) < 2:
                        raise Exception("Insufficient sequences available after dereplication to create a HMM. One solution might be to modify the dereplication level to do less dereplication.")
                    logging.info("Removing %i sequences from the search HMM that are redundant at the %i rank in the taxonomy file" \
                                            % (len(dereplicated_sequence_ids),
                                               dereplication_level))
                    self._remove_sequences_from_alignment(dereplicated_sequence_ids, sequences, temporary_dereplicated)
                    self._align_sequences(temporary_dereplicated, temporary_alignment, threads)
                else:
                    logging.debug("Skipping dereplication step and using all sequences to build search HMM")
                    self._align_sequences(sequences, temporary_alignment, threads)

                self._get_hmm_from_alignment(temporary_alignment, search_hmm, temporary_alignment)


    def _create_dmnd_database(self, unaligned_sequences_path, daa_output):
        '''
        Build a diamond database using diamond makedb

        Parameters
        ----------
        unaligned_sequences_path: str
            path to a FASTA file containing unaligned sequences
        daa_output: str
            Name of output database.
        '''
        logging.debug("Building diamond database")

        cmd = "diamond makedb --in '%s' -d '%s'" % (unaligned_sequences_path, daa_output)
        extern.run(cmd)

    def _align_and_create_hmm(self, sequences, alignment, user_hmm,
                              output_align_hmm, output_alignment, threads):

        # align sequences to HMM (and potentially build hmm from alignment)
        if user_hmm:
            output_hmm = user_hmm
            ptype, _ = self._pipe_type(output_hmm)
            if alignment:
                logging.info("Using pre-computed alignment, and assuming it \
                    was generated by aligning it to the user-specified HMM ..")
                output_alignment = alignment
            else:
                self._align_sequences_to_hmm(output_hmm, sequences, output_alignment)
        else:
            if not alignment:
                self._align_sequences(sequences, output_alignment, threads)
            else:
                output_alignment = alignment
            ptype = self._get_hmm_from_alignment(output_alignment,
                                                 output_align_hmm,
                                                 output_alignment)
        return ptype, output_alignment

    def _mask_strange_sequence_letters(self, sequences, package_type):
        '''Replace strange characters like selenocysteine (U) with X or N for
        protein and nucleotide sequences, respectively. Sequences are
        modified in place.

        Parameters
        ----------
        aligned_sequences: array of Sequence objects
            sequences to mask
        package_type: _PROTEIN_PACKAGE_TYPE or _NUCLEOTIDE_PACKAGE_TYPE
            type of sequences these are

        Returns
        -------
        None
        '''
        if package_type == Create._PROTEIN_PACKAGE_TYPE:
            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
        elif taxtastic_seqinfo and taxtastic_taxonomy:
            logging.info("Reading taxonomy from taxtastic taxonomy and seqinfo files")
            with open(taxtastic_taxonomy) as taxf:
                with open(taxtastic_seqinfo) as seqf:
                    taxonomy_definition = \
                        gtns.read_taxtastic_taxonomy_and_seqinfo(taxf, seqf)
        else:
            raise Exception("Taxonomy is required somehow e.g. by --taxonomy or --rerooted_annotated_tree")

        # Check for duplicates
        logging.info("Checking for duplicate sequence names")
        dup = self._check_for_duplicate_sequence_names(sequences)
        if dup:
            raise Exception("Found duplicate sequence name '%s' in sequences input file" % dup)
        output_alignment_fh = tempfile.NamedTemporaryFile(prefix='graftm', suffix='.aln.faa')
        tempfiles_to_close.append(output_alignment_fh)
        output_alignment = output_alignment_fh.name
        if user_hmm:
            align_hmm = user_hmm
        else:
            align_hmm_fh = tempfile.NamedTemporaryFile(prefix='graftm', suffix='_align.hmm')
            tempfiles_to_close.append(align_hmm_fh)
            align_hmm = align_hmm_fh.name

        if alignment:
            dup = self._check_for_duplicate_sequence_names(alignment)
            
            output_alignment_tmpfile = tempfile.NamedTemporaryFile(prefix='graftm-create-align',suffix='.fa')
            tempfiles_to_close.append(output_alignment_tmpfile)
            output_alignment = output_alignment_tmpfile.name

            if dup:
                raise Exception("Found duplicate sequence name '%s' in alignment input file" % dup)
            if user_hmm:
                ptype, _ = self._pipe_type(align_hmm)
            else:
                ptype = self._get_hmm_from_alignment(alignment,
                                                    align_hmm,
                                                    output_alignment)
        else:
            logging.info("Aligning sequences to create aligned FASTA file")
            ptype, output_alignment = self._align_and_create_hmm(sequences, alignment, user_hmm,
                                               align_hmm, output_alignment, threads)

        logging.info("Checking for incorrect or fragmented reads")
        with open(output_alignment) as f:
            insufficiently_aligned_sequences = self._check_reads_hit(
                f, min_aligned_percent)
        while len(insufficiently_aligned_sequences) > 0:
            logging.warning("One or more alignments do not span > %.2f %% of HMM" % (min_aligned_percent*100))
            for s in insufficiently_aligned_sequences:
                logging.warning("Insufficient alignment of %s, not including this sequence" % s)

            sequences2_fh = tempfile.NamedTemporaryFile(prefix='graftm', suffix='.faa')
            tempfiles_to_close.append(sequences2_fh)
            sequences2 = sequences2_fh.name
            num_sequences = self._remove_sequences_from_alignment(insufficiently_aligned_sequences,
                                                                  sequences,
                                                                  sequences2)
            sequences = sequences2

            if alignment:
                alignment2_fh = tempfile.NamedTemporaryFile(prefix='graftm', suffix='.aln.faa')
                tempfiles_to_close.append(alignment2_fh)
                alignment2 = alignment2_fh.name
                num_sequences = self._remove_sequences_from_alignment(insufficiently_aligned_sequences,
                                                                      alignment,
                                                                      alignment2)
                alignment = alignment2
                for name in insufficiently_aligned_sequences:
                    if rerooted_tree or rerooted_annotated_tree:
                        logging.warning('''Sequence %s in provided alignment does not meet the --min_aligned_percent cutoff. This sequence will be removed from the tree
in the final GraftM package. If you are sure these sequences are correct, turn off the --min_aligned_percent cutoff, provide it with a 0 (e.g. --min_aligned_percent 0) ''' % name)
                    removed_sequence_names.append(name)


            logging.info("After removing %i insufficiently aligned sequences, left with %i sequences" % (len(insufficiently_aligned_sequences), num_sequences))
            if num_sequences < 4:
                raise Exception("Too few sequences remaining in alignment after removing insufficiently aligned sequences: %i" % num_sequences)
            else:
                logging.info("Reconstructing the alignment and HMM from remaining sequences")
                output_alignment_fh = tempfile.NamedTemporaryFile(prefix='graftm', suffix='.aln.faa')
                tempfiles_to_close.append(output_alignment_fh)
                output_alignment = output_alignment_fh.name
                if not user_hmm:
                    align_hmm_fh = tempfile.NamedTemporaryFile(prefix='graftm', suffix='.hmm')
                    tempfiles_to_close.append(align_hmm_fh)
                    align_hmm = align_hmm_fh.name
                ptype, output_alignment= self._align_and_create_hmm(sequences, alignment, user_hmm,
                                                   align_hmm, output_alignment, threads)
                logging.info("Checking for incorrect or fragmented reads")
                insufficiently_aligned_sequences = self._check_reads_hit(open(output_alignment),
                                                                         min_aligned_percent)
        if not search_hmm_files:
            search_hmm_fh = tempfile.NamedTemporaryFile(prefix='graftm', suffix='_search.hmm')
            tempfiles_to_close.append(search_hmm_fh)
            search_hmm = search_hmm_fh.name
            self._create_search_hmm(sequences, taxonomy_definition, search_hmm, dereplication_level, threads)
            search_hmm_files = [search_hmm]

        # Make sure each sequence has been assigned a taxonomy:
        aligned_sequence_objects = seqio.read_fasta_file(output_alignment)
        unannotated = []
        for s in aligned_sequence_objects:
            if s.name not in taxonomy_definition:
                unannotated.append(s.name)
        if len(unannotated) > 0:
            for s in unannotated:
                logging.error("Unable to find sequence '%s' in the taxonomy definition" % s)
            raise Exception("All sequences must be assigned a taxonomy, cannot continue")


        logging.debug("Looking for non-standard characters in aligned sequences")
        self._mask_strange_sequence_letters(aligned_sequence_objects, ptype)

        # Deduplicate sequences - pplacer cannot handle these
        logging.info("Deduplicating sequences")
        dedup = Deduplicator()
        deduplicated_arrays = dedup.deduplicate(aligned_sequence_objects)
        deduplicated_taxonomy = dedup.lca_taxonomy(deduplicated_arrays, taxonomy_definition)
        deduplicated_taxonomy_hash = {}
        for i, tax in enumerate(deduplicated_taxonomy):
            deduplicated_taxonomy_hash[deduplicated_arrays[i][0].name] = tax
        deduplicated_alignment_file = base+"_deduplicated_aligned.fasta"
        seqio.write_fasta_file([seqs[0] for seqs in deduplicated_arrays],
                               deduplicated_alignment_file)

        logging.info("Removed %i sequences as duplicates, leaving %i non-identical sequences"\
                     % ((len(aligned_sequence_objects)-len(deduplicated_arrays)),
                        len(deduplicated_arrays)))

        # Get corresponding unaligned sequences
        filtered_names=[]
        for list in [x for x in [x[1:] for x in deduplicated_arrays] if x]:
            for seq in list:
                filtered_names.append(seq.name)
        sequences2_fh = tempfile.NamedTemporaryFile(prefix='graftm', suffix='.faa')
        tempfiles_to_close.append(sequences2_fh)
        sequences2 = sequences2_fh.name


        # Create tree unless one was provided
        if not rerooted_tree and not rerooted_annotated_tree and not unrooted_tree and not no_tree:
            logging.debug("No tree provided")
            logging.info("Building tree")
            log_file, tre_file = self._build_tree(deduplicated_alignment_file,
                                                  base, ptype,
                                                  self.fasttree)
            no_reroot = False
        elif no_tree:
            logging.info("Tree-less package requested")
        else:
            if rerooted_tree:
                logging.debug("Found unannotated pre-rerooted tree file %s" % rerooted_tree)
                tre_file=rerooted_tree
                no_reroot = True
            elif rerooted_annotated_tree:
                logging.debug("Found annotated pre-rerooted tree file %s" % rerooted_tree)
                tre_file=rerooted_annotated_tree
                no_reroot = True
            elif unrooted_tree:
                logging.info("Using input unrooted tree")
                tre_file = unrooted_tree
                no_reroot = False
            else:
                raise

            # Remove any sequences from the tree that are duplicates
            logging.info("Removing duplicates from tree")
            cleaner = DendropyTreeCleaner()
            tree = Tree.get(path=tre_file, schema='newick')
            for group in deduplicated_arrays:
                [removed_sequence_names.append(s.name) for s in group[1:]]
            cleaner.remove_sequences(tree, removed_sequence_names)

            # Ensure there is nothing amiss now as a user-interface thing
            cleaner.match_alignment_and_tree_sequence_ids(\
                [g[0].name for g in deduplicated_arrays], tree)

            if tree_log:
                # User specified a log file, go with that
                logging.debug("Using user-specified log file %s" % tree_log)
                log_file = tree_log
            else:
                logging.info("Generating log file")
                log_file_tempfile = tempfile.NamedTemporaryFile(suffix='.tree_log', prefix='graftm')
                tempfiles_to_close.append(log_file_tempfile)
                log_file = log_file_tempfile.name
                tre_file_tempfile = tempfile.NamedTemporaryFile(suffix='.tree', prefix='graftm')
                tempfiles_to_close.append(tre_file_tempfile)
                tre_file = tre_file_tempfile.name
                with tempfile.NamedTemporaryFile(suffix='.tree', prefix='graftm') as f:
                    # Make the newick file simple (ie. un-arb it) for fasttree.
                    cleaner.write_fasttree_newick(tree, f)
                    f.flush()
                    self._generate_tree_log_file(f.name, deduplicated_alignment_file,
                                                 tre_file, log_file, ptype, self.fasttree)

        # Create tax and seqinfo .csv files
        taxonomy_to_keep=[
                          seq.name for seq in
                                [x for x in [x[0] for x in deduplicated_arrays]
                                 if x]
                          ]
        refpkg = "%s.refpkg" % output_gpkg_path
        self.the_trash.append(refpkg)
        if taxtastic_taxonomy and taxtastic_seqinfo:
            logging.info("Creating reference package")
            if no_tree:
                refpkg = self._taxit_create_no_tree(base, deduplicated_alignment_file,
                                            taxtastic_taxonomy, taxtastic_seqinfo,
                                            refpkg)
            else:
                refpkg = self._taxit_create(base, deduplicated_alignment_file,
                                            tre_file, log_file, taxtastic_taxonomy,
                                            taxtastic_seqinfo, refpkg, no_reroot)
        else:
            gtns = Getaxnseq()
            seq = base+"_seqinfo.csv"
            tax = base+"_taxonomy.csv"
            self.the_trash += [seq, tax]
            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
            else:
                raise Exception("Programming error: Taxonomy is required somehow e.g. by --taxonomy or --rerooted_annotated_tree")

            taxonomy_definition = {x:taxonomy_definition[x]
                                   for x in taxonomy_definition
                                   if x in taxonomy_to_keep}

            gtns.write_taxonomy_and_seqinfo_files(taxonomy_definition,
                                                  tax,
                                                  seq)

            # Create the reference package
            logging.info("Creating reference package")
            if no_tree:
                refpkg = self._taxit_create_no_tree(base, deduplicated_alignment_file,
                                                    tax, seq, refpkg)
            else:
                refpkg = self._taxit_create(base, deduplicated_alignment_file,
                                            tre_file, log_file, tax, seq, refpkg,
                                            no_reroot)
        if sequences:
            # Run diamond makedb
            logging.info("Creating diamond database")
            if ptype == Create._PROTEIN_PACKAGE_TYPE:
                cmd = "diamond makedb --in '%s' -d '%s'" % (sequences, base)
                extern.run(cmd)
                diamondb = '%s.dmnd' % base
            elif ptype == Create._NUCLEOTIDE_PACKAGE_TYPE:
                diamondb = None
            else: raise Exception("Programming error")
        else:
            diamondb = None

        if sequences:
            # Get range
            max_range = self._define_range(sequences)
        else:
            max_range = self._define_range(alignment)

        # Compile the gpkg
        logging.info("Compiling gpkg")

        GraftMPackageVersion3.compile(output_gpkg_path, refpkg, align_hmm, diamondb,
                                      max_range, sequences, search_hmm_files=search_hmm_files)

        logging.info("Cleaning up")
        self._cleanup(self.the_trash)
        for tf in tempfiles_to_close:
            tf.close()

        # Test out the gpkg just to be sure.
        #
        # TODO: Use graftM through internal means rather than via extern. This
        # requires some refactoring so that graft() can be called easily with
        # sane defaults.
        if not no_tree:
            logging.info("Testing gpkg package works")
            self._test_package(output_gpkg_path)

        logging.info("Finished\n") class Update(Create):
    def _concatenate_file(self, file_list, output):
        '''
        Call unix "cat" to concatenate a list of files

        Parameters
        ----------
        file_list: list
            List of strings, each leading to a file. These files are the ones to
            be concatenate together. E.g.:
                ["/path/to/file1", "/path/to/file2"]
        output: str
            Path to file to which to the files in file_list will be concatenated
            into.

        '''
        to_cat = ' '.join(file_list)
        logging.debug("Concatenating files: %s" % (to_cat))
        cmd = "cat %s > %s" % (to_cat, output)
        extern.run(cmd)

    def update(self, **kwargs):
        '''
        Update an existing GraftM package with new sequences and taxonomy. If no
        taxonomy is provided, attempt to decorate the new sequences with
        pre-existing taxonomy.

        Parameters
        ----------
        input_sequence_path: str
            Path to FASTA file containing sequences to add to the update GraftM
            package
        input_taxonomy_path: str
            Taxonomy corresponding to the sequences in input_sequence_path. If None,
            then attempt to assign taxonomy by decorating the tree made out of all
            sequences.
        input_graftm_package_path: str
            Path to the directory of the GraftM package that is to be updated
        output_graftm_package_path: str
            Path to the directory to which the new GraftM package will be
            written to
        '''
        input_sequence_path = kwargs.pop('input_sequence_path')
        input_taxonomy_path = kwargs.pop('input_taxonomy_path', None)
        input_graftm_package_path = kwargs.pop('input_graftm_package_path')
        output_graftm_package_path = kwargs.pop('output_graftm_package_path')
        threads = kwargs.pop('threads', UpdateDefaultOptions.threads) #TODO: add to user options
        if len(kwargs) > 0:
            raise Exception("Unexpected arguments detected: %s" % kwargs)

        logging.info("Reading previous GraftM package")
        old_gpkg = GraftMPackage.acquire(input_graftm_package_path)
        min_input_version = 3
        if old_gpkg.version < min_input_version:
            raise InsufficientGraftMPackageVersion(
                "GraftM below version %s cannot be updated using the update function." % min_input_version +
                " Unaligned sequences are not included in these packages, therefore no new"
                " alignment/HMM/Tree can be created")

        new_gpkg = UpdatedGraftMPackage()
        new_gpkg.output = output_graftm_package_path
        new_gpkg.name = output_graftm_package_path.replace(".gpkg", "")

        #######################################
        ### Collect all unaligned sequences ###
        logging.info("Concatenating unaligned sequence files")
        new_gpkg.unaligned_sequences = "%s_sequences.fa" % (new_gpkg.name) #TODO: replace hard-coded paths like this with tempfiles
        self._concatenate_file([old_gpkg.unaligned_sequence_database_path(),
                                input_sequence_path],
                               new_gpkg.unaligned_sequences)

        #########################################################
        ### Parse taxonomy info up front so errors come early ###
        if input_taxonomy_path:
            logging.info("Reading new taxonomy information")
            input_taxonomy = GreenGenesTaxonomy.read_file(input_taxonomy_path)
            original_taxonomy_hash = old_gpkg.taxonomy_hash()
            total_taxonomy_hash = original_taxonomy_hash.copy()
            total_taxonomy_hash.update(input_taxonomy.taxonomy)
            num_duplicate_taxonomies = len(total_taxonomy_hash) - \
                                       len(input_taxonomy.taxonomy) - \
                                       len(original_taxonomy_hash)
            logging.debug("Found %i taxonomic definitions in common between the previous and updated taxonomies" % num_duplicate_taxonomies)
            if num_duplicate_taxonomies > 0:
                logging.warn("Found %i taxonomic definitions in common between the previous and updated taxonomies. Using the updated taxonomy in each case." % num_duplicate_taxonomies)

        ###############################
        ### Re-construct alignments ###
        logging.info("Multiple sequence aligning all sequences")
        new_gpkg.aligned_sequences = "%s_mafft_alignment.fa" % (new_gpkg.name)
        self._align_sequences(new_gpkg.unaligned_sequences, new_gpkg.aligned_sequences, threads)

        ########################
        ### Re-construct HMM ###
        logging.info("Creating HMM from alignment")
        new_gpkg.hmm = "%s.hmm" % (new_gpkg.name)
        new_gpkg.hmm_alignment = "%s_hmm_alignment.fa" % (new_gpkg.name)
        self._get_hmm_from_alignment(new_gpkg.aligned_sequences, new_gpkg.hmm, new_gpkg.hmm_alignment)

        #########################
        ### Re-construct tree ###
        logging.info("Generating phylogenetic tree")
        new_gpkg.unrooted_tree = "%s.tre" % (new_gpkg.name)
        new_gpkg.unrooted_tree_log = "%s.tre.log" % (new_gpkg.name)
        new_gpkg.package_type, new_gpkg.hmm_length = self._pipe_type(old_gpkg.alignment_hmm_path())
        new_gpkg.unrooted_gpkg_tree_log, new_gpkg.unrooted_gpkg_tree = \
            self._build_tree(new_gpkg.hmm_alignment, new_gpkg.name,
                             new_gpkg.package_type, self.fasttree)

        ##############################################
        ### Re-root and decorate tree if necessary ###
        if input_taxonomy_path:
            new_gpkg.gpkg_tree_log = new_gpkg.unrooted_tree_log
            new_gpkg.gpkg_tree = new_gpkg.unrooted_gpkg_tree
        else:
            logging.info("Finding taxonomy for new sequences")
            rerooter = Rerooter()
            
            old_tree = Tree.get(path=old_gpkg.reference_package_tree_path(),
                                schema='newick')
            new_tree = Tree.get(path=new_gpkg.unrooted_gpkg_tree,
                                schema='newick')
            old_tree = rerooter.reroot(old_tree)
            new_tree = rerooter.reroot(new_tree)
            # TODO: Shouldn't call an underscore method, eventually use
            # Rerooter instead.
            rerooted_tree = rerooter.reroot_by_tree(old_tree, new_tree)
            new_gpkg.gpkg_tree = "%s_gpkg.tree" % new_gpkg.name
            td = TreeDecorator(
                rerooted_tree,
                old_gpkg.taxtastic_taxonomy_path(),
                old_gpkg.taxtastic_seqinfo_path())
            
            with tempfile.NamedTemporaryFile(suffix='tsv') as taxonomy:
                td.decorate(new_gpkg.gpkg_tree, taxonomy.name, True) 
                total_taxonomy_hash = GreenGenesTaxonomy.read_file(taxonomy.name).taxonomy

            ################################
            ### Generating tree log file ###
            logging.info("Generating phylogenetic tree log file")
            new_gpkg.gpkg_tree = "%s_gpkg.tree" % new_gpkg.name
            new_gpkg.gpkg_tree_log = "%s_gpkg.tree.log" % new_gpkg.name
            self._generate_tree_log_file(new_gpkg.unrooted_tree,
                                         new_gpkg.hmm_alignment,
                                         new_gpkg.gpkg_tree,
                                         new_gpkg.gpkg_tree_log,
                                         new_gpkg.package_type,
                                         self.fasttree)

        ################################
        ### Creating taxtastic files ###
        logging.info("Writing new taxonomy files")
        new_gpkg.tt_seqinfo = "%s_seqinfo.csv" % new_gpkg.name
        new_gpkg.tt_taxonomy = "%s_taxonomy.csv" % new_gpkg.name
        gtns = Getaxnseq()

        gtns.write_taxonomy_and_seqinfo_files(
            total_taxonomy_hash,
            new_gpkg.tt_taxonomy,
            new_gpkg.tt_seqinfo)
        
        ######################
        ### Compile refpkg ###
        logging.info("Compiling pplacer refpkg")
        new_gpkg.refpkg = "%s.refpkg" % (new_gpkg.name)
        refpkg = self._taxit_create(new_gpkg.name,
                                    new_gpkg.hmm_alignment,
                                    new_gpkg.gpkg_tree,
                                    new_gpkg.gpkg_tree_log,
                                    new_gpkg.tt_taxonomy,
                                    new_gpkg.tt_seqinfo,
                                    new_gpkg.refpkg,
                                    True)

        #####################################
        ### Re-construct diamond database ###
        logging.info("Recreating DIAMOND DB")
        new_gpkg.diamond_database = "%s.dmnd" % (new_gpkg.name)
        self._create_dmnd_database(new_gpkg.unaligned_sequences, new_gpkg.name)

        ####################
        ### Compile gpkg ###
        logging.info("Compiling GraftM package")
        new_gpkg.name = "%s.gpkg" % new_gpkg.name
        GraftMPackageVersion3.compile(new_gpkg.name, new_gpkg.refpkg,
                                      new_gpkg.hmm, new_gpkg.diamond_database,
                                      self._define_range(new_gpkg.unaligned_sequences),
                                      new_gpkg.unaligned_sequences,
                                      search_hmm_files=old_gpkg.search_hmm_paths())

        ###################
        ### Test it out ###
        logging.info("Testing newly updated GraftM package works")
        self._test_package(new_gpkg.name)

        logging.info("Finished") class UnpackRawReads:
    class UnexpectedFileFormatException(Exception): pass

    FORMAT_FASTA    = "FORMAT_FASTA"
    FORMAT_FASTQ    = "FORMAT_FASTQ"
    FORMAT_FASTQ_GZ = "FORMAT_FASTQ_GZ"
    FORMAT_FASTA_GZ = "FORMAT_FASTA_GZ"

    PROTEIN_SEQUENCE_TYPE = 'aminoacid'
    NUCLEOTIDE_SEQUENCE_TYPE = 'nucleotide'

    _EXTENSION_TO_FILE_TYPE = {'.fa': FORMAT_FASTA,
                               '.faa': FORMAT_FASTA,
                               '.fna': FORMAT_FASTA,
                               '.fasta': FORMAT_FASTA,

                               '.fq': FORMAT_FASTQ,
                               '.fastq': FORMAT_FASTQ,

                               '.fq.gz': FORMAT_FASTQ_GZ,
                               '.fastq.gz': FORMAT_FASTQ_GZ,

                               '.fa.gz': FORMAT_FASTA_GZ,
                               '.faa.gz': FORMAT_FASTA_GZ,
                               '.fna.gz': FORMAT_FASTA_GZ,
                               '.fasta.gz': FORMAT_FASTA_GZ,
                               }

    def __init__(self, read_file, known_sequence_type=None, interleaved=False):
        '''New object from a read file.

        read_file: str
            path to input reads
        sequence_type: str:
            PROTEIN_SEQUENCE_TYPE, NUCLEOTIDE_SEQUENCE_TYPE or None
            Whether input is nucleotide, amino acid, or should be guessed by
        peeking at the input sequence file.

        '''
        logging.debug("Loading %s, type %s, interleaved %s", read_file,
                      known_sequence_type, interleaved)
        self.read_file = read_file
        self.known_sequence_type = known_sequence_type
        self.interleaved = interleaved

    def _guess_sequence_type_from_string(self, seq):
        '''Return 'protein' if there is >10% amino acid residues in the
        (string) seq parameter, else 'nucleotide'. Raise Exception if a
        non-standard character is encountered'''
        # Define expected residues for each sequence type
        aa_chars = ['P','V','L','I','M','F','Y','W','H','K','R','Q','N','E','D','S','X','*']
        aas = set(itertools.chain(aa_chars, [str.lower(a) for a in aa_chars]))
        na_chars = ['A', 'T', 'G', 'C', 'N', 'U']
        nas = set(itertools.chain(na_chars, [str.lower(a) for a in na_chars]))

        num_nucleotide = 0
        num_protein = 0
        count = 0
        for residue in seq:
            if residue in nas:
                num_nucleotide += 1
            elif residue in aas:
                num_protein += 1
            else:
                raise Exception('Encountered unexpected character when attempting to guess sequence type: %s' % (residue))
            count += 1
            if count >300: break
        if float(num_protein) / (num_protein+num_nucleotide) > 0.1:
            return self.PROTEIN_SEQUENCE_TYPE
        else:
            return self.NUCLEOTIDE_SEQUENCE_TYPE

    def sequence_type(self):
        '''Guess the type of input sequence provided to graftM (i.e. nucleotide
        or amino acid) and return'''
        if self.known_sequence_type is not None:
            return self.known_sequence_type
        else:
            # If its Gzipped and fastq make a small sample of the sequence to be
            # read
            cmd='%s | head -n 2' % (self.command_line())
            # We cannot use extern here because the zcat has non-zero
            # exitstatus because it is head'd and extern uses bash -o pipefail
            first_seq = subprocess.run(['bash','-c',cmd], stdout=subprocess.PIPE, stderr=subprocess.PIPE).stdout.decode('UTF-8')
            seq = first_seq.split('\n')[1]
            self.type = self._guess_sequence_type_from_string(seq)
            logging.debug("Detected sequence type as %s" % self.type)
            return self.type

    def guess_sequence_input_file_format(self, sequence_file_path):
        '''Given a sequence file, guess the format and return. Raise an
        exception if it cannot be guessed'''
        return self._EXTENSION_TO_FILE_TYPE[self._get_extension(sequence_file_path)]

    def format(self):
        return self.guess_sequence_input_file_format(self.read_file)

    def is_zcattable(self):
        return self.guess_sequence_input_file_format(self.read_file) in \
            (self.FORMAT_FASTA_GZ, self.FORMAT_FASTQ_GZ)

    def basename(self):
        '''Return the name of the file with the '.fasta' or 'fq.gz' etc
        removed'''
        return os.path.basename(self.read_file)[:-len(self._get_extension(self.read_file))]

    def _get_extension(self, sequence_file_path):
        for ext in list(self._EXTENSION_TO_FILE_TYPE.keys()):
            if sequence_file_path.endswith(ext): return ext
        raise self.UnexpectedFileFormatException("Unable to guess file format of sequence file: %s" % sequence_file_path)

    def get_interleaved_cmd(self):
        """ return cmd snippet to make sure interleaved reads have unique
        names """
        return r""" | perl -pe 'if (m/^>/) {$i++; if ($i % 2 == 1) { if (m/^(\S+)(?<!\/1)(\s+\S.*)?(\s*)$/) { $_ = "$1/1$2$3" }} elsif ($i % 2 == 0) { if (m/^(\S+)(?<!\/2)(\s+\S.*)?(\s*)$/) { $_ = "$1/2$2$3" }}}'"""

    def command_line(self):
        '''Return a string to open read files with'''
        file_format=self.guess_sequence_input_file_format(self.read_file)
        logging.debug("Detected file format %s" % file_format)
        if file_format == self.FORMAT_FASTA:
            cmd="""cat '%s'""" % (self.read_file)
        elif file_format == self.FORMAT_FASTQ_GZ:
            cmd="""zcat '%s' | awk '{print ">" substr($0,2);getline;print;getline;getline}' -""" % (self.read_file)
        elif file_format == self.FORMAT_FASTA_GZ:
            cmd="""zcat '%s'""" % (self.read_file)
        elif file_format == self.FORMAT_FASTQ:
            cmd="""awk '{print ">" substr($0,2);getline;print;getline;getline}' '%s'""" % (self.read_file)
        if self.interleaved:
            cmd+=self.get_interleaved_cmd()
        logging.debug("raw read unpacking command chunk: %s" % cmd)
        return cmd

    def get_file_as_process(self):
        return "<(%s)" % (self.command_line()) 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 ExpandSearcher:

    # Define constants.
    DIAMOND_SEARCH_METHOD = "diamond"
    HMM_SEARCH_METHOD = "hmmsearch"

    def __init__(self, **kwargs):
        '''
        Parameters
        ----------
        search_hmm_files: list of str
            list of HMM files to search with
        graftm_package: GraftMPackage
            use the search HMMs from this graftm package in addition to
            those specified in search_hmm_files
        threads, evalue, min_orf_length, restrict_read_length:
            as per hmmer.search_and_extract_orfs_matching_protein_database
        '''
        self.search_hmm_files = kwargs.pop('search_hmm_files',[])
        self.maximum_range = kwargs.pop('maximum_range',None)
        self.threads = kwargs.pop('threads',None)
        self.evalue = kwargs.pop('evalue',None)
        self.min_orf_length = kwargs.pop('min_orf_length',None)
        graftm_package = kwargs.pop('graftm_package',None)
        if len(kwargs) > 0:
            raise Exception("Unexpected arguments detected: %s" % kwargs)

        if graftm_package:
            self.diamond_database = graftm_package.diamond_database_path()
            self.unaligned_sequence_database = graftm_package.unaligned_sequence_database_path()
            if self.search_hmm_files is None:
                self.search_hmm_files = []
            for h in graftm_package.search_hmm_paths():
                if h not in self.search_hmm_files:
                    self.search_hmm_files.append(h)
            if self.maximum_range is None:
                self.maximum_range = graftm_package.maximum_range()
        else:
            self.diamond_database = None
            self.unaligned_sequence_database = None

        if len(kwargs) > 0:
            raise Exception("Unexpected arguments detected: %s" % kwargs)


    def generate_expand_search_database_from_contigs(self, contig_files, output_database_file,
                                                 search_method):
        '''Given a collection of search_hmm_files, search the contigs in
        contig_files, and generate an HMM from the resulting hits, outputting
        it as output_database_file.

        Parameters
        ----------
        contig_files: list of str
            list of files to search
        output_database_file: str
            path to output file
        search_method: str
            "diamond" or "hmmsearch", to specify search method to use and what
            type of database to build.

        Returns
        -------
        True if genes were recovered, else False'''

        ss = SequenceSearcher(self.search_hmm_files)
        seqio = SequenceIO()
        if search_method == self.DIAMOND_SEARCH_METHOD:
            if self.diamond_database == None or self.unaligned_sequence_database == None:
                logging.warning("Cannot expand_search continue with no diamond database or unaligned sequences.")
                return False

        with tempfile.NamedTemporaryFile(prefix='graftm_expand_search_orfs') as orfs:
            logging.info("Finding expand_search hits in provided contigs..")
            for contig_file in contig_files:
                logging.debug("Finding expand_search hits in %s.." % contig_file)
                unpack = UnpackRawReads(contig_file)

                with tempfile.NamedTemporaryFile(prefix='graftm_expand_search') as \
                                                        hit_reads_orfs_fasta:
                    # search and extract matching ORFs
                    with tempfile.NamedTemporaryFile(prefix='graftm_expand_search2') as \
                                                        hmmsearch_output_table:
                        with tempfile.NamedTemporaryFile(prefix='graftm_expand_search3') as \
                                                        hit_reads_fasta:
                            ss.search_and_extract_orfs_matching_protein_database(\
                                    unpack,
                                    search_method,
                                    self.maximum_range,
                                    self.threads,
                                    self.evalue,
                                    self.min_orf_length,
                                    None,
                                    (self.diamond_database if self.diamond_database else None),
                                    hmmsearch_output_table.name,
                                    hit_reads_fasta.name,
                                    hit_reads_orfs_fasta.name,
                                    '')
                    # Append to the file
                    with open(hit_reads_orfs_fasta.name, 'rb') as f:
                        orfs.write(f.read())

            # Now have a fasta file of ORFs.
            # Check to make sure the file is not zero-length
            orfs.flush()



            with tempfile.NamedTemporaryFile(prefix="graftm_expand_search_aln") as aln:

                if search_method == self.HMM_SEARCH_METHOD:

                    # Check that there is more than one sequence to align.
                    if len(seqio.read_fasta_file(orfs.name)) <= 1:# Just to build on this, you need to check if there is > 1 hit
                                                                  # otherwise mafft will fail to align, causing a crash when hmmbuild is
                                                                  # run on an empty file.
                        logging.warn("Failed to find two or more matching ORFs in the expand_search contigs")
                        return False

                    # Run mafft to align them
                    cmd = "mafft --auto %s >%s" % (orfs.name, aln.name)
                    logging.info("Aligning expand_search hits..")
                    extern.run(cmd)

                    # Run hmmbuild to create an HMM
                    cmd = "hmmbuild --amino %s %s >/dev/null" % (output_database_file, aln.name)
                    logging.info("Building HMM from expand_search hits..")

                    extern.run(cmd)

                elif search_method == self.DIAMOND_SEARCH_METHOD:

                    # Concatenate database with existing database
                    with tempfile.NamedTemporaryFile(prefix="concatenated_database") as databasefile:
                        for f in [orfs.name, self.unaligned_sequence_database]:
                            for line in open(f,'rb'):
                                databasefile.write(line)
                        databasefile.flush()

                        # Run diamond make to create a diamond database
                        cmd = "diamond makedb --in '%s' -d '%s'" % (databasefile.name, output_database_file)
                        logging.info("Building a diamond database from expand_search hits..")
                        extern.run(cmd)

                else:
                    raise Exception("Search method not recognised: %s" % search_method)
                    return False

                return True class Diamond:
    def __init__(self, database, threads=None, evalue=None):
        self._database = database
        self._threads = threads
        self._evalue = evalue

    def run(self, input_sequence_file, input_sequence_type, daa_file_basename=None, extra_args=''):
        '''Run input sequences in either blastp or blastx mode against the
        database specified in __init__.

        Parameters
        ----------
        input_sequence_file: str
            path to query sequences
        input_sequence_type: either 'nucleotide' or 'protein'
            the input_sequences are this kind of sequence

        Returns
        -------
        DiamondSearchResult
        '''

        cmd_list = ["diamond"]
        if input_sequence_type == UnpackRawReads.PROTEIN_SEQUENCE_TYPE:
            cmd_list.append('blastp')
        elif input_sequence_type == UnpackRawReads.NUCLEOTIDE_SEQUENCE_TYPE:
            cmd_list.append('blastx')
        else:
            raise Exception("Programming error")

        basename = daa_file_basename
        if basename is None:
            with tempfile.NamedTemporaryFile(prefix='graftm_diamond') as t:
                # we are just stealing the name, don't need the file itself
                basename = t.name

        for c in ['-k 1',
                  "-d",
                    self._database,
                    "-q",
                    "%s" % input_sequence_file,
                    "-a",
                    basename,
                    extra_args]:
            cmd_list.append(c)
        if self._threads:
            cmd_list.append("--threads")
            cmd_list.append(str(self._threads))
        if self._evalue:
            cmd_list.append("--evalue")
            cmd_list.append(str(self._evalue))

        cmd = ' '.join(cmd_list)
        extern.run(cmd)

        daa_name = "%s.daa" % basename
        res = DiamondSearchResult.import_from_daa_file(daa_name)

        if daa_file_basename is None:
            # Diamond makes an extra file, need to remove this
            os.remove(daa_name)

        return res 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 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 Timer:
    
    def __init__(self): pass

    def timeit(self, method):
        def timed(*args, **kw):
            ts = time.time()
            result = method(*args, **kw)
            te = time.time()
            return round(te-ts, 2), result
        return timed class Clusterer:

    def __init__(self):
        self.clust = Deduplicator()
        self.seqio = SequenceIO()
        self.seq_library = {}

        self.orfm_regex = OrfM.regular_expression()

    def uncluster_annotations(self, input_annotations, reverse_pipe):
        '''
        Update the annotations hash provided by pplacer to include all
        representatives within each cluster

        Parameters
        ----------
        input_annotations : hash
            Classifications for each representative sequence of the clusters.
            each key being the sequence name, and the entry being the taxonomy
            string as a list.
        reverse_pipe : bool
            True/False, whether the reverse reads pipeline is being followed.

        Returns
        -------
        output_annotations : hash
            An updated version of the above, which includes all reads from
            each cluster
        '''
        output_annotations = {}
        for placed_alignment_file_path, clusters in self.seq_library.items():

            if reverse_pipe and placed_alignment_file_path.endswith("_reverse_clustered.fa"): continue
            placed_alignment_file = os.path.basename(placed_alignment_file_path)
            cluster_classifications = input_annotations[placed_alignment_file]

            if reverse_pipe:
                placed_alignment_base = placed_alignment_file.replace('_forward_clustered.fa', '')
            else:
                placed_alignment_base = placed_alignment_file.replace('_clustered.fa', '')
            output_annotations[placed_alignment_base] = {}
            for rep_read_name, rep_read_taxonomy in cluster_classifications.items():

                if reverse_pipe:
                    orfm_regex = OrfM.regular_expression()
                    clusters={(orfm_regex.match(key).groups(0)[0] if orfm_regex.match(key) else key): item for key, item in iter(clusters.items())}
                for read in clusters[rep_read_name]:
                    output_annotations[placed_alignment_base][read.name] = rep_read_taxonomy

        return output_annotations

    def cluster(self, input_fasta_list, reverse_pipe):
        '''
        cluster - Clusters reads at 100% identity level and  writes them to
        file. Resets the input_fasta variable as the FASTA file containing the
        clusters.

        Parameters
        ----------
        input_fasta_list : list
            list of strings, each a path to input fasta files to be clustered.
        reverse_pipe : bool
            True/False, whether the reverse reads pipeline is being followed.
        Returns
        -------
        output_fasta_list : list
            list of strings, each a path to the output fasta file to which
            clusters were written to.
        '''
        output_fasta_list = []
        for input_fasta in input_fasta_list:
            output_path  = input_fasta.replace('_hits.aln.fa', '_clustered.fa')
            cluster_dict = {}

            logging.debug('Clustering reads')
            if os.path.exists(input_fasta):
                reads=self.seqio.read_fasta_file(input_fasta) # Read in FASTA records
                logging.debug('Found %i reads' % len(reads)) # Report number found
                clusters=self.clust.deduplicate(reads) # Cluster redundant sequences
                logging.debug('Clustered to %s groups' % len(clusters)) # Report number of clusters
                logging.debug('Writing representative sequences of each cluster to: %s' % output_path) # Report the name of the file
            else:
                logging.debug("Found no reads to be clustered")
                clusters = []

            self.seqio.write_fasta_file(
                                        [x[0] for x in clusters],
                                        output_path
                                        ) # Choose the first sequence to write to file as representative (all the same anyway)
            for cluster in clusters:
                cluster_dict[cluster[0].name]=cluster # assign the cluster to the dictionary
            self.seq_library[output_path]= cluster_dict

            output_fasta_list.append(output_path)

        return output_fasta_list 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 ExternalProgramSuite:

    _FASTTREE="FastTreeMP"

    package_urls = {'taxit': 'https://github.com/fhcrc/taxtastic',
             'FastTree': 'http://www.microbesonline.org/fasttree/',
             'orfm': 'https://github.com/wwood/OrfM',
            'mfqe': 'https://github.com/wwood/mfqe',
            'pplacer': 'http://matsen.fhcrc.org/pplacer/',
            'krona': 'http://sourceforge.net/p/krona/home/krona/',
            'mafft': 'http://mafft.cbrc.jp/alignment/software/',
            'diamond': 'https://github.com/bbuchfink/diamond',
            'hmmer': 'http://hmmer.janelia.org/'}

    programs_to_packages = {'hmmalign': 'hmmer',
                            'hmmsearch': 'hmmer',
                            'nhmmer': 'hmmer',
                            'ktImportText': 'krona',
                            'FastTreeMP': 'FastTree'}

    programs_to_possibilities = {
                                 _FASTTREE: ["FastTreeMP",
                                              "fasttree",
                                              "fasttreeMP",
                                              "FastTree"]
                                 }

    def __init__(self, program_list):
        '''Given a list of executable names, check that they are available
        on the PATH, raising an exception otherwise

        Parameters
        ----------
        program_list: iterable of str
            list of program names to check
        '''

        uninstalled_programs = []
        for program in program_list:
            if program in ExternalProgramSuite.programs_to_possibilities:
                check = [p for p in ExternalProgramSuite\
                                        .programs_to_possibilities[program]
                         if extern.which(p)]
                if len(check)>1:
                    logging.warning("Program found with multiple commands: %s. \
Arbitrarily selecting %s" % (' '.join(program_list), check[0]))

                for key, item in ExternalProgramSuite\
                                        .programs_to_possibilities.items():
                    if program in item:
                        program = key

                if len(check) > 0:
                    if key == self._FASTTREE:
                        self.fasttree = check[0]
                else:
                    uninstalled_programs.append(program)
            else:
                if not extern.which(program):
                    uninstalled_programs.append(program)

        if any(uninstalled_programs):
            msg = "The following programs appear to be missing, and need to be updated or installed before GraftM can continue:"
            logging.error(msg)
            for program in sorted(uninstalled_programs):
                try:
                    package_name = ExternalProgramSuite.programs_to_packages[program]
                    print_name = '/'.join((program, package_name))
                except KeyError:
                    package_name = program
                    print_name = program

                logging.error('\t%25s\t%s' % (print_name, ExternalProgramSuite.package_urls[package_name]))
            exit(1) class Archive:
    def create(self, input_package_path, output_package_path, **kwargs):
        """Create an archived GraftM package

        Parameters
        ----------
        input_package_path: str
            path to gpkg to be archived
        output_pacakge_path: str
            output package path
        kwargs:
            force: bool
                overwrite an existing directory
        """
        force = kwargs.pop('force',ArchiveDefaultOptions.force)
        if len(kwargs) > 0:
            raise Exception("Unexpected arguments detected: %s" % kwargs)

        logging.info("Archiving GraftM package '%s' as '%s'" % (input_package_path, output_package_path))
        gpkg = GraftMPackage.acquire(input_package_path)
        if gpkg.version != 3:
            raise Exception("Archiving GraftM packages only works with format 3 packages")

        self._setup_output(output_package_path, force)

        logging.debug("Compressing contents for archive")
        with tarfile.open(output_package_path, 'w:gz') as tar:
            for path in itertools.chain([gpkg.contents_file_path(),
                                         gpkg.alignment_hmm_path(),
                                         gpkg.reference_package_path(),
                                         gpkg.unaligned_sequence_database_path()],
                                        [hmm for hmm in gpkg.search_hmm_paths() if hmm != gpkg.alignment_hmm_path()]):
                logging.debug("Compressing '%s'" % path)
                # Put the gpkg folder itself in the archive so as not to tar bomb.
                tar.add(path,
                        os.path.join(os.path.basename(os.path.abspath(gpkg._base_directory)),
                                     os.path.basename(path)))
            
        logging.info("Archive successfully created")

    def extract(self, archive_path, output_package_path, **kwargs):
        '''Extract an archived GraftM package.

        Parameters
        ----------
        archive_path: str
            path to archive
        output_package_path: str
            path to where to put the extracted file
        kwargs:
            force: bool
                overwrite an existing directory
        '''
        force = kwargs.pop('force', ArchiveDefaultOptions.force)
        if len(kwargs) > 0:
            raise Exception("Unexpected arguments detected: %s" % kwargs)
        
        logging.info("Un-archiving GraftM package '%s' from '%s'" % (output_package_path, archive_path))
        archive = os.path.abspath(archive_path)
        output = os.path.abspath(output_package_path)
        
        self._setup_output(output, force)
        
        with in_tempdir():
            with tarfile.open(archive) as tar:
                tar.extractall()
                for tarinfo in tar:
                    folder = os.path.dirname(tarinfo.name)
                    break

            # recreate diamond file
            gpkg = GraftMPackage.acquire(folder)
            if gpkg.version != 3:
                raise Exception("Encountered an archived gpkg of unexpected version: %d" % gpkg.version)
            if gpkg.diamond_database_path():
                logging.debug("Creating diamond DB")
                gpkg.create_diamond_db()

            shutil.move(folder, output)

        logging.info("Archive successfully extracted")

    def _setup_output(self, path, force):
        '''Clear the way for an output to be placed at path'''
        # Allow for special case of output being a pipe
        if os.path.isdir(path) or os.path.isfile(path):
            if force:
                logging.warn("Deleting previous file/directory '%s'" % path)
                if os.path.isfile(path):
                    os.remove(path)
                else:
                    shutil.rmtree(path)
            else:
                raise Exception("Cowardly refusing to overwrite already existing path at %s" % path) class DecoyFilter:
    def __init__(self, proper_hits_diamond, decoy_diamond=None):
        '''Create a new DecoyFilter.

        Parameters
        ----------
        proper_hits_diamond: Diamond
            Diamond object containing real sequences.
        decoy_diamond: Diamond or None.
            Diamond object with db containing decoy sequences. If None, no
            searching against a decoy database is carried out.

        '''
        self._decoy_diamond = decoy_diamond
        self._proper_hits_diamond = proper_hits_diamond
        
    def filter(self, candidate_sequences_fasta_path, filtered_output_fasta_path):
        '''Filter the fasta file by only keeping sequences that hit a proper sequence
        better than the decoy database.
            
        Parameters
        ----------
        candidate_sequences_fasta_path: str
            path to candidate sequences fasta file i.e. sequences to be filtered.
        filtered_output_fasta_path: str
            path to file which, after this method is run, will contain only
            those sequences which pass the filter.
            
        Returns
        -------
        False if no sequences remain after filtering, else True.
        '''
        # Run query sequences against the proper database
        seq_ids_and_bitscores = {}
        logging.debug("Running diamond against the non-decoy sequences")
        pd = self._proper_hits_diamond.run(
            candidate_sequences_fasta_path,
            UnpackRawReads.PROTEIN_SEQUENCE_TYPE)
        for res in pd.each([SequenceSearchResult.QUERY_ID_FIELD,
                            SequenceSearchResult.ALIGNMENT_BIT_SCORE]):
            seq = res[0]
            score = float(res[1])
            # Possible a single sequence gets 2 split up hits (maybe), so take
            # the highest bitscore.
            if seq not in seq_ids_and_bitscores or seq_ids_and_bitscores[seq] < score:
                seq_ids_and_bitscores[seq] = score

        num_before_decoy_removal = len(seq_ids_and_bitscores)
        logging.info("Found %i sequences which hit the non-decoy sequences" %\
                     num_before_decoy_removal)

        if self._decoy_diamond is None:
            logging.debug("Not running against the decoy database")
        else:
            # Run the query sequences against the decoy database, removing from the
            # list any sequences which hit better the decoy DB.
            logging.debug("Running diamond against decoy sequences")
            pd = self._decoy_diamond.run(
                candidate_sequences_fasta_path,
                UnpackRawReads.PROTEIN_SEQUENCE_TYPE)
            for res in pd.each([SequenceSearchResult.QUERY_ID_FIELD,
                                SequenceSearchResult.ALIGNMENT_BIT_SCORE]):
                seq = res[0]
                score = float(res[1])
                if seq in seq_ids_and_bitscores and seq_ids_and_bitscores[seq] < score:
                    logging.debug("Removing sequence with better hit to the decoy database: %s" % seq)
                    del seq_ids_and_bitscores[seq]

            logging.info("Removed %i"
                         " sequences which hit the decoy sequences better"
                         " than the non-decoy sequences" %\
                         (num_before_decoy_removal-len(seq_ids_and_bitscores)))
            if len(seq_ids_and_bitscores) == 0:
                return False

        # Extract the found sequences into the output file
        logging.debug("Extracting query sequences")
        SequenceExtractor().extract(seq_ids_and_bitscores.keys(),
                                    candidate_sequences_fasta_path,
                                    filtered_output_fasta_path)
        return True
import os
import logging
import tempfile
import shutil
from graftm.sequence_search_results import SequenceSearchResult
from graftm.graftm_output_paths import GraftMFiles
from graftm.search_table import SearchTableWriter
from graftm.sequence_searcher import SequenceSearcher
from graftm.hmmsearcher import NoInputSequencesException
from graftm.housekeeping import HouseKeeping
from graftm.summarise import Stats_And_Summary
from graftm.pplacer import Pplacer
from graftm.create import Create
from graftm.update import Update
from graftm.unpack_sequences import UnpackRawReads
from graftm.graftm_package import GraftMPackage
from graftm.expand_searcher import ExpandSearcher
from graftm.diamond import Diamond
from graftm.getaxnseq import Getaxnseq
from graftm.sequence_io import SequenceIO
from graftm.timeit import Timer
from graftm.clusterer import Clusterer
from graftm.decorator import Decorator
from graftm.external_program_suite import ExternalProgramSuite
from graftm.archive import Archive
from graftm.decoy_filter import DecoyFilter
from biom.util import biom_open
            else:
                if self.args.search_method==self.hk.HMMSEARCH_SEARCH_METHOD:
                    if not self.args.search_only:
                        maximum_range = self.hk.get_maximum_range(self.args.aln_hmm_file)
                    else:
                        logging.debug("Running search only pipeline. maximum_range not configured.")
                        maximum_range = None
                else:
                    logging.warning('Cannot determine maximum range when using %s pipeline and with no GraftM package specified' % self.args.search_method)
                    logging.warning('Setting maximum_range to None (linked hits will not be detected)')
                    maximum_range = None
            if self.args.search_diamond_file:
                diamond_db = self.args.search_diamond_file
            else:
                if self.args.search_method == self.hk.HMMSEARCH_SEARCH_METHOD:
                    diamond_db = None
                else:
                    logging.error("%s search method selected, but no gpkg or diamond database selected" % self.args.search_method)

        if self.args.assignment_method == Run.DIAMOND_TAXONOMIC_ASSIGNMENT:
            if self.args.reverse:
                logging.warn("--reverse reads specified with --assignment_method diamond. Reverse reads will be ignored.")
                self.args.reverse = None


        # If merge reads is specified, check that there are reverse reads to merge with
        if self.args.merge_reads and not hasattr(self.args, 'reverse'):
            raise Exception("Programming error")

        # Set the output directory if not specified and create that directory
        logging.debug('Creating working directory: %s' % self.args.output_directory)
        self.hk.make_working_directory(self.args.output_directory,
                                       self.args.force)

        # Set pipeline and evalue by checking HMM format
        if self.args.search_only:
            if self.args.search_method == self.hk.HMMSEARCH_SEARCH_METHOD:
                hmm_type, hmm_tc = self.hk.setpipe(self.args.search_hmm_files[0])
                logging.debug("HMM type: %s Trusted Cutoff: %s" % (hmm_type, hmm_tc))
        else:
            hmm_type, hmm_tc = self.hk.setpipe(self.args.aln_hmm_file)
            logging.debug("HMM type: %s Trusted Cutoff: %s" % (hmm_type, hmm_tc))

        if self.args.search_method == self.hk.HMMSEARCH_SEARCH_METHOD:
            setattr(self.args, 'type', hmm_type)
            if hmm_tc:
                setattr(self.args, 'evalue', '--cut_tc')
        else:
            setattr(self.args, 'type', self.PIPELINE_AA)

        if self.args.filter_minimum is not None:
            filter_minimum = self.args.filter_minimum
        else:
            if self.args.type == self.PIPELINE_NT:
                filter_minimum = Run.MIN_ALIGNED_FILTER_FOR_NUCLEOTIDE_PACKAGES
            else:
                filter_minimum = Run.MIN_ALIGNED_FILTER_FOR_AMINO_ACID_PACKAGES

        # Generate expand_search database if required
        if self.args.expand_search_contigs:
            if self.args.graftm_package:
                pkg = GraftMPackage.acquire(self.args.graftm_package)
            else:
                pkg = None
            boots = ExpandSearcher(
                search_hmm_files = self.args.search_hmm_files,
                maximum_range = self.args.maximum_range,
                threads = self.args.threads,
                evalue = self.args.evalue,
                min_orf_length = self.args.min_orf_length,
                graftm_package = pkg)

            # this is a hack, it should really use GraftMFiles but that class isn't currently flexible enough
            new_database = (os.path.join(self.args.output_directory, "expand_search.hmm") \
                            if self.args.search_method == self.hk.HMMSEARCH_SEARCH_METHOD \
                            else os.path.join(self.args.output_directory, "expand_search")
                            )

            if boots.generate_expand_search_database_from_contigs(
                                     self.args.expand_search_contigs,
                                     new_database,
                                     self.args.search_method):
                if self.args.search_method == self.hk.HMMSEARCH_SEARCH_METHOD:
                    self.ss.search_hmm.append(new_database)
                else:
                    diamond_db = new_database

        first_search_method = self.args.search_method
        if self.args.decoy_database:
            decoy_filter = DecoyFilter(Diamond(diamond_db, threads=self.args.threads),
                                       Diamond(self.args.decoy_database,
                                               threads=self.args.threads))
            doing_decoy_search = True
        elif self.args.search_method == self.hk.HMMSEARCH_AND_DIAMOND_SEARCH_METHOD:
            decoy_filter = DecoyFilter(Diamond(diamond_db, threads=self.args.threads))
            doing_decoy_search = True
            first_search_method = self.hk.HMMSEARCH_SEARCH_METHOD
        else:
            doing_decoy_search = False

        # For each pair (or single file passed to GraftM)
        logging.debug('Working with %i file(s)' % len(self.sequence_pair_list))
        for pair in self.sequence_pair_list:
            # Guess the sequence file type, if not already specified to GraftM
            unpack = UnpackRawReads(pair[0],
                                    self.args.input_sequence_type,
                                    INTERLEAVED)

            # Set the basename, and make an entry to the summary table.
            base = unpack.basename()
            pair_direction = ['forward', 'reverse']
            logging.info("Working on %s" % base)

            # Make the working base subdirectory
            self.hk.make_working_directory(os.path.join(self.args.output_directory,
                                                        base),
                                           self.args.force)

            # for each of the paired end read files
            for read_file in pair:

                unpack = UnpackRawReads(read_file,