python-cogent 1.5.1-2 / usr / share / pyshared / cogent / app / infernal.py

#!/usr/bin/env python
"""
Provides an application controller for the commandline version of:
Infernal 1.0 and 1.0.2 only.
"""
from cogent.app.parameters import FlagParameter, ValuedParameter, FilePath
from cogent.app.util import CommandLineApplication, ResultPath, get_tmp_filename
from cogent.parse.fasta import MinimalFastaParser
from cogent.parse.rfam import MinimalRfamParser, ChangedSequence, \
    ChangedRnaSequence, ChangedDnaSequence
from cogent.parse.infernal import CmsearchParser
from cogent.core.moltype import DNA, RNA
from cogent.core.alignment import SequenceCollection, Alignment, DataError
from cogent.format.stockholm import stockholm_from_alignment
from cogent.struct.rna2d import ViennaStructure, wuss_to_vienna
from os import remove

__author__ = "Jeremy Widmann"
__copyright__ = "Copyright 2007-2011, The Cogent Project"
__credits__ = ["Jeremy Widmann"]
__license__ = "GPL"
__version__ = "1.5.1"
__maintainer__ = "Jeremy Widmann"
__email__ = "jeremy.widmann@colorado.edu"
__status__ = "Development"

MOLTYPE_MAP = {'DNA':'--dna',\
                DNA:'--dna',\
               'RNA':'--rna',\
                RNA:'--rna',\
               }
               
SEQ_CONSTRUCTOR_MAP = {'DNA':ChangedDnaSequence,\
                        DNA:ChangedDnaSequence,\
                       'RNA':ChangedRnaSequence,\
                        RNA:ChangedRnaSequence,\
                       }

class Cmalign(CommandLineApplication):
    """cmalign application controller."""
    _options = {
 
    # -o <f> Save the alignment in Stockholm format to a file <f>. The default
    #   is to write it to standard output.
    '-o':ValuedParameter(Prefix='-',Name='o',Delimiter=' '),\
    
    # -l Turn on the local alignment algorithm. Default is global.
    '-l':FlagParameter(Prefix='-',Name='l'),\
    
    # -p Annotate the alignment with posterior probabilities calculated using
    #   the Inside and Outside algorithms.
    '-p':FlagParameter(Prefix='-',Name='p'),\
    
    # -q Quiet; suppress the verbose banner, and only print the resulting
    #   alignment to stdout.
    '-q':FlagParameter(Prefix='-',Name='q'),\
    
    # --informat <s> Assert that the input seqfile is in format <s>. Do not run
    #   Babelfish format autodection. Acceptable formats are: FASTA, EMBL,
    #   UNIPROT, GENBANK, and DDBJ. <s> is case-insensitive.
    '--informat':ValuedParameter(Prefix='--',Name='informat',Delimiter=' '),\
    
    # --mpi Run as an MPI parallel program.  (see User's Guide for details).
    '--mpi':FlagParameter(Prefix='--',Name='mpi'),\
    
    # Expert Options 
    
    # --optacc Align sequences using the Durbin/Holmes optimal accuracy 
    #   algorithm. This is default behavior, so this option is probably useless. 
    '--optacc':FlagParameter(Prefix='--',Name='optacc'),\
    
    # --cyk Do not use the Durbin/Holmes optimal accuracy alignment to align the 
    #   sequences, instead use the CYK algorithm which determines the optimally
    #   scoring alignment of the sequence to the model. 
    '--cyk':FlagParameter(Prefix='--',Name='cyk'),\
    
    # --sample Sample an alignment from the posterior distribution of
    #   alignments.
    '--sample':FlagParameter(Prefix='--',Name='sample'),\
    
    # -s <n> Set the random number generator seed to <n>, where <n> is a 
    #   positive integer. This option can only be used in combination with 
    #   --sample. The default is to use time() to generate a different seed for
    #   each run, which means that two different runs of cmalign --sample on the
    #   same alignment will give slightly different results. You can use this
    #   option to generate reproducible results.
    '-s':ValuedParameter(Prefix='-',Name='s',Delimiter=' '),\
    
    # --viterbi Do not use the CM to align the sequences, instead use the HMM
    #   Viterbi algorithm to align with a CM Plan 9 HMM.
    '--viterbi':FlagParameter(Prefix='--',Name='viterbi'),\
    
    # --sub Turn on the sub model construction and alignment procedure.
    '--sub':FlagParameter(Prefix='--',Name='sub'),\
    
    # --small Use the divide and conquer CYK alignment algorithm described in 
    #   SR Eddy, BMC Bioinformatics 3:18, 2002.
    '--small':FlagParameter(Prefix='--',Name='small'),\
    
    # --hbanded This option is turned on by default. Accelerate alignment by
    #   pruning away regions of the CM DP matrix that are deemed negligible by
    #   an HMM.
    '--hbanded':FlagParameter(Prefix='--',Name='hbanded'),\
    
    # --nonbanded Turns off HMM banding.
    '--nonbanded':FlagParameter(Prefix='--',Name='nonbanded'),\
    
    # --tau <x> Set the tail loss probability used during HMM band calculation
    #   to <x>.
    '--tau':ValuedParameter(Prefix='--',Name='tau',Delimiter=' '),\
    
    # --mxsize <x> Set the maximum allowable DP matrix size to <x> megabytes.
    '--mxsize':ValuedParameter(Prefix='--',Name='mxsize',Delimiter=' '),\

    # --rna Output the alignments as RNA sequence alignments. This is true by
    #   default.
    '--rna':FlagParameter(Prefix='--',Name='rna'),\
    
    # --dna Output the alignments as DNA sequence alignments.
    '--dna':FlagParameter(Prefix='--',Name='dna'),\
    
    # --matchonly Only include match columns in the output alignment, do not
    #   include any insertions relative to the consensus model.
    '--matchonly':FlagParameter(Prefix='--',Name='matchonly'),\
    
    # --resonly Only include match columns in the output alignment that have at
    #   least 1 residue (non-gap character) in them.
    '--resonly':FlagParameter(Prefix='--',Name='resonly'),\
    
    # --fins Change the behavior of how insert emissions are placed in the 
    #   alignment.
    '--fins':FlagParameter(Prefix='--',Name='fins'),\
    
    # --onepost Modifies behavior of the -p option. Use only one character
    #   instead of two to annotate the posterior probability of each aligned
    #   residue.
    '--onepost':FlagParameter(Prefix='--',Name='onepost'),\
    
    # --withali <f> Reads an alignment from file <f> and aligns it as a single
    #   object to the CM; e.g. the alignment in <f> is held fixed.
    '--withali':ValuedParameter(Prefix='--',Name='withali',Delimiter=' '),\
    
    # --withpknots Must be used in combination with --withali <f>. Propogate
    #   structural information for any pseudoknots that exist in <f> to the
    #   output alignment.
    '--withpknots':FlagParameter(Prefix='--',Name='withpknots'),\
    
    # --rf Must be used in combination with --withali <f>. Specify that the
    #   alignment in <f> has the same "#=GC RF" annotation as the alignment file
    #   the CM was built from using cmbuild and further that the --rf option was 
    #   supplied to cmbuild when the CM was constructed.
    '--rf':FlagParameter(Prefix='--',Name='rf'),\
    
    # --gapthresh <x> Must be used in combination with --withali <f>. Specify
    #   that the --gapthresh <x> option was supplied to cmbuild when the CM was
    #   constructed from the alignment file <f>.
    '--gapthresh':ValuedParameter(Prefix='--',Name='gapthresh',Delimiter=' '),\
    
    # --tfile <f> Dump tabular sequence tracebacks for each individual sequence
    #   to a file <f>. Primarily useful for debugging.
    '--tfile':ValuedParameter(Prefix='--',Name='tfile',Delimiter=' '),\
    
    
    }
    _parameters = {}
    _parameters.update(_options)
    _command = "cmalign"
    _suppress_stderr=True

    def getHelp(self):
        """Method that points to the Infernal documentation."""
        
        help_str = \
        """
        See Infernal documentation at:
        http://infernal.janelia.org/
        """
        return help_str
    
    def _tempfile_as_multiline_string(self, data):
        """Write a multiline string to a temp file and return the filename.

            data: a multiline string to be written to a file.

           * Note: the result will be the filename as a FilePath object 
            (which is a string subclass).

        """
        filename = FilePath(self.getTmpFilename(self.TmpDir))
        data_file = open(filename,'w')
        data_file.write(data)
        data_file.close()
        return filename

    def _alignment_out_filename(self):
        
        if self.Parameters['-o'].isOn():
            refined_filename = self._absolute(str(\
                self.Parameters['-o'].Value))
        else:
            raise ValueError, 'No alignment output file specified.'
        return refined_filename

    def _get_result_paths(self,data):
        result = {}
        if self.Parameters['-o'].isOn():
            out_name = self._alignment_out_filename()
            result['Alignment'] = ResultPath(Path=out_name,IsWritten=True)
        
        return result

