def hmmer_search(input, hmmer, output):
    """Blast the fragments against all marker genes+16S sequences, return output
  """
    os.system(
        "%s --noali -E 10000 --cpu %d -o %s %s %s"
        % (options().__getattribute__("hmmsearch").path, options().cpu, output, hmmer, input)
    )

 def generate_backbone(self):
     _LOG.info("Reading input sequences: %s" %(self.options.sequence_file))
     sequences = MutableAlignment()
     sequences.read_file_object(self.options.sequence_file)
     if (options().backbone_size is None):            
         options().backbone_size = min(100,int(.20*sequences.get_num_taxa()))
         _LOG.info("Backbone size set to: %d" %(options().backbone_size))
     backbone_sequences = sequences.get_hard_sub_alignment(random.sample(sequences.keys(), options().backbone_size))        
     [sequences.pop(i) for i in backbone_sequences.keys()]
     
     _LOG.info("Writing query and backbone set. ")
     query = get_temp_file("query", "backbone", ".fas")
     backbone = get_temp_file("backbone", "backbone", ".fas")
     _write_fasta(sequences, query)
     _write_fasta(backbone_sequences, backbone)
             
     _LOG.info("Generating sate backbone alignment and tree. ")
     satealignJob = SateAlignJob()
     moleculeType = options().molecule
     if (options().molecule == 'amino'):
         moleculeType =  'protein'
     satealignJob.setup(backbone,options().backbone_size,self.options.outdir,moleculeType,options().cpu)
     satealignJob.run()
     satealignJob.read_results()
     
     options().placement_size = self.options.backbone_size
     options().alignment_file = open(self.options.outdir + "/sate.fasta")
     options().tree_file = open(self.options.outdir + "/sate.fasttree")
     _LOG.info("Backbone alignment written to %s.\nBackbone tree written to %s" % (options().alignment_file, options().tree_file))
     options().fragment_file = query

def build_profile(input,output_directory):  
  global taxon_map,level_map,key_map,levels
  temp_dir=tempfile.mkdtemp(dir=options().__getattribute__('tempdir'))
  binned_fragments=bin_to_markers(input,temp_dir)    
  
  #load up taxonomy for 30 marker genes
  (taxon_map, level_map, key_map) = load_taxonomy(options().__getattribute__('reference').path + 'refpkg/rpsB.refpkg/all_taxon.taxonomy')
    
  #all classifications stored here  
  classifications = {}

  #Now run TIPP on each fragment    
  for (gene,frags) in binned_fragments.items():    
    #Get size of each marker
    total_taxa = 0
    with open(options().__getattribute__('reference').path + 'refpkg/%s.refpkg/sate.size'%gene, 'r') as f:
      total_taxa = int(f.readline().strip())
    decomp_size = options().alignment_size
    if (decomp_size > total_taxa):
      decomp_size = int(total_taxa/2)
    cpus = options().cpu
    if (len(frags.keys()) < cpus):
      cpus = len(frags.keys())
    os.system('run_tipp.py -c %s --cpu %s -m %s -f %s -t %s -adt %s -a %s -r %s -tx %s -txm %s -at %0.2f -pt %0.2f -A %d -P %d -p %s -o %s -d %s' % (options().config_file.name, cpus, options().molecule, temp_dir+"/%s.frags.fas.fixed" % gene,options().__getattribute__('reference').path + 'refpkg/%s.refpkg/sate.taxonomy'%gene,options().__getattribute__('reference').path + 'refpkg/%s.refpkg/sate.tree'%gene,options().__getattribute__('reference').path + 'refpkg/%s.refpkg/sate.fasta'%gene,options().__getattribute__('reference').path + 'refpkg/%s.refpkg/sate.taxonomy.RAxML_info'%gene,options().__getattribute__('reference').path + 'refpkg/%s.refpkg/all_taxon.taxonomy'%gene,options().__getattribute__('reference').path + 'refpkg/%s.refpkg/species.mapping'%gene,options().alignment_threshold,options().placement_threshold,decomp_size,total_taxa,temp_dir+"/temp_file","tipp_%s" % gene,output_directory+"/markers/"))
    if (not os.path.exists(output_directory+"/markers/tipp_%s_classification.txt" % gene)):
      continue

    gene_classification = generate_classification(output_directory+"/markers/tipp_%s_classification.txt" % gene,0)

    #Now write individual classification and also pool classifications    
    write_classification(gene_classification, output_directory+"/markers/tipp_%s.classification.0" % gene)    
    classifications.update(gene_classification)    
  remove_unclassified_level(classifications)
  write_classification(classifications, output_directory+"/markers/all.classification.0")
  write_abundance(classifications,output_directory)

 def check_options(self, supply=[]):
     if (options().reference_pkg is not None):
         self.load_reference(os.path.join(options().reference.path, 'refpkg/%s.refpkg/' % options().reference_pkg))                  
     if (options().taxonomy_file is None):
         supply = supply + ["taxonomy file"]
     if (options().taxonomy_name_mapping_file is None):
         supply = supply + ["taxonomy name mapping file"]
     ExhaustiveAlgorithm.check_options(self, supply)

