<!DOCTYPE article  PUBLIC "-//OASIS//DTD DocBook V4.1//EN">
  <article lang="en">
    <articleinfo>
      <title>Web-based SimpleAnalysis HOWTO</title>
      <author>
	<!-- one of (CONTRIB AUTHORBLURB AFFILIATION OTHERNAME LINEAGE
	SURNAME FIRSTNAME HONORIFIC) -->
	<firstname>Richard</firstname>
	<surname>Adams</surname>
	<authorblurb>
	  <para>Bioperl Core Developer</para>
	</authorblurb>
	<affiliation>
	  <orgname>Dept. Medical Genetics, University of Edinburgh</orgname>
	  <address><email>radams_at_staffmail.ed.ac.uk</email></address>
	</affiliation>
      </author>
      <pubdate>2003-11-14</pubdate>
      <revhistory>
	<revision>
		<revnumber>0.1</revnumber>
		<date>2003-11-14</date>
		<authorinitials>RA</authorinitials>
		<revremark>First version</revremark>
	</revision>
       </revhistory>

      <legalnotice>	
	<para>This document is copyright Richard Adams, 2003.  It can
        be copied and distributed under the terms of the Perl
        Artistic License.
	</para>
      </legalnotice>

      <abstract>
	<para>
	  This HOWTO tries to teach you to run a web based sequence
	  analysis program using basic SimpleAnalysis rules.
	</para>
      </abstract>
    </articleinfo>



    <section id="introduction">
      <title>Introduction</title>
      <para>
        There are several situations where it would be useful to run a web
        based sequence analysis program via a Perl script rather than using
        the web interface manually or downloading the program. Firstly, the
        source code or binaries may be unavailable or unavailable on your
        platform. Secondly, the analysis may depend on a regularly updated or
        large database that you don't wish to store or keep updated
        yourself. Thirdly, novel analyses are frequently available via a web
        server before being available for download.
      </para>
      
      <para>
        The aim of the Bio::Tools::Analysis modules is to allow automatic
        submission of sequences to prediction servers to facilitate sequence
        annotation for low to medium numbers of sequences (perhaps tens to
        hundreds). The modules both submit the sequences and parse the results
        into a number of useful formats, usually including Bio::SeqFeature
        objects, Bio::Seq::MetaI sequence objects and a basic Perl data
        structure, as well as the raw output text.
      </para>

      <para>
        At present the following prediction servers are supported, mainly
        reflecting my own research interests. Under current development
        are modules wrapping remote analyses using HMMER, Pfam, the ELM 
        peptide motif database and SIFT for prediction of untolerated 
        non-synonymous SNPs. 
       <table>
        <title>Supported analyses</title>
        <tgroup cols="3">
          <thead>
             <row>
                <entry>Analysis</entry>
                <entry>Purpose</entry>
                <entry>Reference</entry>
             </row>
         </thead>
         <tbody>
              <row>
                 <entry>Domcut</entry>
                 <entry>Protein domain boundaries</entry>
                 <entry>Bioinformatics 19, 673-674 (2003)</entry>
              </row>
              <row>
                 <entry>NetPhos</entry>
                 <entry>Protein phosphorylation sites</entry>
                 <entry>J Mol Biol 294, 1351-1362 (1999)</entry>
              </row>
              <row>
                 <entry>GOR4</entry>
                 <entry>Protein Secondary structure</entry>
                 <entry>Meths. Enzymology 266, 540-553 (1996)</entry>
              </row>
              <row>
                 <entry>HNN</entry>
                 <entry>Protein Secondary structure</entry>
                 <entry>Bioinformatics 15,413-421 (1999)</entry>
              </row>
              <row>
                 <entry>Sopma</entry>
                 <entry>Protein Secondary structure</entry>
                 <entry>Comput Appl Biosci 11, 681-684 (1995)</entry>
              </row>
              <row>
                 <entry>Mitoprot</entry>
                 <entry>Mitochondrial cleavage site prediction</entry>
                 <entry>Eur J Biochem 241, 779-786 (1996)</entry>
              </row>
              <row>
                <entry>ESEfinder</entry>
                <entry>Exonic splice site enhancers</entry>
                <entry>NAR 31, 3568-3571 (2003)</entry>
             </row>
          </tbody>
         </tgroup>
        </table>
      </para>
</section>