class Cmbuild(CommandLineApplication):
    """cmbuild application controller."""
    _options = {
    
    # -n <s> Name the covariance model <s>. (Does not work if alifile contains
    #   more than one alignment).
    '-n':ValuedParameter(Prefix='-',Name='n',Delimiter=' '),\
    
    # -A Append the CM to cmfile, if cmfile already exists.
    '-A':FlagParameter(Prefix='-',Name='A'),\
    
    # -F Allow cmfile to be overwritten. Normally, if cmfile already exists,
    #   cmbuild exits with an error unless the -A or -F option is set.
    '-F':FlagParameter(Prefix='-',Name='F'),\
    
    # -v Run in verbose output mode instead of using the default single line
    #   tabular format. This output format is similar to that used by older
    #   versions of Infernal.
    '-v':FlagParameter(Prefix='-',Name='v'),\
    
    # --iins Allow informative insert emissions for the CM. By default, all CM
    #   insert emission scores are set to 0.0 bits.
    '--iins':FlagParameter(Prefix='--',Name='iins'),\
    
    # --Wbeta<x> Set the beta tail loss probability for query-dependent banding
    #   (QDB) to <x> The QDB algorithm is used to determine the maximium length
    #   of a hit to the model. For more information on QDB see (Nawrocki and
    #   Eddy, PLoS Computational Biology 3(3): e56).
    '--Wbeta':ValuedParameter(Prefix='--',Name='Wbeta',Delimiter=' '),\
    
    # Expert Options
    
    # --rsearch <f> Parameterize emission scores a la RSEARCH, using the
    #   RIBOSUM matrix in file <f>. For more information see the RSEARCH 
    #   publication (Klein and Eddy, BMC Bioinformatics 4:44, 2003). Actually,
    #   the emission scores will not exactly With --rsearch enabled, all
    #   alignments in alifile must contain exactly one sequence or the --call
    #   option must also be enabled.
    '--rsearch':ValuedParameter(Prefix='--',Name='rsearch',Delimiter=' '),\
     
    # --binary Save the model in a compact binary format. The default is a more
    #   readable ASCII text format.
    '--binary':FlagParameter(Prefix='--',Name='binary'),\
    
    # --rf Use reference coordinate annotation (#=GC RF line, in Stockholm) to
    #   determine which columns are consensus, and which are inserts.
    '--rf':FlagParameter(Prefix='--',Name='rf'),\
    
    # --gapthresh <x> Set the gap threshold (used for determining which columns
    #   are insertions versus consensus; see --rf above) to <x>. The default is
    #   0.5.
    '--gapthresh':ValuedParameter(Prefix='--',Name='gapthresh',Delimiter=' '),\
    
    # --ignorant Strip all base pair secondary structure information from all
    #   input alignments in alifile before building the CM(s).
    '--ignorant':FlagParameter(Prefix='--',Name='ignorant'),\
    
    # --wgsc Use the Gerstein/Sonnhammer/Chothia (GSC) weighting algorithm.
    #   This is the default unless the number of sequences in the alignment
    #   exceeds a cutoff (see --pbswitch), in which case the default becomes
    #   the faster Henikoff position-based weighting scheme.
    '--wgsc':FlagParameter(Prefix='--',Name='wgsc'),\
    
    # --wblosum Use the BLOSUM filtering algorithm to weight the sequences,
    #   instead of the default GSC weighting.
    '--wblosum':FlagParameter(Prefix='--',Name='wblosum'),\
    
    # --wpb Use the Henikoff position-based weighting scheme. This weighting
    #   scheme is automatically used (overriding --wgsc and --wblosum) if the
    #   number of sequences in the alignment exceeds a cutoff (see --pbswitch).
    '--wpb':FlagParameter(Prefix='--',Name='wpb'),\
    
    # --wnone Turn sequence weighting off; e.g. explicitly set all sequence
    #   weights to 1.0.
    '--wnone':FlagParameter(Prefix='--',Name='wnone'),\
    
    # --wgiven Use sequence weights as given in annotation in the input
    #   alignment file. If no weights were given, assume they are all 1.0.
    #   The default is to determine new sequence weights by the Gerstein/
    #   Sonnhammer/Chothia algorithm, ignoring any annotated weights.
    '--wgiven':FlagParameter(Prefix='--',Name='wgiven'),\
    
    # --pbswitch <n> Set the cutoff for automatically switching the weighting
    #   method to the Henikoff position-based weighting scheme to <n>. If the
    #   number of sequences in the alignment exceeds <n> Henikoff weighting is
    #   used. By default <n> is 5000.
    '--pbswitch':ValuedParameter(Prefix='--',Name='pbswitch',Delimiter=' '),\
    
    # --wid <x> Controls the behavior of the --wblosum weighting option by
    #   setting the percent identity for clustering the alignment to <x>.
    '--wid':ValuedParameter(Prefix='--',Name='wid',Delimiter=' '),\
    
    # --eent Use the entropy weighting strategy to determine the effective
    #   sequence number that gives a target mean match state relative entropy.
    '--wgiven':FlagParameter(Prefix='--',Name='wgiven'),\
    
    # --enone Turn off the entropy weighting strategy. The effective sequence
    #   number is just the number of sequences in the alignment.
    '--wgiven':FlagParameter(Prefix='--',Name='wgiven'),\
    
    # --ere <x> Set the target mean match state entropy as <x>. By default the
    #   target entropy 1.46 bits.
    '--ere':ValuedParameter(Prefix='--',Name='ere',Delimiter=' '),\
    
    # --null <f> Read a null model from <f>. The null model defines the
    #   probability of each RNA nucleotide in background sequence, the default
    #   is to use 0.25 for each nucleotide.
    '--null':ValuedParameter(Prefix='--',Name='null',Delimiter=' '),\
    
    # --prior <f> Read a Dirichlet prior from <f>, replacing the default mixture 
    #   Dirichlet.
    '--prior':ValuedParameter(Prefix='--',Name='prior',Delimiter=' '),\
    
    # --ctarget <n> Cluster each alignment in alifile by percent identity.
    #   find a cutoff percent id threshold that gives exactly <n> clusters and
    #   build a separate CM from each cluster. If <n> is greater than the number 
    #   of sequences in the alignment the program will not complain, and each
    #   sequence in the alignment will be its own cluster. Each CM will have a
    #   positive integer appended to its name indicating the order in which it
    #   was built.
    '--ctarget':ValuedParameter(Prefix='--',Name='ctarget',Delimiter=' '),\
    
    # --cmaxid <x> Cluster each sequence alignment in alifile by percent
    #   identity. Define clusters at the cutoff fractional id similarity of <x>
    #   and build a separate CM from each cluster.
    '--cmaxid':ValuedParameter(Prefix='--',Name='cmaxid',Delimiter=' '),\
    
    # --call Build a separate CM from each sequence in each alignment in
    #   alifile. Naming of CMs takes place as described above for --ctarget.
    '--call':FlagParameter(Prefix='--',Name='call'),\
    
    # --corig After building multiple CMs using --ctarget, --cmindiff or --call
    #   as described above, build a final CM using the complete original
    #   alignment from alifile.
    '--corig':FlagParameter(Prefix='--',Name='corig'),\
    
    # --cdump<f> Dump the multiple alignments of each cluster to <f> in
    #   Stockholm format. This option only works in combination with --ctarget,
    #   --cmindiff or --call.
    '--cdump':ValuedParameter(Prefix='--',Name='cdump',Delimiter=' '),\
    
    # --refine <f> Attempt to refine the alignment before building the CM using
    #   expectation-maximization (EM). The final alignment (the alignment used
    #   to build the CM that gets written to cmfile) is written to <f>.
    '--refine':ValuedParameter(Prefix='--',Name='refine',Delimiter=' '),\
    
    # --gibbs Modifies the behavior of --refine so Gibbs sampling is used
    #   instead of EM.
    '--gibbs':FlagParameter(Prefix='--',Name='gibbs'),\
    
    # -s <n> Set the random seed to <n>, where <n> is a positive integer.
    #   This option can only be used in combination with --gibbs. The default is 
    #   to use time() to generate a different seed for each run, which means
    #   that two different runs of cmbuild --refine <f> --gibbs on the same
    #   alignment will give slightly different results. You can use this option
    #   to generate reproducible results.
    '-s':ValuedParameter(Prefix='-',Name='s',Delimiter=' '),\
    
    # -l With --refine, turn on the local alignment algorithm, which allows the
    #   alignment to span two or more subsequences if necessary (e.g. if the
    #   structures of the query model and target sequence are only partially
    #   shared), allowing certain large insertions and deletions in the
    #   structure to be penalized differently than normal indels. The default is 
    #   to globally align the query model to the target sequences.
    '-l':ValuedParameter(Prefix='-',Name='l',Delimiter=' '),\
    
    # -a With --refine, print the scores of each individual sequence alignment.
    '-a':ValuedParameter(Prefix='-',Name='a',Delimiter=' '),\
    
    # --cyk With --refine, align with the CYK algorithm.
    '--cyk':FlagParameter(Prefix='--',Name='cyk'),\
    
    # --sub With --refine, turn on the sub model construction and alignment
    #   procedure.
    '--sub':FlagParameter(Prefix='--',Name='sub'),\
    
    # --nonbanded With --refine, do not use HMM bands to accelerate alignment.
    #   Use the full CYK algorithm which is guaranteed to give the optimal
    #   alignment. This will slow down the run significantly, especially for
    #   large models.
    '--nonbanded':FlagParameter(Prefix='--',Name='nonbanded'),\
    
    # --tau <x> With --refine, set the tail loss probability used during HMM
    #   band calculation to <f>. This is the amount of probability mass within
    #   the HMM posterior probabilities that is considered negligible. The
    #   default value is 1E-7. In general, higher values will result in greater
    #   acceleration, but increase the chance of missing the optimal alignment
    #   due to the HMM bands.
    '--tau':ValuedParameter(Prefix='--',Name='tau',Delimiter=' '),\
    
    # --fins With --refine, change the behavior of how insert emissions are
    #   placed in the alignment.
    '--fins':FlagParameter(Prefix='--',Name='fins'),\
    
    # --mxsize <x> With --refine, set the maximum allowable matrix size for
    #   alignment to <x> megabytes.
    '--mxsize':ValuedParameter(Prefix='--',Name='mxsize',Delimiter=' '),\
    
    # --rdump<x> With --refine, output the intermediate alignments at each
    #   iteration of the refinement procedure (as described above for --refine )
    #   to file <f>.
    '--rdump':ValuedParameter(Prefix='--',Name='rdump',Delimiter=' '),\
    
    }
    _parameters = {}
    _parameters.update(_options)
    _command = "cmbuild"
    _suppress_stderr=True
    
    def getHelp(self):
        """Method that points to the Infernal documentation."""
        
        help_str = \
        """
        See Infernal documentation at:
        http://infernal.janelia.org/
        """
        return help_str
    
    def _refine_out_filename(self):
        
        if self.Parameters['--refine'].isOn():
            refined_filename = self._absolute(str(\
                self.Parameters['--refine'].Value))
        else:
            raise ValueError, 'No refine output file specified.'
        return refined_filename
    
    def _cm_out_filename(self):
        
        if self.Parameters['-n'].isOn():
            refined_filename = self._absolute(str(\
                self.Parameters['-n'].Value))
        else:
            raise ValueError, 'No cm output file specified.'
        return refined_filename
    
    def _tempfile_as_multiline_string(self, data):
        """Write a multiline string to a temp file and return the filename.

            data: a multiline string to be written to a file.

           * Note: the result will be the filename as a FilePath object 
            (which is a string subclass).

        """
        filename = FilePath(self.getTmpFilename(self.TmpDir))
        data_file = open(filename,'w')
        data_file.write(data)
        data_file.close()
        return filename
    
    def _get_result_paths(self,data):
        result = {}
        if self.Parameters['--refine'].isOn():
            out_name = self._refine_out_filename()
            result['Refined'] = ResultPath(Path=out_name,IsWritten=True)
        if self.Parameters['-n'].isOn():
            cm_name = self._cm_out_filename()
            result['CmFile'] = ResultPath(Path=cm_name,IsWritten=True)
        
        return result

    