def main():
    augment_parser()
    sepp.config._options_singelton = sepp.config._parse_options()
    if (options().alignment_size is None):
        options().alignment_size = 100
    input = options().fragment_file.name
    output_directory = options().outdir
    build_profile(input, output_directory)

    def check_options(self):
        options().info_file = "A_dummy_value"

        if options().tree_file is None or options().alignment_file is None:
            _LOG.error("Specify the backbone alignment and tree and query sequences")
            exit(-1)
        sequences = MutableAlignment()
        sequences.read_file_object(open(self.options.alignment_file.name))
        return ExhaustiveAlgorithm.check_options(self)

 def load_reference(self, reference_pkg):
     file = open(reference_pkg + 'CONTENTS.json')
     result=json.load(file)
     file.close()
     options().taxonomy_name_mapping_file = open(reference_pkg + result['files']['seq_info'])
     options().taxonomy_file = open(reference_pkg + result['files']['taxonomy'])
     options().alignment_file = open(reference_pkg + result['files']['aln_fasta'])
     options().tree_file = open(reference_pkg + result['files']['tree'])
     options().info_file = reference_pkg + result['files']['tree_stats']

def hmmer_to_markers(input, temp_dir):
    global marker_genes
    fragments = MutableAlignment()
    fragments.read_filepath(input)

    reverse = dict([(name+'_rev', reverse_sequence(seq))
                    for (name, seq) in fragments.items()])
    all_frags = MutableAlignment()
    all_frags.set_alignment(fragments)
    all_frags.set_alignment(reverse)
    frag_file = temp_dir+"/frags.fas"
    _write_fasta(all_frags, frag_file)

    # Now bin the fragments
    frag_scores = dict([(name, [-10000, 'NA', 'NA'])
                        for name in fragments.keys()])
    gene_set = marker_genes
    align_name = 'sate'
    if (options().genes == 'cogs'):
        gene_set = cog_genes
        align_name = 'pasta'
    for gene in gene_set:
        # Now run HMMER search
        hmmer_search(
            frag_file,
            os.path.join(
                options().__getattribute__('reference').path,
                'refpkg/%s.refpkg/%.profile' % (gene, align_name)),
            temp_dir + "/%s.out" % gene)
        results = read_hmmsearch_results(temp_dir + "/%s.out" % gene)

        # Now select best direction for each frag
        for name, value in results.items():
            bitscore = value[1]
            direction = 'forward'
            true_name = name
            if (name.find('_rev') != -1):
                true_name = true_name.replace('_rev', '')
                direction = 'reverse'
            if frag_scores[true_name][0] < bitscore:
                frag_scores[true_name] = [bitscore, gene, direction]

    # Now bin the fragments
    genes = dict([])
    for name, val in frag_scores.items():
        if (val[1] not in genes):
            genes[val[1]] = {}
        if (val[2] == 'forward'):
            genes[val[1]][name] = fragments[name]
        else:
            genes[val[1]][name] = reverse_sequence(fragments[name])
    genes.pop("NA", None)
    for gene, seq in genes.items():
        gene_file = temp_dir + "/%s.frags.fas" % gene
        _write_fasta(seq, gene_file + ".fixed")
    return genes

 def __init__(self):
     AbstractAlgorithm.__init__(self)
     self.place_nomatch_fragments = False
     ''' Hardcoded E-Lim for hmmsearch ''' #TODO: what to do with this
     self.elim = 99999999
     self.filters = False
     self.strategy = options().exhaustive.strategy
     self.minsubsetsize = int(options().exhaustive.minsubsetsize)
     #Temp fix for now, 
     self.molecule = self.options.molecule

    def testConfigFile(self):
        # Just to make different test cases independent of each other.
        config._options_singelton = None
        # Diasable main config path for this test
        config.main_config_path = self.fp_config

        sys.argv = [
            sys.argv[0], "-A", "2",
            "-c", get_data_path("configs/test.config"),
            "--outdir", "dir_form_commandline"]

        assert options().alignment_size == 2, \
            "Commandline option -A not read properly"

        assert isinstance(options().config_file, filetypes) and \
            options().config_file.name.endswith("data/configs/test.config"), \
            "Commandline option -c not read properly"

        assert (options().pplacer is not None and
                options().pplacer.path == "pplacer"), \
            "config file options not read properly"

        assert options().placement_size == 10, \
            "Config file option placementSize not read properly"

        assert options().outdir.endswith("dir_form_commandline"), \
            "Config file value outdir is not properly overwritten:%s " % \
            options().outdir

        assert options().tempdir is not None, \
            "Default value not properly set for tempfile attribute"

