#! /usr/bin/env python # -*- coding: utf-8 -*- ############################################################################## ## DendroPy Phylogenetic Computing Library. ## ## Copyright 2010-2015 Jeet Sukumaran and Mark T. Holder. ## All rights reserved. ## ## See "LICENSE.rst" for terms and conditions of usage. ## ## If you use this work or any portion thereof in published work, ## please cite it as: ## ## Sukumaran, J. and M. T. Holder. 2010. DendroPy: a Python library ## for phylogenetic computing. Bioinformatics 26: 1569-1571. ## ############################################################################## """ Wrapper around calls to RAxML. """ import sys import os import subprocess import tempfile import dendropy import random from dendropy.utility.messaging import ConsoleMessenger from dendropy.utility import processio def get_messenger(verbosity=1): if verbosity == 0: messaging_level = ConsoleMessenger.ERROR_MESSAGING_LEVEL else: messaging_level = ConsoleMessenger.INFO_MESSAGING_LEVEL messenger = ConsoleMessenger(name="raxml-map-bipartitions", messaging_level=messaging_level) return messenger ############################################################################### ## RAXML WRAPPER # raxmlHPC[-SSE3|-PTHREADS|-PTHREADS-SSE3|-HYBRID|-HYBRID-SSE3] # -s sequenceFileName -n outputFileName -m substitutionModel # [-a weightFileName] [-A secondaryStructureSubstModel] # [-b bootstrapRandomNumberSeed] [-B wcCriterionThreshold] # [-c numberOfCategories] [-C] [-d] [-D] # [-e likelihoodEpsilon] [-E excludeFileName] # [-f a|A|b|c|d|e|E|F|g|h|i|I|j|J|m|n|o|p|r|s|S|t|u|v|w|x|y] [-F] # [-g groupingFileName] [-G placementThreshold] [-h] # [-i initialRearrangementSetting] [-I autoFC|autoMR|autoMRE|autoMRE_IGN] # [-j] [-J MR|MR_DROP|MRE|STRICT|STRICT_DROP] [-k] [-K] [-M] # [-o outGroupName1[,outGroupName2[,...]]] # [-p parsimonyRandomSeed] [-P proteinModel] # [-q multipleModelFileName] [-r binaryConstraintTree] # [-R binaryModelParamFile] [-S secondaryStructureFile] [-t userStartingTree] # [-T numberOfThreads] [-U] [-v] [-w outputDirectory] [-W slidingWindowSize] # [-x rapidBootstrapRandomNumberSeed] [-X] [-y] # [-z multipleTreesFile] [-#|-N numberOfRuns|autoFC|autoMR|autoMRE|autoMRE_IGN] # -a Specify a column weight file name to assign individual weights to each column of # the alignment. Those weights must be integers separated by any type and number # of whitespaces whithin a separate file, see file "example_weights" for an example. # -A Specify one of the secondary structure substitution models implemented in RAxML. # The same nomenclature as in the PHASE manual is used, available models: # S6A, S6B, S6C, S6D, S6E, S7A, S7B, S7C, S7D, S7E, S7F, S16, S16A, S16B # DEFAULT: 16-state GTR model (S16) # -b Specify an integer number (random seed) and turn on bootstrapping # DEFAULT: OFF # -B specify a floating point number between 0.0 and 1.0 that will be used as cutoff threshold # for the MR-based bootstopping criteria. The recommended setting is 0.03. # DEFAULT: 0.03 (recommended empirically determined setting) # -c Specify number of distinct rate catgories for RAxML when modelOfEvolution # is set to GTRCAT or GTRMIX # Individual per-site rates are categorized into numberOfCategories rate # categories to accelerate computations. # DEFAULT: 25 # -C Conduct model parameter optimization on gappy, partitioned multi-gene alignments with per-partition # branch length estimates (-M enabled) using the fast method with pointer meshes described in: # Stamatakis and Ott: "Efficient computation of the phylogenetic likelihood function on multi-gene alignments and multi-core processors" # WARNING: We can not conduct useful tree searches using this method yet! Does not work with Pthreads version. # -d start ML optimization from random starting tree # DEFAULT: OFF # -D ML search convergence criterion. This will break off ML searches if the relative # Robinson-Foulds distance between the trees obtained from two consecutive lazy SPR cycles # is smaller or equal to 1%. Usage recommended for very large datasets in terms of taxa. # On trees with more than 500 taxa this will yield execution time improvements of approximately 50% # While yielding only slightly worse trees. # DEFAULT: OFF # -e set model optimization precision in log likelihood units for final # optimization of tree topology under MIX/MIXI or GAMMA/GAMMAI # DEFAULT: 0.1 for models not using proportion of invariant sites estimate # 0.001 for models using proportion of invariant sites estimate # -E specify an exclude file name, that contains a specification of alignment positions you wish to exclude. # Format is similar to Nexus, the file shall contain entries like "100-200 300-400", to exclude a # single column write, e.g., "100-100", if you use a mixed model, an appropriatly adapted model file # will be written. # -f select algorithm: # "-f a": rapid Bootstrap analysis and search for best-scoring ML tree in one program run # "-f A": compute marginal ancestral states on a ROOTED reference tree provided with "t" # "-f b": draw bipartition information on a tree provided with "-t" based on multiple trees # (e.g., from a bootstrap) in a file specifed by "-z" # "-f c": check if the alignment can be properly read by RAxML # "-f d": new rapid hill-climbing # DEFAULT: ON # "-f e": optimize model+branch lengths for given input tree under GAMMA/GAMMAI only # "-f E": execute very fast experimental tree search, at present only for testing # "-f F": execute fast experimental tree search, at present only for testing # "-f g": compute per site log Likelihoods for one ore more trees passed via # "-z" and write them to a file that can be read by CONSEL # "-f h": compute log likelihood test (SH-test) between best tree passed via "-t" # and a bunch of other trees passed via "-z" # "-f i": EXPERIMENTAL do not use for real tree inferences: conducts a single cycle of fast lazy SPR moves # on a given input tree, to be used in combination with -C and -M # "-f I": EXPERIMENTAL do not use for real tree inferences: conducts a single cycle of thorough lazy SPR moves # on a given input tree, to be used in combination with -C and -M # "-f j": generate a bunch of bootstrapped alignment files from an original alignemnt file. # You need to specify a seed with "-b" and the number of replicates with "-#" # "-f J": Compute SH-like support values on a given tree passed via "-t". # "-f m": compare bipartitions between two bunches of trees passed via "-t" and "-z" # respectively. This will return the Pearson correlation between all bipartitions found # in the two tree files. A file called RAxML_bipartitionFrequencies.outpuFileName # will be printed that contains the pair-wise bipartition frequencies of the two sets # "-f n": compute the log likelihood score of all trees contained in a tree file provided by # "-z" under GAMMA or GAMMA+P-Invar # "-f o": old and slower rapid hill-climbing without heuristic cutoff # "-f p": perform pure stepwise MP addition of new sequences to an incomplete starting tree and exit # "-f r": compute pairwise Robinson-Foulds (RF) distances between all pairs of trees in a tree file passed via "-z" # if the trees have node labales represented as integer support values the program will also compute two flavors of # the weighted Robinson-Foulds (WRF) distance # "-f s": split up a multi-gene partitioned alignment into the respective subalignments # "-f S": compute site-specific placement bias using a leave one out test inspired by the evolutionary placement algorithm # "-f t": do randomized tree searches on one fixed starting tree # "-f u": execute morphological weight calibration using maximum likelihood, this will return a weight vector. # you need to provide a morphological alignment and a reference tree via "-t" # "-f v": classify a bunch of environmental sequences into a reference tree using thorough read insertions # you will need to start RAxML with a non-comprehensive reference tree and an alignment containing all sequences (reference + query) # "-f w": compute ELW test on a bunch of trees passed via "-z" # "-f x": compute pair-wise ML distances, ML model parameters will be estimated on an MP # starting tree or a user-defined tree passed via "-t", only allowed for GAMMA-based # models of rate heterogeneity # "-f y": classify a bunch of environmental sequences into a reference tree using parsimony # you will need to start RAxML with a non-comprehensive reference tree and an alignment containing all sequences (reference + query) # DEFAULT for "-f": new rapid hill climbing # -F enable ML tree searches under CAT model for very large trees without switching to # GAMMA in the end (saves memory). # This option can also be used with the GAMMA models in order to avoid the thorough optimization # of the best-scoring ML tree in the end. # DEFAULT: OFF # -g specify the file name of a multifurcating constraint tree # this tree does not need to be comprehensive, i.e. must not contain all taxa # -G enable the ML-based evolutionary placement algorithm heuristics # by specifiyng a threshold value (fraction of insertion branches to be evaluated # using slow insertions under ML). # -h Display this help message. # -i Initial rearrangement setting for the subsequent application of topological # changes phase # -I a posteriori bootstopping analysis. Use: # "-I autoFC" for the frequency-based criterion # "-I autoMR" for the majority-rule consensus tree criterion # "-I autoMRE" for the extended majority-rule consensus tree criterion # "-I autoMRE_IGN" for metrics similar to MRE, but include bipartitions under the threshold whether they are compatible # or not. This emulates MRE but is faster to compute. # You also need to pass a tree file containg several bootstrap replicates via "-z" # -j Specifies that intermediate tree files shall be written to file during the standard ML and BS tree searches. # DEFAULT: OFF # -J Compute majority rule consensus tree with "-J MR" or extended majority rule consensus tree with "-J MRE" # or strict consensus tree with "-J STRICT". # Options "-J STRICT_DROP" and "-J MR_DROP" will execute an algorithm that identifies dropsets which contain # rogue taxa as proposed by Pattengale et al. in the paper "Uncovering hidden phylogenetic consensus". # You will also need to provide a tree file containing several UNROOTED trees via "-z" # -k Specifies that bootstrapped trees should be printed with branch lengths. # The bootstraps will run a bit longer, because model parameters will be optimized # at the end of each run under GAMMA or GAMMA+P-Invar respectively. # DEFAULT: OFF # -K Specify one of the multi-state substitution models (max 32 states) implemented in RAxML. # Available models are: ORDERED, MK, GTR # DEFAULT: GTR model # -m Model of Binary (Morphological), Nucleotide, Multi-State, or Amino Acid Substitution: # BINARY: # "-m BINCAT" : Optimization of site-specific # evolutionary rates which are categorized into numberOfCategories distinct # rate categories for greater computational efficiency. Final tree might be evaluated # automatically under BINGAMMA, depending on the tree search option # "-m BINCATI" : Optimization of site-specific # evolutionary rates which are categorized into numberOfCategories distinct # rate categories for greater computational efficiency. Final tree might be evaluated # automatically under BINGAMMAI, depending on the tree search option # "-m BINGAMMA" : GAMMA model of rate # heterogeneity (alpha parameter will be estimated) # "-m BINGAMMAI" : Same as BINGAMMA, but with estimate of proportion of invariable sites # NUCLEOTIDES: # "-m GTRCAT" : GTR + Optimization of substitution rates + Optimization of site-specific # evolutionary rates which are categorized into numberOfCategories distinct # rate categories for greater computational efficiency. Final tree might be evaluated # under GTRGAMMA, depending on the tree search option # "-m GTRCATI" : GTR + Optimization of substitution rates + Optimization of site-specific # evolutionary rates which are categorized into numberOfCategories distinct # rate categories for greater computational efficiency. Final tree might be evaluated # under GTRGAMMAI, depending on the tree search option # "-m GTRGAMMA" : GTR + Optimization of substitution rates + GAMMA model of rate # heterogeneity (alpha parameter will be estimated) # "-m GTRGAMMAI" : Same as GTRGAMMA, but with estimate of proportion of invariable sites # MULTI-STATE: # "-m MULTICAT" : Optimization of site-specific # evolutionary rates which are categorized into numberOfCategories distinct # rate categories for greater computational efficiency. Final tree might be evaluated # automatically under MULTIGAMMA, depending on the tree search option # "-m MULTICATI" : Optimization of site-specific # evolutionary rates which are categorized into numberOfCategories distinct # rate categories for greater computational efficiency. Final tree might be evaluated # automatically under MULTIGAMMAI, depending on the tree search option # "-m MULTIGAMMA" : GAMMA model of rate # heterogeneity (alpha parameter will be estimated) # "-m MULTIGAMMAI" : Same as MULTIGAMMA, but with estimate of proportion of invariable sites # You can use up to 32 distinct character states to encode multi-state regions, they must be used in the following order: # 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, A, B, C, D, E, F, G, H, I, J, K, L, M, N, O, P, Q, R, S, T, U, V # i.e., if you have 6 distinct character states you would use 0, 1, 2, 3, 4, 5 to encode these. # The substitution model for the multi-state regions can be selected via the "-K" option # AMINO ACIDS: # "-m PROTCATmatrixName[F]" : specified AA matrix + Optimization of substitution rates + Optimization of site-specific # evolutionary rates which are categorized into numberOfCategories distinct # rate categories for greater computational efficiency. Final tree might be evaluated # automatically under PROTGAMMAmatrixName[f], depending on the tree search option # "-m PROTCATImatrixName[F]" : specified AA matrix + Optimization of substitution rates + Optimization of site-specific # evolutionary rates which are categorized into numberOfCategories distinct # rate categories for greater computational efficiency. Final tree might be evaluated # automatically under PROTGAMMAImatrixName[f], depending on the tree search option # "-m PROTGAMMAmatrixName[F]" : specified AA matrix + Optimization of substitution rates + GAMMA model of rate # heterogeneity (alpha parameter will be estimated) # "-m PROTGAMMAImatrixName[F]" : Same as PROTGAMMAmatrixName[F], but with estimate of proportion of invariable sites # Available AA substitution models: DAYHOFF, DCMUT, JTT, MTREV, WAG, RTREV, CPREV, VT, BLOSUM62, MTMAM, LG, MTART, MTZOA, PMB, HIVB, HIVW, JTTDCMUT, FLU, GTR # With the optional "F" appendix you can specify if you want to use empirical base frequencies # Please note that for mixed models you can in addition specify the per-gene AA model in # the mixed model file (see manual for details). Also note that if you estimate AA GTR parameters on a partitioned # dataset, they will be linked (estimated jointly) across all partitions to avoid over-parametrization # -M Switch on estimation of individual per-partition branch lengths. Only has effect when used in combination with "-q" # Branch lengths for individual partitions will be printed to separate files # A weighted average of the branch lengths is computed by using the respective partition lengths # DEFAULT: OFF # -n Specifies the name of the output file. # -o Specify the name of a single outgrpoup or a comma-separated list of outgroups, eg "-o Rat" # or "-o Rat,Mouse", in case that multiple outgroups are not monophyletic the first name # in the list will be selected as outgroup, don't leave spaces between taxon names! # -p Specify a random number seed for the parsimony inferences. This allows you to reproduce your results # and will help me debug the program. # -P Specify the file name of a user-defined AA (Protein) substitution model. This file must contain # 420 entries, the first 400 being the AA substitution rates (this must be a symmetric matrix) and the # last 20 are the empirical base frequencies # -q Specify the file name which contains the assignment of models to alignment # partitions for multiple models of substitution. For the syntax of this file # please consult the manual. # -r Specify the file name of a binary constraint tree. # this tree does not need to be comprehensive, i.e. must not contain all taxa # -R Specify the file name of a binary model parameter file that has previously been generated # with RAxML using the -f e tree evaluation option. The file name should be: # RAxML_binaryModelParameters.runID # -s Specify the name of the alignment data file in PHYLIP format # -S Specify the name of a secondary structure file. The file can contain "." for # alignment columns that do not form part of a stem and characters "()<>[]{}" to define # stem regions and pseudoknots # -t Specify a user starting tree file name in Newick format # -T PTHREADS VERSION ONLY! Specify the number of threads you want to run. # Make sure to set "-T" to at most the number of CPUs you have on your machine, # otherwise, there will be a huge performance decrease! # -U Try to save memory by using SEV-based implementation for gap columns on large gappy alignments # WARNING: this will only work for DNA under GTRGAMMA and is still in an experimental state. # -v Display version information # -w FULL (!) path to the directory into which RAxML shall write its output files # DEFAULT: current directory # -W Sliding window size for leave-one-out site-specific placement bias algorithm # only effective when used in combination with "-f S" # DEFAULT: 100 sites # -x Specify an integer number (random seed) and turn on rapid bootstrapping # CAUTION: unlike in version 7.0.4 RAxML will conduct rapid BS replicates under # the model of rate heterogeneity you specified via "-m" and not by default under CAT # -X EXPERIMENTAL OPTION: This option will do a per-site estimate of protein substitution models # by looping over all given, fixed models LG, WAG, JTT, etc and using their respective base frequencies to independently # assign a prot subst. model to each site via ML optimization # At present this option only works with the GTR+GAMMA model, unpartitioned datasets, and in the sequential # version only. # DEFAULT: OFF # -y If you want to only compute a parsimony starting tree with RAxML specify "-y", # the program will exit after computation of the starting tree # DEFAULT: OFF # -z Specify the file name of a file containing multiple trees e.g. from a bootstrap # that shall be used to draw bipartition values onto a tree provided with "-t", # It can also be used to compute per site log likelihoods in combination with "-f g" # and to read a bunch of trees for a couple of other options ("-f h", "-f m", "-f n"). # -#|-N Specify the number of alternative runs on distinct starting trees # In combination with the "-b" option, this will invoke a multiple boostrap analysis # Note that "-N" has been added as an alternative since "-#" sometimes caused problems # with certain MPI job submission systems, since "-#" is often used to start comments. # If you want to use the bootstopping criteria specify "-# autoMR" or "-# autoMRE" or "-# autoMRE_IGN" # for the majority-rule tree based criteria (see -I option) or "-# autoFC" for the frequency-based criterion. # Bootstopping will only work in combination with "-x" or "-b" # DEFAULT: 1 single analysis class RaxmlRunner(object): def __init__(self, working_dir_path=None, replace=None, postclean=None, name=None, verbosity=1, raxml_path="raxmlHPC"): self.dirs_to_clean = [] self.files_to_clean = [] if working_dir_path is None: self.working_dir_path = None self.replace = False self.postclean = postclean if postclean is not None else True else: self.working_dir_path = working_dir_path self.replace = replace if replace is not None else True self.postclean = postclean if postclean is not None else False self._name = name self.verbosity = verbosity self.messenger = get_messenger(self.verbosity) self.input_format = "nexus" self.output_format = "nexus" self.raxml_path = raxml_path self.taxon_label_map = {} @property def name(self): if self._name is None: self._name = "dendropy_raxml" return self._name def _compose_fname(self, s): return "RAxML_%s.%s" % (s, self.name) @property def input_seq_fname(self): return "%s.seqs" % self.name @property def best_tree_fname(self): return self._compose_fname("bestTree") @property def bipartitions_fname(self): return self._compose_fname("bipartitions") @property def info_fname(self): return self._compose_fname("info") @property def log_fname(self): return self._compose_fname("log") @property def parsimony_tree_fname(self): return self._compose_fname("parsimonyTree") @property def result_fname(self): return self._compose_fname("result") def _raxml_output_filenames(self): return [self.input_seq_fname, self.best_tree_fname, self.bipartitions_fname, self.info_fname, self.log_fname, self.parsimony_tree_fname, self.result_fname] def _raxml_output_filepaths(self): return [os.path.join(self.working_dir_path, fp) for fp in self._raxml_output_filenames()] def _get_trees(self, tree_filepath, tree_list=None, **kwargs): if tree_list is None: tree_list = dendropy.TreeList() tree_list.read_from_path(tree_filepath, self.input_format, **kwargs) return tree_list def _expand_path(self, path): return os.path.expanduser(os.path.expandvars(path)) def _check_overwrite(self, path): if os.path.exists(path) and not self.replace: ok = input("Overwrite existing file '{}'? (y/n/all [n])? ".format(path)) if not ok: return False ok = ok[0].lower() if ok == "a": self.replace = True return True if ok == "y": return True return False else: return True # def _send_info(self, msg): # self.messenger.send_info(msg, wrap=False) # def _send_warning(self, msg): # self.messenger.send_warning(msg, wrap=False) # def _send_error(self, msg): # self.messenger.send_info(msg, wrap=False) def _write_dummy_seqs(self, taxon_namespace, out): nchar = 19 out.write("{} {}\n".format(len(taxon_namespace), nchar)) bases = ["A", "C", "G", "T"] for idx, taxon in enumerate(taxon_namespace): base_seq = [random.choice(bases) for x in range(nchar)] out.write("{} {}\n".format(taxon.label, "".join(base_seq))) def _remap_taxon_labels(self, taxa): for idx, taxon in enumerate(taxa): label = "T{}".format(idx) self.taxon_label_map[label] = taxon.label taxon.label = label def _create_working_dir(self): if self.working_dir_path is None: self.working_dir_path = tempfile.mkdtemp() self.dirs_to_clean.append(self.working_dir_path) if not os.path.exists(self.working_dir_path): # self._send_info("Creating work directory: {}".format(self.working_dir_path)) os.makedirs(self.working_dir_path) def _clean_working_files(self): for fpath in self._raxml_output_filepaths(): if not self._check_overwrite(fpath): sys.exit(0) if os.path.exists(fpath): os.remove(fpath) def _preclean_working_dir(self): self._clean_working_files() def _postclean_working_dir(self): if self.postclean: # self._send_info("Cleaning up run files") for fpath in self.files_to_clean: # self._send_info("Deleting file: {}".format(fpath)) try: os.remove(fpath) except OSError as e: pass for dir_path in self.dirs_to_clean: # self._send_info("Deleting directory: {}".format(dir_path)) try: os.rmdir(dir_path) except OSError as e: pass def estimate_tree(self, char_matrix, raxml_args=None): # set up taxa taxa = char_matrix.taxon_namespace # create working directory self._create_working_dir() # remap taxon labels self.taxon_label_map = {} self._remap_taxon_labels(taxa) # clean working directory of previous runs self._preclean_working_dir() # write input sequences raxml_seqs_filepath = os.path.join(self.working_dir_path, self.input_seq_fname) # self._send_info("Creating RAxML dummy sequences file: {}".format(raxml_seqs_filepath)) # if not self._check_overwrite(raxml_seqs_filepath): # sys.exit(0) raxml_seqs_filepath_out = open(raxml_seqs_filepath, "w") char_matrix.write_to_stream(raxml_seqs_filepath_out, "phylip") raxml_seqs_filepath_out.flush() raxml_seqs_filepath_out.close() self.files_to_clean.append(raxml_seqs_filepath) self.files_to_clean.append(raxml_seqs_filepath + ".reduced") # run RAxML if raxml_args is None: raxml_args = [] cmd = [self.raxml_path, '-m', 'GTRCAT', '-s', raxml_seqs_filepath, '-n', self.name, '-p', str(random.randint(0, sys.maxsize))] + raxml_args # self._send_info("Executing: {}".format(" ".join(cmd))) if self.verbosity >= 2: stdout_pipe = None stderr_pipe = None else: stdout_pipe = subprocess.PIPE stderr_pipe = subprocess.PIPE p = subprocess.Popen(cmd, stdout=stdout_pipe, stderr=stderr_pipe, cwd=self.working_dir_path) stdout, stderr = processio.communicate(p) if p.returncode != 0: sys.stderr.write("[RAxML run failed]:\n\n%s\n\n" % (" ".join(cmd))) sys.stdout.write(stdout) sys.stderr.write(stderr) sys.exit(p.returncode) # # read result raxml_best_tree_fpath = os.path.join(self.working_dir_path, self.best_tree_fname) if not os.path.exists(raxml_best_tree_fpath): self._send_error("RAxML result not found: {}".format(raxml_best_tree_fpath)) sys.exit(1) best_tree = dendropy.Tree.get_from_path(raxml_best_tree_fpath, "newick", taxon_namespace=taxa) # remap labels for taxon in best_tree.taxon_namespace: taxon.label = self.taxon_label_map[taxon.label] # # write results # mapped_tree.write_to_stream(self.output_dest, self.output_format) # clean-up self._postclean_working_dir() # # return result return best_tree def map_bipartitions(self, target_tree_fpath, bootstrap_trees_fpaths): # set up taxa taxa = dendropy.TaxonNamespace() taxon_label_map = {} # read target tree target_tree_fpath = self._expand_path(target_tree_fpath) # self._send_info("Reading target tree file: {}".format(target_tree_fpath)) target_tree = self._get_trees(target_tree_fpath, taxon_namespace=taxa)[0] # read boostrap trees boot_trees = dendropy.TreeList() for fpath in bootstrap_trees_fpaths: fpath = self._expand_path(fpath) # self._send_info("Reading bootstrap tree file: {}".format(fpath)) self._get_trees(tree_filepath=fpath, tree_list=boot_trees, taxon_namespace=taxa) # self._send_info("Read: {} taxa, {} bootstrap trees".format(len(taxa), len(boot_trees))) # create working directory self._create_working_dir() # remap taxon labels self.taxon_label_map = {} self._remap_taxon_labels(taxa) # write input target tree raxml_target_tree_filepath = os.path.join(self.working_dir_path, "{}.target_tree".format(self.name)) # self._send_info("Creating RAxML target tree file: {}".format(raxml_target_tree_filepath)) if not self._check_overwrite(raxml_target_tree_filepath): sys.exit(0) target_tree.write_to_path(raxml_target_tree_filepath, "newick") self.files_to_clean.append(raxml_target_tree_filepath) # write input bootstrap trees raxml_bootstrap_trees_filepath = os.path.join(self.working_dir_path, "{}.boot_trees".format(self.name)) # self._send_info("Creating RAxML bootstrap tree file: {}".format(raxml_bootstrap_trees_filepath)) if not self._check_overwrite(raxml_bootstrap_trees_filepath): sys.exit(0) boot_trees.write_to_path(raxml_bootstrap_trees_filepath, "newick") self.files_to_clean.append(raxml_bootstrap_trees_filepath) # write input (dummy) sequences raxml_seqs_filepath = os.path.join(self.working_dir_path, "{}.seqs".format(self.name)) # self._send_info("Creating RAxML dummy sequences file: {}".format(raxml_seqs_filepath)) if not self._check_overwrite(raxml_seqs_filepath): sys.exit(0) raxml_seqs_filepath_out = open(raxml_seqs_filepath, "w") self._write_dummy_seqs(taxa, raxml_seqs_filepath_out) raxml_seqs_filepath_out.flush() raxml_seqs_filepath_out.close() self.files_to_clean.append(raxml_seqs_filepath) # clean working directory of previous runs self._preclean_working_dir() # run RAxML cmd = [self.raxml_path, '-f', 'b', '-t', os.path.basename(raxml_target_tree_filepath), '-z', os.path.basename(raxml_bootstrap_trees_filepath), '-s', os.path.basename(raxml_seqs_filepath), '-m', 'GTRCAT', '-n', self.name] # self._send_info("Executing: {}".format(" ".join(cmd))) if self.verbosity >= 2: stdout_pipe = None stderr_pipe = None else: stdout_pipe = subprocess.PIPE stderr_pipe = subprocess.PIPE p = subprocess.Popen(cmd, stdout=stdout_pipe, stderr=stderr_pipe, cwd=self.working_dir_path) stdout, stderr = processio.communicate(p) if p.returncode != 0: self._send_error("RAxML run failed") if self.verbosity < 2: sys.stdout.write(stdout) sys.stderr.write(stderr) sys.exit(p.returncode) # read result raxml_mapped_tree_fpath = os.path.join(self.working_dir_path, self.bipartitions_fname) if not os.path.exists(raxml_mapped_tree_fpath): self._send_error("RAxML result not found: {}".format(raxml_mapped_tree_fpath)) sys.exit(1) mapped_tree = dendropy.Tree.get_from_path(raxml_mapped_tree_fpath, "newick") # remap labels for taxon in mapped_tree.taxon_namespace: taxon.label = taxon_label_map[taxon.label] # # write results # mapped_tree.write_to_stream(self.output_dest, self.output_format) # clean-up self.files_to_clean.append(raxml_mapped_tree_fpath) self.files_to_clean.append(self.info_fname) self._postclean_working_dir() # return result return mapped_tree # def estimate_tree(self, # char_matrix, # raxml_args=None): # command = [self.raxml_path] + self.raxml_args + raxml_args # for dup_cmd in ['-n', '-s']: # if dup_cmd in command: # raise Exception("Cannot specify '%s' option when running through wrapper" % dup_cmd) # # if '-w' not in command: # # work_dir = tempfile.mkdtemp() # # clean_up_dir = True # # command.a # # command.extend(['-w', work_dir]) # # else: # # arg_idx = command(command.index('-w')+1) # # work_dir = os.path.abspath(arg_idx) # # command[arg_idx] = work_dir # # clean_up_dir = False # work_dir = os.path.abspath('.') # input_data_path = os.path.join(work_dir, "input.phy") # char_matrix.write_to_path(input_data_path, "phylip") # command.extend(['-n', 'dendropy.raxml']) # command.extend(['-s', input_data_path]) # run = subprocess.Popen(command, # stdin=subprocess.PIPE, # stdout=subprocess.PIPE, # stderr=subprocess.PIPE) # stdout, stderr = processio.communicate(run) # results = stdout.split("\n") # if run.returncode: # sys.stderr.write("\n*** ERROR FROM RAxML:\n") # sys.stderr.write(stdout) # sys.stderr.write("\n\n*** COMMAND SENT TO RAxML:\n ") # sys.stderr.write(' '.join(command)) # sys.stderr.write("\n") # sys.exit(1) # return results