class Cmcalibrate(CommandLineApplication):
    """cmcalibrate application controller."""
    _options = {
    
    # -s <n> Set the random number generator seed to <n>, where <n> is a
    #   positive integer. The default is to use time() to generate a different
    #   seed for each run, which means that two different runs of cmcalibrate on 
    #   the same CM will give slightly different E-value and HMM filter
    #   threshold parameters. You can use this option to generate reproducible
    #   results.
    '-s':ValuedParameter(Prefix='-',Name='s',Delimiter=' '),\
    
    # --forecast <n> Predict the running time of the calibration for cmfile and
    #   provided options and exit, DO NOT perform the calibration.
    '--forecast':ValuedParameter(Prefix='--',Name='forecast',Delimiter=' '),\
    
    # --mpi Run as an MPI parallel program.
    '--mpi':FlagParameter(Prefix='--',Name='mpi'),\
    
    # Expert Options
    
    # --exp-cmL-glc <x> Set the length of random sequence to search for the CM
    #   glocal exponential tail fits to <x> megabases (Mb).
    '--exp-cmL-glc':ValuedParameter(Prefix='--',Name='exp-cmL-glc',\
        Delimiter=' '),\
    
    # --exp-cmL-loc <x> Set the length of random sequence to search for the CM
    #   local exponential tail fits to <x> megabases (Mb).
    '--exp-cmL-loc':ValuedParameter(Prefix='--',Name='exp-cmL-loc',\
        Delimiter=' '),\
    
    # --exp-hmmLn-glc <x> Set the minimum random sequence length to search for
    #   the HMM glocal exponential tail fits to <x> megabases (Mb).
    '--exp-hmmLn-glc':ValuedParameter(Prefix='--',Name='exp-hmmLn-glc',\
        Delimiter=' '),\
    
    # --exp-hmmLn-loc <x> Set the minimum random sequence length to search for
    #   the HMM local exponential tail fits to <x> megabases (Mb).
    '--exp-hmmLn-loc':ValuedParameter(Prefix='--',Name='exp-hmmLn-loc',\
        Delimiter=' '),\
    
    # --exp-hmmLx <x> Set the maximum random sequence length to search when
    #   determining HMM E-values to <x> megabases (Mb).
    '--exp-hmmLx':ValuedParameter(Prefix='--',Name='exp-hmmLx',Delimiter=' '),\
    
    # --exp-fract <x> Set the HMM/CM fraction of dynamic programming
    #   calculations to <x>.
    '--exp-fract':ValuedParameter(Prefix='--',Name='exp-fract',Delimiter=' '),\
    
    # --exp-tailn-cglc <x> During E-value calibration of glocal CM search modes
    #   fit the exponential tail to the high scores in the histogram tail that
    #   includes <x> hits per Mb searched.
    '--exp-tailn-cglc':ValuedParameter(Prefix='--',Name='exp-tailn-cglc',\
        Delimiter=' '),\
    
    # --exp-tailn-cloc <x> During E-value calibration of local CM search modes
    #   fit the exponential tail to the high scores in the histogram tail that
    #   includes <x> hits per Mb searched.
    '--exp-tailn-cloc':ValuedParameter(Prefix='--',Name='exp-tailn-cloc',\
        Delimiter=' '),\
    
    # --exp-tailn-hglc <x> During E-value calibration of glocal HMM search modes
    #   fit the exponential tail to the high scores in the histogram tail that
    #   includes <x> hits per Mb searched.
    '--exp-tailn-hglc':ValuedParameter(Prefix='--',Name='exp-tailn-hglc',\
        Delimiter=' '),\
    
    # --exp-tailn-hloc <x> During E-value calibration of local HMM search modes
    #   fit the exponential tail to the high scores in the histogram tail that
    #   includes <x> hits per Mb searched.
    '--exp-tailn-hloc':ValuedParameter(Prefix='--',Name='exp-tailn-hloc',\
        Delimiter=' '),\
    
    # --exp-tailp <x> Ignore the --exp-tailn prefixed options and fit the <x>
    #   fraction right tail of the histogram to exponential tails, for all
    #   search modes.
    '--exp-tailp':ValuedParameter(Prefix='--',Name='exp-tailp',Delimiter=' '),\
    
    # --exp-tailxn <n> With --exp-tailp enforce that the maximum number of hits
    #   in the tail that is fit is <n>.
    '--exp-tailxn':ValuedParameter(Prefix='--',Name='exp-tailxn',\
        Delimiter=' '),\
    
    # --exp-beta <x> During E-value calibration, by default query-dependent
    #   banding (QDB) is used to accelerate the CM search algorithms with a beta
    #   tail loss probability of 1E-15.
    '--exp-beta':ValuedParameter(Prefix='--',Name='exp-beta',Delimiter=' '),\
    
    # --exp-no-qdb Turn of QDB during E-value calibration. This will slow down
    #   calibration, and is not recommended unless you plan on using --no-qdb in
    #   cmsearch.
    '--exp-no-qdb':FlagParameter(Prefix='--',Name='exp-no-qdb'),\
    
    # --exp-hfile <f> Save the histograms fit for the E-value calibration to
    #   file <f>. The format of this file is two tab delimited columns.
    '--exp-hfile':ValuedParameter(Prefix='--',Name='exp-hfile',Delimiter=' '),\
    
    # --exp-sfile <f> Save a survival plot for the E-value calibration to file
    #   <f>. The format of this file is two tab delimited columns.
    '--exp-sfile':ValuedParameter(Prefix='--',Name='exp-sfile',Delimiter=' '),\
    
    # --exp-qqfile <f> Save a quantile-quantile plot for the E-value calibration
    #   to file <f>. The format of this file is two tab delimited columns.
    '--exp-qqfile':ValuedParameter(Prefix='--',Name='exp-qqfile',\
        Delimiter=' '),\
    
    # --exp-ffile <f> Save statistics on the exponential tail statistics to file
    #   <f>. The file will contain the lambda and mu values for exponential
    #   tails fit to tails of different sizes.
    '--exp-ffile':ValuedParameter(Prefix='--',Name='exp-ffile',Delimiter=' '),\
    
    # --fil-N <n> Set the number of sequences sampled and searched for the HMM
    #   filter threshold calibration to <n>. By default, <n> is 10,000.
    '--fil-N':ValuedParameter(Prefix='--',Name='fil-N',Delimiter=' '),\
    
    # --fil-F <x> Set the fraction of sample sequences the HMM filter must be
    #   able to recognize, and allow to survive, to <x>, where <x> is a positive
    #   real number less than or equal to 1.0. By default, <x> is 0.995.
    '--fil-F':ValuedParameter(Prefix='--',Name='fil-F',Delimiter=' '),\
    
    # --fil-xhmm <x> Set the target number of dynamic programming calculations
    #   for a HMM filtered CM QDB search with beta = 1E-7 to <x> times the
    #   number of calculations required to do an HMM search. By default, <x> is
    #   2.0.
    '--fil-xhmm':ValuedParameter(Prefix='--',Name='fil-xhmm',Delimiter=' '),\
    
    # --fil-tau <x> Set the tail loss probability during HMM band calculation
    #   for HMM filter threshold calibration to <x>.
    '--fil-tau':ValuedParameter(Prefix='--',Name='fil-tau',Delimiter=' '),\
    
    # --fil-gemit During HMM filter calibration, always sample sequences from a
    #   globally configured CM, even when calibrating local modes.
    '--fil-gemit':FlagParameter(Prefix='--',Name='fil-gemit'),\
    
    # --fil-dfile <f> Save statistics on filter threshold calibration, including
    #   HMM and CM scores for all sampled sequences, to file <f>.
    '--fil-dfile':ValuedParameter(Prefix='--',Name='fil-dfile',Delimiter=' '),\
    
    # --mxsize <x> Set the maximum allowable DP matrix size to <x> megabytes.
    '--mxsize':ValuedParameter(Prefix='--',Name='mxsize',Delimiter=' '),\

    }
    
    _parameters = {}
    _parameters.update(_options)
    _command = "cmcalibrate"
    _suppress_stderr=True
    
    def getHelp(self):
        """Method that points to the Infernal documentation."""
        
        help_str = \
        """
        See Infernal documentation at:
        http://infernal.janelia.org/
        """
        return help_str