def hmmer_to_markers(input, temp_dir):
    global marker_genes
    fragments = MutableAlignment()
    fragments.read_filepath(input)

    reverse = dict([(name + "_rev", reverse_sequence(seq)) for (name, seq) in fragments.items()])
    all_frags = MutableAlignment()
    all_frags.set_alignment(fragments)
    all_frags.set_alignment(reverse)
    frag_file = temp_dir + "/frags.fas"
    _write_fasta(all_frags, frag_file)

    # Now bin the fragments
    frag_scores = dict([(name, [-10000, "NA", "NA"]) for name in fragments.keys()])
    gene_set = marker_genes
    align_name = "sate"
    if options().genes == "cogs":
        gene_set = cog_genes
        align_name = "pasta"
    for gene in gene_set:
        # Now run HMMER search
        hmmer_search(
            frag_file,
            os.path.join(
                options().__getattribute__("reference").path, "refpkg/%s.refpkg/%.profile" % (gene, align_name)
            ),
            temp_dir + "/%s.out" % gene,
        )
        results = read_hmmsearch_results(temp_dir + "/%s.out" % gene)

        # Now select best direction for each frag
        for name in results.keys():
            bitscore = results[name][1]
            direction = "forward"
            true_name = name
            if name.find("_rev") != -1:
                true_name = true_name.replace("_rev", "")
                direction = "reverse"
            if frag_scores[true_name][0] < bitscore:
                frag_scores[true_name] = [bitscore, gene, direction]

    # Now bin the fragments
    genes = dict([])
    for name in frag_scores.keys():
        if frag_scores[name][1] not in genes:
            genes[frag_scores[name][1]] = {}
        if frag_scores[name][2] == "forward":
            genes[frag_scores[name][1]][name] = fragments[name]
        else:
            genes[frag_scores[name][1]][name] = reverse_sequence(fragments[name])
    genes.pop("NA", None)
    for gene in genes.keys():
        gene_file = temp_dir + "/%s.frags.fas" % gene
        _write_fasta(genes[gene], gene_file + ".fixed")
    return genes

def blast_fragments(input, output):
    '''Blast the fragments against all marker genes+16S sequences, return
    output'''
    os.system(
        ('%s -db %s -outfmt 6 -query %s -out %s -num_threads %d '
         '-max_target_seqs 1 ') %
        (options().__getattribute__('blast').path,
         os.path.join(
            options().__getattribute__('reference').path,
            "blast/%s/alignment.fasta.db" % options().genes),
         input, output, options().cpu))

 def testCpuCount(self):
     config._options_singelton = None # Just to make different test cases independent of each other.
     back = config.main_config_path
     config.main_config_path = os.path.expanduser("~/.sepp/main.config.notexistentfile") # Diasable main config path for this test
     sys.argv = [sys.argv[0], "-x" ,"7"]
       
     assert options().cpu == 7, "Commandline option -x not read properly"
     
     print options()
     
     config.main_config_path = back        