<section id="Simple">
      <title>Simple Examples</title>


      <para>
        The script below runs multiple analyses on a single sequence
        and parses the results into standard BioPerl sequence feature
        objects (Bio::SeqFeature::Generic objects to be precise). 
       <programlisting>
     # load up the modules
     
     use Bio::Tools::Analysis::Protein::HNN;
     use Bio::Tools::Analysis::Protein::Domcut;
     use Bio::Tools::Analysis::Protein::MitoProt;
     use Bio::Tools::Analysis::Protein::NetPhos;
     
     
     our @ISA = qw(Bio::Tools::Analysis::SimpleAnalysisBase);
     
     my $seq;# a Bio::Seq object
     for my $method ( qw(Domcut MitoProt NetPhos HNN)) {
     
        #analyses need a Bio::PrimarySeq, not a Bio::Seq;
        my $tool = Bio::Tools::Analysis::Protein::$method->new(
             -seq => $seq->primary_seq);
        $tool->run();
        my @fts = $tool->result('Bio::SeqFeatureI');
        $seq->add_SeqFeature(@fts);	
     }        
       </programlisting>
      </para>
      
      <para>
        The above script runs several analyses using default parameters. All
        analyses have such defaults and  in general only a sequence of the 
        appropriate type is needed for the analysis to be submitted. A sequence
        can be added either in the constructor, as shown above, or by the seq()
        method. 
       <programlisting>
     my $primary_seq = new Bio::PrimarySeq(-seq=>$seq_as_string);
     my $tool =  new Bio::Tools::Analysis::Protein::NetPhos(); 
     $tool->seq($primary_seq);
       </programlisting>
        Note that the only valid sequence format is a Bio::PrimarySeq object.
		 This is in order to support multiple methods of retrieving the 
        results. If you initially have  a Bio::Seq object or Bio::RichSeq
        object (e.g., from a GenBank file) you can call its primary_seq() 
        method to obtain a Bio::PrimarySeq object.
      </para>
</section>


<section id="retrieving_results">
    <title>Retrieving results</title>
    <para>
		If the run() method executes successfully, the raw output (stripped of
        HTML) is now stored in the Analysis object. 
        All modules should return the raw report by a call to result() 
        with no parameters. e.g., 
     <programlisting>
     my $raw_report = $tool->result();
     </programlisting>
        A second way is to retrieve a ready-parsed data structure:
    <programlisting>
     my $parsed_report = $tool->result('parsed'); 
    </programlisting>
        The data structure returned is described in the $RESULT_SPEC->{'raw'}
        hash reference and should always be a native Perl data structure.
   </para>

   <para>
        A third way is to retrieve an array of sequence features:
    <programlisting>
     my @fts = $tool->result('Bio::SeqFeatureI'); 
     $seq->add_SeqFeature(@fts); # add features to sequence.
    </programlisting>

        These are Bio::SequenceFeature::Generic features. Sometimes a
        module might use some code to judge the significance of a
        result prior to listing it as a feature, for example in the
        secondary structure prediction modules.  The rules governing this
        are described in  individual modules. For example, I have put in a rule 
        that only runs of a minimum of 4 consecutive residues predicted to be
        beta sheet or alpha helix can be annotated as features - it makes no sense 
        for a single isolated residue to be annotated as being in a helix. 
		However you might want to filter the raw results yourself, in which 
        case retrieving the results as Bio::Seq::Meta::Array type objects 
        might be better. 

   </para>
   </section>

   <section id="Metasequences">
   <title>Metasequences</title>
     <para>

        Many analyses produce a list of scores, one for each residue
        in a sequence. For example, the protein secondary structure
        prediction program Sopma returns a list of probabilities for
        each residue being in helix, sheet, turn or coil. These
        modules make this data available in the form of meta
        sequences. These are sequence objects derived from
        Bio::PrimarySeq that have arrays of sequence associated data
        contained in the object. The meta sequence names should be
        listed in the individual module documentation.  To retrieve
        results like this supply the string 'meta' as a parameter
        to the result() method.
		Meta sequence objects can access all the PrimarySeq object 
        methods for accessing and manipulating the protein/DNA sequence
        and also have specific methods for accessing the result data. 

      <programlisting>
     $meta_seq = $analysis->result('meta');
      </programlisting>

        This returns a sequence object with the raw analysis data
        accessible through methods e.g.,

     <programlisting>
     my @scores1_20 = $meta_seq->named_sub_meta('Sopma_turn', 1,20);
     </programlisting>

        returns an array of scores for the first 20 amino acids

     <programlisting>
     my @allscores = $meta_seq->named_meta('Sopma_turn');
     </programlisting>

        returns an array of scores for the whole sequence. The names of 
        individual meta sequence names are listed in the module 
        documentation.
     </para>

     <para>

       Although a full description of how metasequence objects work isn't
       really the purpose of this HOWTO there is excellent documentation in
       Bio::Seq::MetaI and Bio::Seq::Meta::Array modules. 
     </para>