class Cmemit(CommandLineApplication):
    """cmemit application controller."""
    _options = {
    
    # -o <f> Save the synthetic sequences to file <f> rather than writing them
    #   to stdout.
    '-o':ValuedParameter(Prefix='-',Name='o',Delimiter=' '),\
    
    # -n <n> Generate <n> sequences. Default is 10.
    '-n':ValuedParameter(Prefix='-',Name='n',Delimiter=' '),\
    
    # -u Write the generated sequences in unaligned format (FASTA). This is the
    # default, so this option is probably useless.
    '-u':FlagParameter(Prefix='-',Name='u'),\
    
    # -a Write the generated sequences in an aligned format (STOCKHOLM) with
    #   consensus structure annotation rather than FASTA.
    '-a':FlagParameter(Prefix='-',Name='a'),\
    
    # -c Predict a single majority-rule consensus sequence instead of sampling
    #   sequences from the CM's probability distribution.
    '-c':FlagParameter(Prefix='-',Name='c'),\
    
    # -l Configure the CMs into local mode before emitting sequences. See the
    #   User's Guide for more information on locally configured CMs.
    '-l':FlagParameter(Prefix='-',Name='l'),\
    
    # -s <n> Set the random seed to <n>, where <n> is a positive integer. The
    #   default is to use time() to generate a different seed for each run,
    #   which means that two different runs of cmemit on the same CM will give
    #   different results. You can use this option to generate reproducible
    #   results.
    '-s':ValuedParameter(Prefix='-',Name='s',Delimiter=' '),\
    
    # --rna Specify that the emitted sequences be output as RNA sequences. This
    #   is true by default.
    '--rna':FlagParameter(Prefix='--',Name='rna'),\
    
    # --dna Specify that the emitted sequences be output as DNA sequences. By
    #   default, the output alphabet is RNA.
    '--dna':FlagParameter(Prefix='--',Name='dna'),\
    
    # --tfile <f> Dump tabular sequence parsetrees (tracebacks) for each emitted
    #   sequence to file <f>. Primarily useful for debugging.
    '--tfile':ValuedParameter(Prefix='--',Name='tfile',Delimiter=' '),\
    
    # --exp <x> Exponentiate the emission and transition probabilities of the CM
    #   by <x> and then renormalize those distributions before emitting
    #   sequences.
    '--exp':ValuedParameter(Prefix='--',Name='exp',Delimiter=' '),\
    
    # --begin <n> Truncate the resulting alignment by removing all residues
    #   before consensus column <n>, where <n> is a positive integer no greater
    #   than the consensus length of the CM. Must be used in combination with
    #   --end and either -a or --shmm (a developer option).
    '--begin':ValuedParameter(Prefix='--',Name='begin',Delimiter=' '),\
    
    # --end <n> Truncate the resulting alignment by removing all residues after
    #   consensus column <n>, where <n> is a positive integer no greater than
    #   the consensus length of the CM. Must be used in combination with --begin
    #   and either -a or --shmm (a developer option).
    '--end':ValuedParameter(Prefix='--',Name='end',Delimiter=' '),\
    
    }
    _parameters = {}
    _parameters.update(_options)
    _command = "cmemit"
    _suppress_stderr=True
    
    def getHelp(self):
        """Method that points to the Infernal documentation."""
        
        help_str = \
        """
        See Infernal documentation at:
        http://infernal.janelia.org/
        """
        return help_str
    