def blast_fragments(input, output):
    """Blast the fragments against all marker genes+16S sequences, return output
  """
    os.system(
        "%s -db %s -outfmt 6 -query %s -out %s -num_threads %d -max_target_seqs 1 "
        % (
            options().__getattribute__("blast").path,
            os.path.join(options().__getattribute__("reference").path, "blast/%s/alignment.fasta.db" % options().genes),
            input,
            output,
            options().cpu,
        )
    )

 def __init__(self, **kwargs):
     self.job_type = 'jsonmerger'
     ExternalSeppJob.__init__(self, self.job_type, **kwargs)
     self.out_file = None
     self.distribution = False
     self.taxonomy = None
     self.mapping = None
     self.threshold = None
     self.classification_file = None
     self.elim = float(options().hmmsearch.elim)
     if options().hmmsearch.filters.upper() == "TRUE":
         self.filters = True
     else:
         if options().hmmsearch.filters.upper() == "FALSE":
             self.filters = False
         else:
             self.filters = None
     if self.filters is None:
         raise Exception(
             "Expecting true/false for options().hmmsearch.filters")
     self.strategy = options().exhaustive.strategy
     self.minsubsetsize = int(options().exhaustive.minsubsetsize)
     self.alignment_threshold = float(options().alignment_threshold)
     self.molecule = options().molecule
     self.placer = options().exhaustive.__dict__['placer'].lower()
     self.cutoff = 0

    def run(self):
        checkpoint_manager = options().checkpoint
        assert isinstance(checkpoint_manager, CheckPointManager)

        t = time.time()

        if checkpoint_manager.is_recovering:
            checkpoint_manager.restore_checkpoint()
            self.root_problem = \
                checkpoint_manager.checkpoint_state.root_problem
            self.check_outputprefix()
        else:
            '''check input arguments'''
            self.check_options()

            '''build the problem structure'''
            self.root_problem = self.build_subproblems()

            '''build jobs'''
            self.build_jobs()

        '''connect jobs into a DAG'''
        self.connect_jobs()

        '''Queue up first level jobs (i.e. those with no dependency).
        Once these run, they should automatically enqueue the rest of the
        DAG through joins and callbacks '''
        self.enqueue_firstlevel_job()

        '''start the checkpointing (has any effects only in
           checkpointing mode)'''
        checkpoint_manager.start_checkpointing(self.root_problem)

        '''Wait for all jobs to finish'''
        if (not JobPool().wait_for_all_jobs()):
            _LOG.exception(
                "There have been errors in executed jobs. Terminating.")
            sys.exit(1)

        ''' terminate The job pool and release memory'''
        JobPool().terminate()

        ''' Pause Checkpointing'''
        checkpoint_manager.pause_checkpointing()
        # checkpoint_manager.force_checkpoint()

        '''Merge results into final outputs'''
        self.merge_results()

        '''Output final results'''
        self.output_results()

        ''' Pause Checkpointing'''
        checkpoint_manager.stop_checkpointing()

        _LOG.info("Current execution Finished in %d seconds"
                  % (time.time() - t))
        _LOG.info(
            "All checkpointed executions Finished in %d cumulative time" %
            (checkpoint_manager.get_total_time()))

 def __init__(self):
     '''
     Constructor
     '''
     self.root_problem = None
     self.results = None
     self.options = options() # for ease of access
     pass

 def start_checkpointing(self, root_problem):
     if self.is_checkpointing:
         _LOG.info("Checkpoint every %d seconds" %options().checkpoint_interval)
         self.checkpoint_state.root_problem = root_problem 
         self.checkpoint_state.temp_root = get_root_temp_dir()
         if self.checkpoint_state.cumulative_time is None:
             self.checkpoint_state.cumulative_time = 0
         save_checkpoint(self)

def bin_blast_results(input):
    # Map the blast results to the markers
    gene_mapping = read_mapping(
        os.path.join(
            options().__getattribute__('reference').path,
            'blast/%s/seq2marker.tab' % options().genes))

    genes = {}
    with open(input) as f:
        for line in f:
            results = line.split('\t')
            gene = gene_mapping[results[1]][1]
            if gene in genes:
                genes[gene].append(results[0])
            else:
                genes[gene] = [results[0]]
    return genes

 def __init__(self):
     '''
     Constructor
     '''
     self.root_problem = None
     self.results = None
     self.options = options()
     self.outchecked = False  # for ease of access