</section>



   <section id="parameters">
      <title>How to  run the same analysis with varying parameters</title>

   <para>

      You might want to run some analyses with varying parameters in
      order to determine the effect on the prediction.  
      At present only the Sopma module takes alternative parameters i.e.
      arguments other than just the sequence. Any parameter that is settable
      on the web form should have a method of the same name to get/set its
      values, or alternatively it can be set in the constructor.

      <programlisting>
     my $sopma = Bio::Tools::Analysis::Protein::Sopma->new();
     $sopma->seq(seqobj->primary_seq);
     $sopma->window_width(21);
     </programlisting>

     So, let's suppose we want to investigate how varying the window_width 
     parameter affects the secondary structure prediction for a sequence. 
     We can do this inthe following manner:
   
		
    <programlisting>
     my $seq = Bio::Seq->new(-seq => 'ASFATFDATFATFSTFATSFATFSTAF');
     my $tool = Bio::Tools::Analysis::Protein::Sopma->new
       (-seq=>$seq->primary_seq);
     
     for my $w_size(qw(11 13 15 17 19 21 23)) {
     
         $tool->window_width($w_size);  #set new parameter
         $tool->run(); #default parameters
     
         #2nd param is appended to metasequence name
         $tool->result('meta', $w_size);
             
         ## add to sequence features
         $seq->add_SeqFeature(
         $tool->result('Bio::SeqFeatureI'));
         $tool->clear(); #removes raw data from the previous analysis
     }
     # meta seq now has 28 arrays of metadats - 7 analyses and
     # 4 prediction categories,
     
     ## now get meta_sequence
     my $meta_seq = $tool->seq();
    </programlisting>
        Only 3 new points are raised by this program. Firstly, each time
        the analysis is run it needs to be parsed immediately. The raw
        report data is removed from  one analysis to another to reduce
        memory usage. Secondly, in order to distinguish the meta data 
        each analysis needs to be given a specific identifier which is 
        appended to the default name. So in this case the meta data names 
        would be Sopma_meta|11, Sopma_meta|13, etc. 
        Thirdly, the  clear() method should be called between 
        analyses of the same sequence to ensure the internal data fields
        that hold the raw report are removed. So if you want to keep the raw 
        reports you need to store them after each analysis is returned.
   </para>

   <para>

	    So, how are features obtained from multiple runs distinguished? 
	    This information is contained in tags with the same name as the
        parameters, when the settings aren't the default ones. In other words,
        the features retain knowledge of the analysis method's parameters. 
    
	<programlisting>                         
     ## carrying on from previous listing.... ##
     ##
     ## get all secondary structure features
     my @sec_fts = grep{$_->primary_tag eq '2ary'} 
                       $seq->all_SeqFeatures;
     my @ww11 = grep{($_->each_tag_value('method'))[0]  eq 'Sopma' &amp;&amp;
                     ($_->each_tag_value('window_width'))[0] == 11}
                     @sec_fts;
     
     ## now onto comparison..... ##
     </programlisting> 


     </para>
</section>

   <section id="caveats">
      <title>Caveats</title>
      <para>

        The problem with these analyses is their speed of execution, which
        obviously depends on your network speed, the complexity of the
        analysis and the speed of the host. Moreover, it is usually polite to
        leave a second or two between requests to avoid blocking the server,
        which means that to annotate a single protein sequence with all of the
        X methods available may take over a minute. Certainly these modules
        are unsuitable for genome scale analysis and are designed for use with
        smaller numbers of sequences. However, the code in the result()
        methods should still be usable for parsing analyses run locally when
        the local output and webpage output are the same format.

      </para>


   </section>


   <section id="developing">
      <title>Interested in developing your own Analysis module?</title>
      <para>
        Most of the hard work is done by Bio::WebAgent and
        Bio::SimpleAnalysisBase. Any module must inherit from this latter
        module as it contains methods common to all analyses.
      </para>

      <para>
        For a minimal prediction server which takes just sequence as a
        parameter (e.g., Domcut), only 3 methods need to be explicitly
        written.
      </para>

      <para>
        1. _run() which mimics the web form, submits the sequence
        and extracts the raw text data of the result from the HTML.
      </para>

      <para>
        2. result() which parses the raw data into whatever useful format you
        wish. Usually these include SeqFeature objects, standard Perl data
        structures and meta sequences if there is a result for each residue in
        the sequence.
      </para>
      <para>
        3. _init() which imports the analysis specifications into the
        Analysis object
      </para>
      <para>
        As well as these methods, various file scoped lexical hashes need to
        be defined which hold constant data about the analysis, and the
        analysis result. These are useful for reference and are incorporated
        into the analysis objects.
      </para>

      <para>
        
        For more complicated prediction programs with analysis specific
        parameters (e.g., Sopma), get/set methods need to be written. Also, if
        any of these parameters need special error checking then a
        check_parameters() needs to be written as well. The nature of these
        parameters should be listed in a hash referenced by the $INPUT_SPEC
        variable.

      </para>

      <para>
        Alternatively, mail me with your suggestion and I'll try to put one
        together.  It is my intent to keep the modules up to date with new 
        analysis programs as they become available which are not included in
        the usual protein annotations and will be glad to hear of new 
        suggestions.
      </para>

   </section>

  <section id="Acknowledgments">
    <title>Acknowledgments</title>
    <para>
      These modules depend on several recently developed modules and wouldn't
      work at all without them. 
      Bio::Seq::Meta:: modules:
      Chad Matsalla, Aaron Mackey and Heikki Lehvaslaiho
      Bio::WebAgent:
      Heikki Lehvaslaiho
      Bio::SimpleAnalysisI:
      Martin Senger  
    </para>
  </section>

</article>