class Cmscore(CommandLineApplication):
    """cmscore application controller."""
    _options = {
    
    # -n <n> Set the number of sequences to generate and align to <n>. This
    #   option is incompatible with the --infile option.
    '-n':ValuedParameter(Prefix='-',Name='n',Delimiter=' '),\
    
    # -l Turn on the local alignment algorithm, which allows the alignment to
    #   span two or more subsequences if necessary (e.g. if the structures of
    #   the query model and target sequence are only partially shared), allowing
    #   certain large insertions and deletions in the structure to be penalized
    #   differently than normal indels. The default is to globally align the
    #   query model to the target sequences.
    '-l':FlagParameter(Prefix='-',Name='l'),\
    
    # -s <n> Set the random seed to <n>, where <n> is a positive integer. The
    #   default is to use time() to generate a different seed for each run,
    #   which means that two different runs of cmscore on the same CM will give
    #   different results. You can use this option to generate reproducible
    #   results. The random number generator is used to generate sequences to
    #   score, so -s is incompatible with the --infile option which supplies
    #   the sequences to score in an input file.
    '-s':ValuedParameter(Prefix='-',Name='s',Delimiter=' '),\
    
    # -a Print individual timings and score comparisons for each sequence in
    #   seqfile. By default only summary statistics are printed.
    '-a':FlagParameter(Prefix='-',Name='a'),\
    
    # --sub Turn on the sub model construction and alignment procedure.
    '--sub':FlagParameter(Prefix='--',Name='sub'),\
    
    # --mxsize <x> Set the maximum allowable DP matrix size to <x> megabytes.
    '--mxsize':ValuedParameter(Prefix='--',Name='mxsize',Delimiter=' '),\
    
    # --mpi Run as an MPI parallel program.
    '--mpi':FlagParameter(Prefix='--',Name='mpi'),\
    
    # Expert Options
    
    # --emit Generate sequences to score by sampling from the CM.
    '--emit':FlagParameter(Prefix='--',Name='emit'),\
    
    # --random Generate sequences to score by sampling from the CMs null
    #   distribution. This option turns the --emit option off.
    '--random':FlagParameter(Prefix='--',Name='random'),\
    
    # --infile <f> Sequences to score are read from the file <f>. All the
    #   sequences from <f> are read and scored, the -n and -s options are
    #   incompatible with --infile.
    '--infile':ValuedParameter(Prefix='--',Name='infile',Delimiter=' '),\
    
    # --outfile <f> Save generated sequences that are scored to the file <f> in
    #   FASTA format. This option is incompatible with the --infile option.
    '--outfile':ValuedParameter(Prefix='--',Name='outfile',Delimiter=' '),\
    
    # --Lmin <n1> Must be used in combination with --random and --Lmax <n2>.
    '--Lmin':ValuedParameter(Prefix='--',Name='Lmin',Delimiter=' '),\
    
    # --pad Must be used in combination with --emit and --search. Add <n> cm->W
    #   (max hit length) minus L (sequence <x> length) residues to the 5' and 3'
    #   end of each emitted sequence <x>.
    '--pad':FlagParameter(Prefix='--',Name='pad'),\
    
    # --hbanded Specify that the second stage alignment algorithm be HMM banded
    #   CYK. This option is on by default.
    '--hbanded':FlagParameter(Prefix='--',Name='hbanded'),\
    
    # --tau <x> For stage 2 alignment, set the tail loss probability used during
    #   HMM band calculation to <x>.
    '--tau':ValuedParameter(Prefix='--',Name='tau',Delimiter=' '),\
    
    # --aln2bands With --search, when calculating HMM bands, use an HMM
    #   alignment algorithm instead of an HMM search algorithm.
    '--aln2bands':FlagParameter(Prefix='--',Name='aln2bands'),\
    
    # --hsafe For stage 2 HMM banded alignment, realign any sequences with a
    #   negative alignment score using non-banded CYK to guarantee finding the
    #   optimal alignment.
    '--hsafe':FlagParameter(Prefix='--',Name='hsafe'),\
    
    # --nonbanded Specify that the second stage alignment algorithm be standard,
    #   non-banded, non-D&C CYK. When --nonbanded is enabled, the program fails
    #   with a non-zero exit code and prints an error message if the parsetree
    #   score for any sequence from stage 1 D&C alignment and stage 2 alignment
    #   differs by more than 0.01 bits. In theory, this should never happen as
    #   both algorithms are guaranteed to determine the optimal parsetree. For
    #   larger RNAs (more than 300 residues) if memory is limiting, --nonbanded
    #   should be used in combination with --scoreonly.
    '--nonbanded':FlagParameter(Prefix='--',Name='nonbanded'),\
    
    # --scoreonly With --nonbanded during the second stage standard non-banded
    #   CYK alignment, use the "score only" variant of the algorithm to save
    #   memory, and don't recover a parse tree.
    '--scoreonly':FlagParameter(Prefix='--',Name='scoreonly'),\
    
    # --viterbi Specify that the second stage alignment algorithm be Viterbi to
    #   a CM Plan 9 HMM.
    '--viterbi':FlagParameter(Prefix='--',Name='viterbi'),\
    
    # --search Run all algorithms in scanning mode, not alignment mode.
    '--search':FlagParameter(Prefix='--',Name='search'),\
    
    # --inside With --search Compare the non-banded scanning Inside algorithm to
    #   the HMM banded scanning Inside algorith, instead of using CYK versions.
    '--inside':FlagParameter(Prefix='--',Name='inside'),\
    
    # --forward With --search Compare the scanning Forward scoring algorithm
    #   against CYK.
    '--forward':FlagParameter(Prefix='--',Name='forward'),\
    
    # --taus <n> Specify the first alignment algorithm as non-banded D&C CYK,
    #   and multiple stages of HMM banded CYK alignment. The first HMM banded
    #   alignment will use tau=1E-<x>, which will be the highest value of tau
    #   used. Must be used in combination with --taue.
    '--taus':ValuedParameter(Prefix='--',Name='taus',Delimiter=' '),\
    
    # --taue <n> Specify the first alignment algorithm as non-banded D&C CYK,
    #   and multiple stages of HMM banded CYK alignment. The final HMM banded
    #   alignment will use tau=1E-<x>, which will be the lowest value of tau
    #   used. Must be used in combination with --taus.
    '--taue':ValuedParameter(Prefix='--',Name='taue',Delimiter=' '),\
    
    # --tfile <f> Print the parsetrees for each alignment of each sequence to
    #   file <f>.
    '--tfile':ValuedParameter(Prefix='--',Name='tfile',Delimiter=' '),\

    }
    _parameters = {}
    _parameters.update(_options)
    _command = "cmscore"
    _suppress_stderr=True
    
    def getHelp(self):
        """Method that points to the Infernal documentation."""
        
        help_str = \
        """
        See Infernal documentation at:
        http://infernal.janelia.org/
        """
        return help_str

class Cmsearch(CommandLineApplication):
    """cmsearch application controller."""
    _options = {
    
    # -o <f> Save the high-scoring alignments of hits to a file <f>. The default
    #   is to write them to standard output.
    '-o':ValuedParameter(Prefix='-',Name='o',Delimiter=' '),\
    
    # -g <f> Turn on the 'glocal' alignment algorithm, local with respect to the
    #   target database, and global with respect to the model. By default, the
    #   local alignment algorithm is used which is local with respect to both
    #   the target sequence and the model.
    '-g':ValuedParameter(Prefix='-',Name='g',Delimiter=' '),\
    
    # -p Append posterior probabilities to alignments of hits.
    '-p':FlagParameter(Prefix='-',Name='p'),\
    
    # -x Annotate non-compensatory basepairs and basepairs that include a gap in
    #   the left and/or right half of the pair with x's in the alignments of
    #   hits.
    '-x':FlagParameter(Prefix='-',Name='x'),\
    
    # -Z <x> Calculate E-values as if the target database size was <x> megabases
    #   (Mb). Ignore the actual size of the database. This option is only valid
    #   if the CM file has been calibrated. Warning: the predictions for timings
    #   and survival fractions will be calculated as if the database was of size
    #   <x> Mb, which means they will be inaccurate.
    '-Z':ValuedParameter(Prefix='-',Name='Z',Delimiter=' '),\
    
    # --toponly Only search the top (Watson) strand of the sequences in seqfile.
    #   By default, both strands are searched.
    '--toponly':FlagParameter(Prefix='--',Name='toponly'),\
    
    # --bottomonly Only search the bottom (Crick) strand of the sequences in
    #   seqfile. By default, both strands are searched.
    '--bottomonly':FlagParameter(Prefix='--',Name='bottomonly'),\
    
    # --forecast <n> Predict the running time of the search with provided files
    #   and options and exit, DO NOT perform the search. This option is only
    #   available with calibrated CM files.
    '--forecast':ValuedParameter(Prefix='--',Name='forecast',Delimiter=' '),\
    
    # --informat <s> Assert that the input seqfile is in format <s>. Do not run
    #   Babelfish format autodection. This increases the reliability of the
    #   program somewhat, because the Babelfish can make mistakes; particularly
    #   recommended for unattended, high-throughput runs of @PACKAGE@. <s> is
    #   case-insensitive. Acceptable formats are: FASTA, EMBL, UNIPROT, GENBANK,
    #   and DDBJ. <s> is case-insensitive.
    '--informat':ValuedParameter(Prefix='--',Name='informat',Delimiter=' '),\
    
    # --mxsize <x> Set the maximum allowable DP matrix size to <x> megabytes.
    '--mxsize':ValuedParameter(Prefix='--',Name='mxsize',Delimiter=' '),\
    
    # --mpi Run as an MPI parallel program.
    '--mpi':FlagParameter(Prefix='--',Name='mpi'),\
    
    # Expert Options
    
    # --inside Use the Inside algorithm for the final round of searching. This
    #   is true by default.
    '--inside':FlagParameter(Prefix='--',Name='inside'),\
    
    # --cyk Use the CYK algorithm for the final round of searching.
    '--cyk':FlagParameter(Prefix='--',Name='cyk'),\
    
    # --viterbi Search only with an HMM. This is much faster but less sensitive
    #   than a CM search. Use the Viterbi algorithm for the HMM search.
    '--viterbi':FlagParameter(Prefix='--',Name='viterbi'),\
    
    # --forward Search only with an HMM. This is much faster but less sensitive
    #   than a CM search. Use the Forward algorithm for the HMM search.
    '--forward':FlagParameter(Prefix='--',Name='forward'),\
    
    # -E <x> Set the E-value cutoff for the per-sequence/strand ranked hit list
    #   to <x>, where <x> is a positive real number.
    '-E':ValuedParameter(Prefix='-',Name='E',Delimiter=' '),\
    
    # -T <x> Set the bit score cutoff for the per-sequence ranked hit list to
    #   <x>, where <x> is a positive real number.
    '-T':ValuedParameter(Prefix='-',Name='T',Delimiter=' '),\
    
    # --nc Set the bit score cutoff as the NC cutoff value used by Rfam curators
    #   as the noise cutoff score.
    '--nc':FlagParameter(Prefix='--',Name='nc'),\
    
    # --ga Set the bit score cutoff as the GA cutoff value used by Rfam curators
    #   as the gathering threshold.
    '--ga':FlagParameter(Prefix='--',Name='ga'),\
    
    # --tc Set the bit score cutoff as the TC cutoff value used by Rfam curators
    #   as the trusted cutoff.
    '--tc':FlagParameter(Prefix='--',Name='tc'),\
    
    # --no-qdb Do not use query-dependent banding (QDB) for the final round of
    #   search.
    '--no-qdb':FlagParameter(Prefix='--',Name='no-qdb'),\
    
    # --beta " <x>" For query-dependent banding (QDB) during the final round of
    #   search, set the beta parameter to <x> where <x> is any positive real
    #   number less than 1.0.
    '--beta':ValuedParameter(Prefix='--',Name='beta',Delimiter=' '),\
    
    # --hbanded Use HMM bands to accelerate the final round of search.
    #   Constraints for the CM search are derived from posterior probabilities
    #   from an HMM. This is an experimental option and it is not recommended
    #   for use unless you know exactly what you're doing.
    '--hbanded':FlagParameter(Prefix='--',Name='hbanded'),\
    
    # --tau <x> Set the tail loss probability during HMM band calculation to
    #   <x>.
    '--tau':ValuedParameter(Prefix='--',Name='tau',Delimiter=' '),\
    
    # --fil-no-hmm Turn the HMM filter off.
    '--fil-no-hmm':FlagParameter(Prefix='--',Name='fil-no-hmm'),\
    
    # --fil-no-qdb Turn the QDB filter off.
    '--fil-no-qdb':FlagParameter(Prefix='--',Name='fil-no-qdb'),\
    
    # --fil-beta For the QDB filter, set the beta parameter to <x> where <x> is
    #   any positive real number less than 1.0.
    '--fil-beta':FlagParameter(Prefix='--',Name='fil-beta'),\
    
    # --fil-T-qdb <x> Set the bit score cutoff for the QDB filter round to <x>,
    #   where <x> is a positive real number.
    '--fil-T-qdb':ValuedParameter(Prefix='--',Name='fil-T-qdb',Delimiter=' '),\
    
    # --fil-T-hmm <x> Set the bit score cutoff for the HMM filter round to <x>,
    #   where <x> is a positive real number.
    '--fil-T-hmm':ValuedParameter(Prefix='--',Name='fil-T-hmm',Delimiter=' '),\
    
    # --fil-E-qdb <x> Set the E-value cutoff for the QDB filter round. <x>,
    #   where <x> is a positive real number. Hits with E-values better than
    #   (less than) or equal to this threshold will survive and be passed to the
    #   final round. This option is only available if the CM file has been
    #   calibrated.
    '--fil-E-qdb':ValuedParameter(Prefix='--',Name='fil-E-qdb',Delimiter=' '),\
    
    # --fil-E-hmm <x> Set the E-value cutoff for the HMM filter round. <x>,
    #   where <x> is a positive real number. Hits with E-values better than
    #   (less than) or equal to this threshold will survive and be passed to the
    #   next round, either a QDB filter round, or if the QDB filter is disable,
    #   to the final round of search. This option is only available if the CM
    #   file has been calibrated.
    '--fil-E-hmm':ValuedParameter(Prefix='--',Name='fil-E-hmm',Delimiter=' '),\
    
    # --fil-Smax-hmm <x> Set the maximum predicted survival fraction for an HMM
    #   filter as <x>, where <x> is a positive real number less than 1.0.
    '--fil-Smax-hmm':ValuedParameter(Prefix='--',Name='fil-Smax-hmm',\
        Delimiter=' '),\
    
    # --noalign Do not calculate and print alignments of each hit, only print
    #   locations and scores.
    '--noalign':FlagParameter(Prefix='--',Name='noalign'),\
    
    # --aln-hbanded Use HMM bands to accelerate alignment during the hit
    #   alignment stage.
    '--aln-hbanded':FlagParameter(Prefix='--',Name='aln-hbanded'),\
    
    # --aln-optacc Calculate alignments of hits from final round of search using
    #   the optimal accuracy algorithm which computes the alignment that
    #   maximizes the summed posterior probability of all aligned residues given
    #   the model, which can be different from the highest scoring one.
    '--aln-optacc':FlagParameter(Prefix='--',Name='aln-optacc'),\
    
    # --tabfile <f> Create a new output file <f> and print tabular results to
    #   it.
    '--tabfile':ValuedParameter(Prefix='--',Name='tabfile',Delimiter=' '),\
    
    # --gcfile <f> Create a new output file <f> and print statistics of the GC
    #   content of the sequences in seqfile to it.
    '--gcfile':ValuedParameter(Prefix='--',Name='gcfile',Delimiter=' '),\
    
    # --rna Output the hit alignments as RNA sequences alignments. This is true
    #   by default.
    '--rna':FlagParameter(Prefix='--',Name='rna'),\
    
    # --dna Output the hit alignments as DNA sequence alignments.
    '--dna':FlagParameter(Prefix='--',Name='dna'),\

    }
    _parameters = {}
    _parameters.update(_options)
    _command = "cmsearch"
    _suppress_stderr=True
    
    def getHelp(self):
        """Method that points to the Infernal documentation."""
        
        help_str = \
        """
        See Infernal documentation at:
        http://infernal.janelia.org/
        """
        return help_str
    
    def _tabfile_out_filename(self):
        
        if self.Parameters['--tabfile'].isOn():
            tabfile_filename = self._absolute(str(\
                self.Parameters['--tabfile'].Value))
        else:
            raise ValueError, 'No tabfile output file specified.'
        return tabfile_filename
    
    def _tempfile_as_multiline_string(self, data):
        """Write a multiline string to a temp file and return the filename.

            data: a multiline string to be written to a file.

           * Note: the result will be the filename as a FilePath object 
            (which is a string subclass).

        """
        filename = FilePath(self.getTmpFilename(self.TmpDir))
        data_file = open(filename,'w')
        data_file.write(data)
        data_file.close()
        return filename
    
    def _get_result_paths(self,data):
        result = {}
        if self.Parameters['--tabfile'].isOn():
            out_name = self._tabfile_out_filename()
            result['SearchResults'] = ResultPath(Path=out_name,IsWritten=True)
        
        return result

class Cmstat(CommandLineApplication):
    """cmstat application controller."""
    _options = {
    
    # -g Turn on the 'glocal' alignment algorithm, local with respect to the
    #   target database, and global with respect to the model. By default, the
    #   model is configured for local alignment which is local with respect to
    #   both the target sequence and the model.
    '-g':FlagParameter(Prefix='-',Name='g'),\
    
    # -m print general statistics on the models in cmfile and the alignment it
    #   was built from.
    '-m':FlagParameter(Prefix='-',Name='m'),\
    
    # -Z <x> Calculate E-values as if the target database size was <x> megabases
    #   (Mb). Ignore the actual size of the database. This option is only valid
    #   if the CM file has been calibrated.
    '-Z':ValuedParameter(Prefix='-',Name='Z',Delimiter=' '),\
    
    # --all print all available statistics
    '--all':FlagParameter(Prefix='--',Name='all'),\
    
    # --le print local E-value statistics. This option only works if cmfile has
    #   been calibrated with cmcalibrate.
    '--le':FlagParameter(Prefix='--',Name='le'),\
    
    # --ge print glocal E-value statistics. This option only works if cmfile has
    #   been calibrated with cmcalibrate.
    '--ge':FlagParameter(Prefix='--',Name='ge'),\
    
    # --beta <x> With the --search option set the beta parameter for the query-
    #   dependent banding algorithm stages to <x> Beta is the probability mass
    #   considered negligible during band calculation. The default is 1E-7.
    '--beta':ValuedParameter(Prefix='--',Name='beta',Delimiter=' '),\
    
    # --qdbfile <f> Save the query-dependent bands (QDBs) for each state to file
    #   <f>
    '--qdbfile':ValuedParameter(Prefix='--',Name='qdbfile',Delimiter=' '),\
    
    # Expert Options
    
    # --lfi Print the HMM filter thresholds for the range of relevant CM bit
    #   score cutoffs for searches with locally configured models using the
    #   Inside algorithm.
    '--lfi':FlagParameter(Prefix='--',Name='lfi'),\
    
    # --gfi Print the HMM filter thresholds for the range of relevant CM bit
    #   score cutoffs for searches with globally configured models using the
    #   Inside algorithm.
    '--gfi':FlagParameter(Prefix='--',Name='gfi'),\
    
    # --lfc Print the HMM filter thresholds for the range of relevant CM bit
    #   score cutoffs for searches with locally configured models using the CYK
    #   algorithm.
    '--lfc':FlagParameter(Prefix='--',Name='lfc'),\
    
    # --gfc Print the HMM filter thresholds for the range of relevant CM bit
    #   score cutoffs for searches with globally configured models using the CYK
    #   algorithm.
    '--gfc':FlagParameter(Prefix='--',Name='gfc'),\
    
    # -E <x> Print filter threshold statistics for an HMM filter if a final CM
    #   E-value cutoff of <x> were to be used for a run of cmsearch on 1 MB of
    #   sequence.
    '-E':ValuedParameter(Prefix='-',Name='E',Delimiter=' '),\
    
    # -T <x> Print filter threshold statistics for an HMM filter if a final CM
    #   bit score cutoff of <x> were to be used for a run of cmsearch.
    '-T':ValuedParameter(Prefix='-',Name='T',Delimiter=' '),\
    
    # --nc Print filter threshold statistics for an HMM filter if a CM bit score
    #   cutoff equal to the Rfam NC cutoff were to be used for a run of
    #   cmsearch.
    '--nc':FlagParameter(Prefix='--',Name='nc'),\
    
    # --ga Print filter threshold statistics for an HMM filter if a CM bit score
    #   cutoff of Rfam GA cutoff value were to be used for a run of cmsearch.
    '--ga':FlagParameter(Prefix='--',Name='ga'),\
    
    # --tc Print filter threshold statistics for an HMM filter if a CM bit score
    #   cutoff equal to the Rfam TC cutoff value were to be used for a run of
    #   cmsearch.
    '--tc':FlagParameter(Prefix='--',Name='tc'),\
    
    # --seqfile <x> With the -E option, use the database size of the database in
    #   <x> instead of the default database size of 1 MB.
    '--seqfile':ValuedParameter(Prefix='--',Name='seqfile',Delimiter=' '),\
    
    # --toponly In combination with --seqfile <x> option, only consider the top
    #   strand of the database in <x> instead of both strands. --search perform
    #   an experiment to determine how fast the CM(s) can search with different
    #   search algorithms.
    '--toponly':FlagParameter(Prefix='--',Name='toponly'),\
    
    # --cmL <n> With the --search option set the length of sequence to search
    #   with CM algorithms as <n> residues. By default, <n> is 1000.
    '--cmL':ValuedParameter(Prefix='--',Name='cmL',Delimiter=' '),\
    
    # --hmmL <n> With the --search option set the length of sequence to search
    #   with HMM algorithms as <n> residues. By default, <n> is 100,000.
    '--hmmL':ValuedParameter(Prefix='--',Name='hmmL',Delimiter=' '),\
    
    # --efile <f> Save a plot of cmsearch HMM filter E value cutoffs versus CM
    #   E-value cutoffs in xmgrace format to file <f>. This option must be used
    #   in combination with --lfi, --gfi, --lfc or --gfc.
    '--efile':ValuedParameter(Prefix='--',Name='efile',Delimiter=' '),\
    
    # --bfile <f> Save a plot of cmsearch HMM bit score cutoffs versus CM bit
    #   score cutoffs in xmgrace format to file <f>. This option must be used in
    #   combination with --lfi, --gfi, --lfc or --gfc.
    '--bfile':ValuedParameter(Prefix='--',Name='bfile',Delimiter=' '),\
    
    # --sfile <f> Save a plot of cmsearch predicted survival fraction from the
    #   HMM filter versus CM E value cutoff in xmgrace format to file <f>. This
    #   option must be used in combination with --lfi, --gfi, --lfc or --gfc.
    '--sfile':ValuedParameter(Prefix='--',Name='sfile',Delimiter=' '),\
    
    # --xfile <f> Save a plot of 'xhmm' versus CM E value cutoff in xmgrace
    #   format to file <f> 'xhmm' is the ratio of the number of dynamic
    #   programming calculations predicted to be required for the HMM filter and
    #   the CM search of the filter survivors versus the number of dynamic
    #   programming calculations for the filter alone. This option must be
    #   used in combination with --lfi, --gfi, --lfc or --gfc.
    '--xfile':ValuedParameter(Prefix='--',Name='xfile',Delimiter=' '),\
    
    # --afile <f> Save a plot of the predicted acceleration for an HMM filtered
    #   search versus CM E value cutoff in xmgrace format to file <f>. This
    #   option must be used in combination with --lfi, --gfi, --lfc or --gfc.
    '--afile':ValuedParameter(Prefix='--',Name='afile',Delimiter=' '),\
    
    # --bits With --efile, --sfile, --xfile, and --afile use CM bit score
    #   cutoffs instead of CM E value cutoffs for the x-axis values of the plot.
    '--bits':FlagParameter(Prefix='--',Name='bits'),\

    }
    _parameters = {}
    _parameters.update(_options)
    _command = "cmstat"
    _suppress_stderr=True
    
    def getHelp(self):
        """Method that points to the Infernal documentation."""
        
        help_str = \
        """
        See Infernal documentation at:
        http://infernal.janelia.org/
        """
        return help_str

def cmbuild_from_alignment(aln, structure_string, refine=False, \
    return_alignment=False,params=None):
    """Uses cmbuild to build a CM file given an alignment and structure string.
    
        - aln: an Alignment object or something that can be used to construct
            one.  All sequences must be the same length.
        - structure_string: vienna structure string representing the consensus
            stucture for the sequences in aln.  Must be the same length as the
            alignment.
        - refine: refine the alignment and realign before building the cm.
            (Default=False)
        - return_alignment: Return (in Stockholm format) alignment file used to
            construct the CM file.  This will either be the original alignment
            and structure string passed in, or the refined alignment if --refine 
            was used. (Default=False)
            - Note.  This will be a string that can either be written to a file
                or parsed.
    """
    aln = Alignment(aln)
    if len(structure_string) != aln.SeqLen:
        raise ValueError, """Structure string is not same length as alignment.  Structure string is %s long. Alignment is %s long."""%(len(structure_string),\
        aln.SeqLen)
    else:
        struct_dict = {'SS_cons':structure_string}
    #Make new Cmbuild app instance.
    app = Cmbuild(InputHandler='_input_as_paths',WorkingDir='/tmp',\
        params=params)
    
    #turn on refine flag if True.
    if refine:
        app.Parameters['--refine'].on(get_tmp_filename(app.WorkingDir))
        
    #Get alignment in Stockholm format
    aln_file_string = stockholm_from_alignment(aln,GC_annotation=struct_dict)
    
    #get path to alignment filename
    aln_path = app._input_as_multiline_string(aln_file_string)
    cm_path = aln_path.split('.txt')[0]+'.cm'
    app.Parameters['-n'].on(cm_path)
    
    filepaths = [cm_path,aln_path]
    
    res = app(filepaths)
    
    cm_file = res['CmFile'].read()
    
    if return_alignment:
        #If alignment was refined, return refined alignment and structure,
        # otherwise return original alignment and structure.
        if refine:
            aln_file_string = res['Refined'].read()
        res.cleanUp()
        return cm_file, aln_file_string
    #Just return cm_file
    else:
        res.cleanUp()
        return cm_file


def cmbuild_from_file(stockholm_file_path, refine=False,return_alignment=False,\
    params=None):
    """Uses cmbuild to build a CM file given a stockholm file.
    
        - stockholm_file_path: a path to a stockholm file.  This file should
            contain a multiple sequence alignment formated in Stockholm format. 
            This must contain a sequence structure line:
                #=GC SS_cons <structure string>
        - refine: refine the alignment and realign before building the cm.
            (Default=False)
        - return_alignment: Return alignment and structure string used to
            construct the CM file.  This will either be the original alignment
            and structure string passed in, or the refined alignment if
            --refine was used. (Default=False)
    """
    #get alignment and structure string from stockholm file.
    info, aln, structure_string = \
        list(MinimalRfamParser(open(stockholm_file_path,'U'),\
            seq_constructor=ChangedSequence))[0]
    
    #call cmbuild_from_alignment.
    res = cmbuild_from_alignment(aln, structure_string, refine=refine, \
        return_alignment=return_alignment,params=params)
    return res

def cmalign_from_alignment(aln, structure_string, seqs, moltype,\
    include_aln=True,refine=False, return_stdout=False,params=None,\
    cmbuild_params=None):
    """Uses cmbuild to build a CM file, then cmalign to build an alignment.
    
        - aln: an Alignment object or something that can be used to construct
            one.  All sequences must be the same length.
        - structure_string: vienna structure string representing the consensus
            stucture for the sequences in aln.  Must be the same length as the
            alignment.
        - seqs: SequenceCollection object or something that can be used to
            construct one, containing unaligned sequences that are to be aligned 
            to the aligned sequences in aln.
        - moltype: Cogent moltype object.  Must be RNA or DNA.
        - include_aln: Boolean to include sequences in aln in final alignment.
            (Default=True)
        - refine: refine the alignment and realign before building the cm.
            (Default=False)
        - return_stdout: Boolean to return standard output from infernal.  This
            includes alignment and structure bit scores and average
            probabilities for each sequence. (Default=False)
    """
    #NOTE: Must degap seqs or Infernal well seg fault!
    seqs = SequenceCollection(seqs,MolType=moltype).degap()
    #Create mapping between abbreviated IDs and full IDs
    int_map, int_keys = seqs.getIntMap()
    #Create SequenceCollection from int_map.
    int_map = SequenceCollection(int_map,MolType=moltype)
    
    cm_file, aln_file_string = cmbuild_from_alignment(aln, structure_string,\
        refine=refine,return_alignment=True,params=cmbuild_params)
    
    if params is None:
        params = {}    
    params.update({MOLTYPE_MAP[moltype]:True})
    
    app = Cmalign(InputHandler='_input_as_paths',WorkingDir='/tmp',\
        params=params)
    app.Parameters['--informat'].on('FASTA')
    
    #files to remove that aren't cleaned up by ResultPath object
    to_remove = []    
    #turn on --withali flag if True.
    if include_aln:
        app.Parameters['--withali'].on(\
            app._tempfile_as_multiline_string(aln_file_string))
        #remove this file at end
        to_remove.append(app.Parameters['--withali'].Value)
    
    seqs_path = app._input_as_multiline_string(int_map.toFasta())
    cm_path = app._tempfile_as_multiline_string(cm_file)
    
    #add cm_path to to_remove
    to_remove.append(cm_path)
    paths = [cm_path,seqs_path]

    app.Parameters['-o'].on(get_tmp_filename(app.WorkingDir))
    
    res = app(paths)
    
    info, aligned, struct_string = \
        list(MinimalRfamParser(res['Alignment'].readlines(),\
            seq_constructor=SEQ_CONSTRUCTOR_MAP[moltype]))[0]
    
    #Make new dict mapping original IDs
    new_alignment={}
    for k,v in aligned.NamedSeqs.items():
        new_alignment[int_keys.get(k,k)]=v
    #Create an Alignment object from alignment dict
    new_alignment = Alignment(new_alignment,MolType=moltype)
    
    std_out = res['StdOut'].read()
    #clean up files
    res.cleanUp()
    for f in to_remove: remove(f)
    
    if return_stdout:
        return new_alignment, struct_string, std_out
    else:
        return new_alignment, struct_string
    

def cmalign_from_file(cm_file_path, seqs, moltype, alignment_file_path=None,\
    include_aln=False,return_stdout=False,params=None):
    """Uses cmalign to align seqs to alignment in cm_file_path.
        
        - cm_file_path: path to the file created by cmbuild, containing aligned
            sequences. This will be used to align sequences in seqs.
        - seqs: unaligned sequendes that are to be aligned to the sequences in
            cm_file.
        - moltype: cogent.core.moltype object.  Must be DNA or RNA
        - alignment_file_path: path to stockholm alignment file used to create
            cm_file.
            __IMPORTANT__: This MUST be the same file used by cmbuild
            originally.  Only need to pass in this file if include_aln=True.
            This helper function will NOT check if the alignment file is correct
            so you must use it correctly.
        - include_aln: Boolean to include sequences in aln_file in final
            alignment. (Default=False)
        - return_stdout: Boolean to return standard output from infernal.  This
            includes alignment and structure bit scores and average
            probabilities for each sequence. (Default=False)
    """
    #NOTE: Must degap seqs or Infernal well seg fault!
    seqs = SequenceCollection(seqs,MolType=moltype).degap()
    
    #Create mapping between abbreviated IDs and full IDs
    int_map, int_keys = seqs.getIntMap()
    #Create SequenceCollection from int_map.
    int_map = SequenceCollection(int_map,MolType=moltype)
    
    if params is None:
        params = {}
    params.update({MOLTYPE_MAP[moltype]:True})
    
    app = Cmalign(InputHandler='_input_as_paths',WorkingDir='/tmp',\
        params=params)
    app.Parameters['--informat'].on('FASTA')
        
    #turn on --withali flag if True.
    if include_aln:
        if alignment_file_path is None:
            raise DataError, """Must have path to alignment file used to build CM if include_aln=True."""
        else:
            app.Parameters['--withali'].on(alignment_file_path)
                
    seqs_path = app._input_as_multiline_string(int_map.toFasta())
    paths = [cm_file_path,seqs_path]
    
    app.Parameters['-o'].on(get_tmp_filename(app.WorkingDir))
    res = app(paths)
    
    info, aligned, struct_string = \
        list(MinimalRfamParser(res['Alignment'].readlines(),\
            seq_constructor=SEQ_CONSTRUCTOR_MAP[moltype]))[0]
    
    
    #Make new dict mapping original IDs
    new_alignment={}
    for k,v in aligned.items():
        new_alignment[int_keys.get(k,k)]=v
    #Create an Alignment object from alignment dict
    new_alignment = Alignment(new_alignment,MolType=moltype)
    std_out = res['StdOut'].read()
    res.cleanUp()
    if return_stdout:
        return new_alignment, struct_string, std_out
    else:
        return new_alignment, struct_string
    
def cmsearch_from_alignment(aln, structure_string, seqs, moltype, cutoff=0.0,\
    refine=False,params=None):
    """Uses cmbuild to build a CM file, then cmsearch to find homologs.
    
        - aln: an Alignment object or something that can be used to construct
            one.  All sequences must be the same length.
        - structure_string: vienna structure string representing the consensus
            stucture for the sequences in aln.  Must be the same length as the
            alignment.
        - seqs: SequenceCollection object or something that can be used to
            construct one, containing unaligned sequences that are to be
            searched.
        - moltype: cogent.core.moltype object.  Must be DNA or RNA
        - cutoff: bitscore cutoff.  No sequences < cutoff will be kept in
            search results. (Default=0.0).  Infernal documentation suggests
            a cutoff of log2(number nucleotides searching) will give most
            likely true homologs.
        - refine: refine the alignment and realign before building the cm.
            (Default=False)
    """
    #NOTE: Must degap seqs or Infernal well seg fault!
    seqs = SequenceCollection(seqs,MolType=moltype).degap()
    #Create mapping between abbreviated IDs and full IDs
    int_map, int_keys = seqs.getIntMap()
    #Create SequenceCollection from int_map.
    int_map = SequenceCollection(int_map,MolType=moltype)
    
    cm_file, aln_file_string = cmbuild_from_alignment(aln, structure_string,\
        refine=refine,return_alignment=True)
    
    app = Cmsearch(InputHandler='_input_as_paths',WorkingDir='/tmp',\
        params=params)
    app.Parameters['--informat'].on('FASTA')
    app.Parameters['-T'].on(cutoff)
    
    to_remove = []
    
    seqs_path = app._input_as_multiline_string(int_map.toFasta())
    cm_path = app._tempfile_as_multiline_string(cm_file)
    paths = [cm_path,seqs_path]
    to_remove.append(cm_path)
    
    app.Parameters['--tabfile'].on(get_tmp_filename(app.WorkingDir))
    res = app(paths)
    
    search_results = list(CmsearchParser(res['SearchResults'].readlines()))
    if search_results:
        for i,line in enumerate(search_results):
            label = line[1]
            search_results[i][1]=int_keys.get(label,label)
    
    res.cleanUp()
    for f in to_remove:remove(f)
    
    return search_results

def cmsearch_from_file(cm_file_path, seqs, moltype, cutoff=0.0, params=None):
    """Uses cmbuild to build a CM file, then cmsearch to find homologs.
    
        - cm_file_path: path to the file created by cmbuild, containing aligned
            sequences. This will be used to search sequences in seqs.
        - seqs: SequenceCollection object or something that can be used to
            construct one, containing unaligned sequences that are to be
            searched.
        - moltype: cogent.core.moltype object.  Must be DNA or RNA
        - cutoff: bitscore cutoff.  No sequences < cutoff will be kept in
            search results. (Default=0.0).  Infernal documentation suggests
            a cutoff of log2(number nucleotides searching) will give most
            likely true homologs.
    """
    #NOTE: Must degap seqs or Infernal well seg fault!
    seqs = SequenceCollection(seqs,MolType=moltype).degap()
    #Create mapping between abbreviated IDs and full IDs
    int_map, int_keys = seqs.getIntMap()
    #Create SequenceCollection from int_map.
    int_map = SequenceCollection(int_map,MolType=moltype)
    
    app = Cmsearch(InputHandler='_input_as_paths',WorkingDir='/tmp',\
        params=params)
    app.Parameters['--informat'].on('FASTA')
    app.Parameters['-T'].on(cutoff)
    
    seqs_path = app._input_as_multiline_string(int_map.toFasta())

    paths = [cm_file_path,seqs_path]
    
    app.Parameters['--tabfile'].on(get_tmp_filename(app.WorkingDir))
    res = app(paths)
    
    search_results = list(CmsearchParser(res['SearchResults'].readlines()))
    
    if search_results:    
        for i,line in enumerate(search_results):
            label = line[1]
            search_results[i][1]=int_keys.get(label,label)
    
    res.cleanUp()

    return search_results