<!DOCTYPE article PUBLIC "-//OASIS//DTD DocBook XML V4.2//EN">
<article lang="en">
  <articleinfo>
    <title>Bio::SearchIO HOWTO</title>
    <author>
      <!-- one of (CONTRIB AUTHORBLURB AFFILIATION OTHERNAME LINEAGE
      SURNAME FIRSTNAME HONORIFIC) -->
      <firstname>Jason</firstname>
      <surname>Stajich</surname>
      <authorblurb>
	<para>Bioperl Core Developer</para>
      </authorblurb>
      <affiliation>
	<orgname><ulink url="http://www.duke.edu">Duke
	University</ulink></orgname>
	<orgdiv><ulink url="http://upg.duke.edu">University Program in Genetics</ulink></orgdiv>
	<orgdiv><ulink url="http://cgt.genetics.duke.edu">Center for
	Genome Technology</ulink></orgdiv>
	<address>
	  Duke University Medical Center
	  <pob>Box 3568</pob>
	  <city>Durham, </city>
	  <state>North Carolina</state> 
	  <postcode>27710-3568</postcode> 
	  <country>USA</country>
	  <email>jason-at-bioperl.org</email>
        </address>
      </affiliation>
    </author>
    <author>
      <!-- one of (CONTRIB AUTHORBLURB AFFILIATION OTHERNAME LINEAGE
      SURNAME FIRSTNAME HONORIFIC) -->
      <firstname>Brian</firstname>
      <surname>Osborne</surname>
      <authorblurb>
	<para>Bioperl contributor</para>
      </authorblurb>
      <affiliation>
	<orgname><ulink url="http://www.cognia.com">Cognia Corporation</ulink></orgname>
	<address>
	  <city>NYC</city>, <state>NY</state> <postcode>10022</postcode> 
	  <country>USA</country>
	  <email>brian-at-cognia.com</email>
        </address>
      </affiliation>
    </author>
 <pubdate>2002-07-14</pubdate>
 
    <revhistory>
      <revision>
	<revnumber>0.1</revnumber>
	<date>2002-07-14</date>
	<authorinitials>js</authorinitials>
	<revremark>first draft</revremark>
      </revision>
      <revision>
	<revnumber>0.2</revnumber>
	<date>2002-10-11</date>
	<authorinitials>js</authorinitials>
	<revremark>added info on extending Search objects</revremark>
      </revision>
      <revision>
	<revnumber>0.3</revnumber>
	<date>2003-02-13</date>
	<authorinitials>BIO</authorinitials>
	<revremark>added table and text to Parsing section</revremark>
      </revision>
      <revision>
	<revnumber>0.4</revnumber>
	<date>2003-09-10</date>
	<authorinitials>BIO</authorinitials>
	<revremark>updated Parsing section</revremark>
      </revision>
    </revhistory>
    <legalnotice>
      <para>This document is copyright Jason Stajich, 2002.  For
	reproduction other than personal use please contact jason-at-bioperl.org
      </para>      
    </legalnotice>
    
    <abstract>
      <para>
	This is a HOWTO written in DocBook (SGML) for the reasoning behind
	the creation of the Bio::SearchIO system, how to use it, and how one
	goes about writing new adaptors to different output formats.  We will
	also describe how the Bio::SearchIO::Writer modules work for outputting
	various formats from Bio::Search objects. 
      </para>
    </abstract>
  </articleinfo>
  
  <section id="background">
    <title>Background</title>    
    <para>
      One of the most common and necessary tasks in bioinformatics is
      parsing analysis reports so that one can write programs which can
      help interpret the sheer volume of data that can be produced by
      processing many sequences.  To this end the Bioperl project has
      produced a number of parsers for the ubiquitous BLAST report.
      Steve Chervitz wrote one of the first Bioperl modules for BLAST
      called Bio::Tools::Blast.  Ian Korf allowed us to import and
      modify his BPlite (<emphasis>B</emphasis>last
      <emphasis>P</emphasis>arser) Bio::Tools::BPlite module into
      Bioperl.  This is of course in a sea of BLAST parsers that have
      been written by numerous people, but we will only cover the ones
      associated directly with the Bioperl project in this document.
      One of the reasons for writing yet another BLAST parser in the
      form of Bio::SearchIO is that even though both Bio::Tools::Blast
      and Bio::Tools::BPlite did their job correctly, and could parse
      WU-BLAST and NCBI-BLAST output, they did not adequately
      genericize what they were doing.  By this we mean everything was
      written around the BLAST format and was not easily applicable to
      parsing say, FastA alignments or a new alignment format.  One of
      the powerful features of the Object-Oriented framework in Bioperl
      is the ability to read in say, a sequence file, in different
      formats or from different data sources like a database or
      XML-flatfile, and have the program code process the sequences
      objects in the same manner.  We wanted to have this capability in
      place for analysis reports as well and thus the generic design of
      the Bio::SearchIO module.
    </para>
  </section>
  <section id="design">
    <title>Design</title>
    <para>
      The Bio::SearchIO system was designed with the following
      assumptions:  That all reports parsed with it could be separated
      into a hierarchy of components. The Result is the entire
      analysis for a single query sequence, and multiple Results can be
      concatenated together into a single file (i.e. running blastall
      with a fasta database as the input file rather than a single
      sequence).  Each Result is a set of Hits for the query sequence.
      Hits are sequences in the searched database which could be
      aligned to the query sequence and met the minimal search
      parameters, such as e-value threshold.  Each Hit has one or more
      High-scoring segment Pairs (HSPs) which are the alignments of the
      query and hit sequence.  Each Result has a set of one or more
      Hits and each Hit has a set of one or more HSPs, and this
      relationship can be used to describe results from all pairwise
      alignment programs including BLAST, FastA,
      and implementations of the Smith-Waterman and Needleman-Wunsch
      algorithms.
    </para>
    <para>
      A design pattern, called Factory, is utilized in object
      oriented programming to separate the entity which process data
      from objects which will hold the information produced.  In the
      same manner that the Bio::SeqIO module is used to parse different
      file formats and produces objects which are Bio::PrimarySeqI
      compliant, we have written Bio::SearchIO to produce the
      Bio::Search objects.  Sequences are a little less complicated so
      there is only one primary object (Bio::PrimarySeqI) which Search
      results need three main components to represent the data
      processed in a file: Bio::Search::Result::ResultI (top level results),
      Bio::Search::Hit::HitI (hits) and Bio::Search::HSP::HSPI (HSPs).   
      The Bio::SearchIO object is then a factory which produces
      Bio::Search::Result::ResultI objects  and the
      Bio::Search::Result::ResultI objects contain information about
      the query, the database searched, and the full collection of Hits
      found for the query.
    </para>
  </section>
  <section id="functionality">
    <title>New Functionality</title>    
    <para>
      The generality of the SearchIO approach is demonstrated by large number
      of report formats that have appeared since its introduction. These
      formats include AXT format reports (BLAT, BLASTZ), NCBI tabular
      output (-m 8 or -m 9 options), NCBI Blast XML, chadosxpr format flat databases,
      Exonerate output, FASTA output, hmmsearch output (HMMER),
      megablast output, PSL format output (BLAT), sim4 output, WABA
      output, and output from Wise.
    </para>
  </section>
  <section id="use">
    <title>Parsing with Bio::SearchIO</title>
    <para>
      This section is going to describe how to use the SearchIO system
      to process reports.  We'll describe BLAST reports but the
      idea is that once you understand the methods associated with the
      objects you won't need to know anything special about other 
      SearchIO parsers.
    </para>
    <para>
      Before we get into the details we should admit that there is
      some confusion about the names and functions of the objects for
      historical reasons. Both Steve Chervitz and Jason Stajich have implemented
      parsers in this system.  Steve created the psiblast parser (which
      does parse regular BLAST files too) and a host of objects named
      Bio::Search::XXX::BlastXXX where XXX is HSP, Hit, and Result.
      These objects are created by his Bio::SearchIO::psiblast
      implementation.  The objects Jason has created are called
      Bio::Search::XXX::GenericXXX where, again, XXX is HSP, Hit, and
      Result.  Because of some of the assumptions made in Steve's
      implementation and his utilization of what is known as 'lazy
      parsing', it is probably not going to be very easy to maintain his
      system without his help.  On the other hand Jason has tried 
      to make his implementations much easier to follow
      because all the parsing is done in one module.  
      </para>
    <para>
      The important take home message is that you cannot assume that methods in the
      BlastXXX objects are in fact implemented by the GenericHSP objects.
      More likely than not the BlastXXX objects will be deprecated and
      dismantled as their functionality is ported to the GenericHSP
      objects. For this reason we'll only be discussing the Generic*
      objects, though we'll use the terms 'hit', 'HSP', and 'result'. 
    </para>  
    <para>
      Here's example code which processes a BLAST report finding all the
      hits where the HSPs are greater than 100 residues and the percent
      identity is less than 75 percent. This code demonstrates
      that a result, in this case from a BLAST report, contains one or more hits,
      and a hit contains one or HSPs.
      <programlisting>
    use strict;
    use Bio::SearchIO;
   
    my $in = new Bio::SearchIO(-format => 'blast', 
                               -file   => 'report.bls');
    while( my $result = $in->next_result ) {
      while( my $hit = $result->next_hit ) {
        while( my $hsp = $hit->next_hsp ) {
          if( $hsp->length('total') > 100 ) {
            if ( $hsp->percent_identity >= 75 ) {
              print "Hit= ",       $hit->name, 
                   ",Length=",     $hsp->length('total'), 
                   ",Percent_id=", $hsp->percent_identity, "\n";
            }
          }
        }  
      }
    }
      </programlisting>
      The example above shows just a few of the many methods
      available in SearchIO. In order to display all these methods
      and what they return let's use a report as input, a simple BLASTX result: 
    </para>
    <programlisting>      
 BLASTX 2.2.4 [Aug-26-2002]


 Reference: Altschul, Stephen F., Thomas L. Madden, Alejandro A. Schaffer, 
 Jinghui Zhang, Zheng Zhang, Webb Miller, and David J. Lipman (1997), 
 "Gapped BLAST and PSI-BLAST: a new generation of protein database search
 programs",  Nucleic Acids Res. 25:3389-3402.

 Query= gi|20521485|dbj|AP004641.2 Oryza sativa (japonica
 cultivar-group) genomic DNA, chromosome 1, BAC clone:B1147B04, 3785
 bases, 977CE9AF checksum.
          (3059 letters)

 Database: test.fa 
            5 sequences; 1291 total letters



                                                                 Score    E
 Sequences producing significant alignments:                      (bits) Value

 gb|443893|124775 LaForas sequence                                 92   2e-022

 >gb|443893|124775 LaForas sequence
           Length = 331

  Score = 92.0 bits (227), Expect = 2e-022
  Identities = 46/52 (88%), Positives = 48/52 (91%)
  Frame = +1

 Query: 2896 DMGRCSSGCNRYPEPMTPDTMIKLYREKEGLGAYIWMPTPDMSTEGRVQMLP 3051
             D+ + SSGCNRYPEPMTPDTMIKLYRE EGL AYIWMPTPDMSTEGRVQMLP
 Sbjct: 197  DIVQNSSGCNRYPEPMTPDTMIKLYRE-EGL-AYIWMPTPDMSTEGRVQMLP 246


   Database: test.fa
     Posted date:  Feb 12, 2003  9:51 AM
   Number of letters in database: 1291
   Number of sequences in database:  5
  
 Lambda     K      H
    0.318    0.135    0.401 

 Gapped
 Lambda     K      H
    0.267   0.0410    0.140 


 Matrix: BLOSUM62
 Gap Penalties: Existence: 11, Extension: 1
 Number of Hits to DB: 7140
 Number of Sequences: 5
 Number of extensions: 180
 Number of successful extensions: 2
 Number of sequences better than 10.0: 2
 Number of HSP's better than 10.0 without gapping: 1
 Number of HSP's successfully gapped in prelim test: 0
 Number of HSP's that attempted gapping in prelim test: 0
 Number of HSP's gapped (non-prelim): 1
 length of database: 1291
 effective HSP length: 46
 effective length of database: 1061
 effective search space used:  1032353
 frameshift window, decay const: 50,  0.1
 T: 12
 A: 40
 X1: 16 ( 7.3 bits)
 X2: 38 (14.6 bits)
 X3: 64 (24.7 bits)
 S1: 32 (17.6 bits)
    </programlisting>
    <para>
      Table 1 shows all the data returned by methods used by the
      Result, Hit, and HSP objects when the report shown above is used as
      input. Note that many of the methods shown can be used to either
      get or set values, but we're just showing what they get.
    </para>
    <table>
      <title>SearchIO Methods</title>
      <tgroup cols="4">
	<thead>
	  <row>
	    <entry>Object</entry>
	    <entry>Method</entry>
	    <entry>Example</entry>
	    <entry>Description</entry>
	  </row>
	</thead>
	<tbody>
	  <row>
	    <entry>Result</entry>
	    <entry>algorithm</entry>
	    <entry>BLASTX</entry>
	    <entry>algorithm</entry>
	  </row>
	  <row>
	    <entry>Result</entry>
	    <entry>algorithm_version</entry>
	    <entry>2.2.4 [Aug-26-2002]</entry>
	    <entry>algorithm version</entry>
	  </row>
	  <row>
	    <entry>Result</entry>
	    <entry>query_name</entry>
	    <entry>gi|20521485|dbj|AP004641.2</entry>
	    <entry>query name</entry>
	  </row>
	  <row>
	    <entry>Result</entry>
	    <entry>query_accession</entry>
	    <entry>AP004641.2</entry>
	    <entry>query accession</entry>
	  </row>
	  <row>
	    <entry>Result</entry>
	    <entry>query_length</entry>
	    <entry>3059</entry>
	    <entry>query length</entry>
	  </row>
	  <row>
	    <entry>Result</entry>
	    <entry>query_description</entry>
	    <entry>Oryza sativa ... 977CE9AF checksum.</entry>
	    <entry>query description</entry>
	  </row>
	  <row>
	    <entry>Result</entry>
	    <entry>database_name</entry>
	    <entry>test.fa </entry>
	    <entry>database name</entry>
	  </row>
	  <row>
	    <entry>Result</entry>
	    <entry>database_letters</entry>
	    <entry>1291</entry>
	    <entry>number of residues in database</entry>
	  </row>
	  <row>
	    <entry>Result</entry>
	    <entry>database_entries</entry>
	    <entry>5</entry>
	    <entry>number of database entries</entry>
	  </row>
	  <row>
	    <entry>Result</entry>
	    <entry>available_statistics</entry>
	    <entry>effectivespaceused ... dbletters</entry>
	    <entry>statistics used</entry>
	  </row>
	  <row>
	    <entry>Result</entry>
	    <entry>available_parameters</entry>
	    <entry>gapext matrix allowgaps gapopen</entry>
	    <entry>parameters used</entry>
	  </row>
	  <row>
	    <entry>Result</entry>
	    <entry>num_hits</entry>
	    <entry>1</entry>
	    <entry>number of hits</entry>
	  </row>
	  <row>
	    <entry>Result</entry>
	    <entry>hits</entry>
	    <entry> </entry>
	    <entry>Search::Hit::GenericHit object</entry>
	  </row>
	  <row>
	    <entry>Hit</entry>
	    <entry>name</entry>
	    <entry>gb|443893|124775</entry>
	    <entry>hit name</entry>
	  </row>
	  <row>
	    <entry>Hit</entry>
	    <entry>accession</entry>
	    <entry>443893</entry>
	    <entry>accession</entry>
	  </row>
	  <row>
	    <entry>Hit</entry>
	    <entry>description</entry>
	    <entry>LaForas sequence</entry>
	    <entry>hit description</entry>
	  </row>
	  <row>
	    <entry>Hit</entry>
	    <entry>algorithm</entry>
	    <entry>BLASTX</entry>
	    <entry>algorithm</entry>
	  </row>
	  <row>
	    <entry>Hit</entry>
	    <entry>raw_score</entry>
	    <entry>92</entry>
	    <entry>hit raw score</entry>
	  </row>
	  <row>
	    <entry>Hit</entry>
	    <entry>significance</entry>
	    <entry>2e-022</entry>
	    <entry>hit significance</entry>
	  </row>
	  <row>
	    <entry>Hit</entry>
	    <entry>bits</entry>
	    <entry>92.0</entry>
	    <entry>hit bits</entry>
	  </row>
	  <row>
	    <entry>Hit</entry>
	    <entry>hsps</entry>
	    <entry> </entry>
	    <entry>Search::HSP::GenericHSP object</entry>
	  </row>
	  <row>
	    <entry>Hit</entry>
	    <entry>num_hsps</entry>
	    <entry>1</entry>
	    <entry>number of HSPs in hit</entry>
	  </row>
	  <row>
	    <entry>Hit</entry>
	    <entry>locus</entry>
	    <entry>124775</entry>
	    <entry>locus name</entry>
	  </row>
	  <row>
	    <entry>Hit</entry>
	    <entry>accession_number</entry>
	    <entry>443893</entry>
	    <entry>accession number</entry>
	  </row>
	  <row>
	    <entry>HSP</entry>
	    <entry>algorithm</entry>
	    <entry>BLASTX</entry>
	    <entry>algorithm</entry>
	  </row>
	  <row>
	    <entry>HSP</entry>
	    <entry>evalue</entry>
	    <entry>2e-022</entry>
	    <entry>e-value</entry>
	  </row>
	   <row>
	    <entry>HSP</entry>
	    <entry>expect</entry>
	    <entry>2e-022</entry>
	    <entry>alias for evalue()</entry>
	  </row>
	  <row>
	    <entry>HSP</entry>
	    <entry>frac_identical</entry>
	    <entry>0.884615384615385</entry>
	    <entry>Fraction identical</entry>
	  </row>
	  <row>
	    <entry>HSP</entry>
	    <entry>frac_conserved</entry>
	    <entry>0.923076923076923</entry>
	    <entry>desc</entry>
	  </row>
	  <row>
	    <entry>HSP</entry>
	    <entry>gaps</entry>
	    <entry>2</entry>
	    <entry>number of gaps</entry>
	  </row>
	  <row>
	    <entry>HSP</entry>
	    <entry>query_string</entry>
	    <entry>DMGRCSSG ...</entry>
	    <entry>string from alignment</entry>
	  </row>
	  <row>
	    <entry>HSP</entry>
	    <entry>hit_string</entry>
	    <entry>DIVQNSS ...</entry>
	    <entry>string from alignment</entry>
	  </row>
	  <row>
	    <entry>HSP</entry>
	    <entry>homology_string</entry>
	    <entry>D+ + SSGCN ...</entry>
	    <entry>string from alignment</entry>
	  </row>
	  <row>
	    <entry>HSP</entry>
	    <entry>length('total')</entry>
	    <entry>52</entry>
	    <entry>length of HSP</entry>
	  </row>
	  <row>
	    <entry>HSP</entry>
	    <entry>length('hit')</entry>
	    <entry>50</entry>
	    <entry>length of hit minus gaps</entry>
	  </row>
	  <row>
	    <entry>HSP</entry>
	    <entry>length('query')</entry>
	    <entry>156</entry>
	    <entry>length of query minus gaps</entry>
	  </row>
	  <row>
	    <entry>HSP</entry>
	    <entry>hsp_length</entry>
	    <entry>52</entry>
	    <entry>desc</entry>
	  </row>
	  <row>
	    <entry>HSP</entry>
	    <entry>frame</entry>
	    <entry>0</entry>
	    <entry>frame, GFF convention</entry>
	  </row>
	  <row>
	    <entry>HSP</entry>
	    <entry>num_conserved</entry>
	    <entry>48</entry>
	    <entry>number of conserved residues</entry>
	  </row>
	  <row>
	    <entry>HSP</entry>
	    <entry>num_identical</entry>
	    <entry>46</entry>
	    <entry>number of identical residues</entry>
	  </row>
	  <row>
	    <entry>HSP</entry>
	    <entry>rank</entry>
	    <entry>1</entry>
	    <entry>rank of HSP</entry>
	  </row>
	  <row>
	    <entry>HSP</entry>
	    <entry>seq_inds('query','identical')</entry>
	    <entry>(966,971,972,973,974,975 ...)</entry>
	    <entry>identical positions as array</entry>
	  </row>
	  <row>
	    <entry>HSP</entry>
	    <entry>seq_inds('query','conserved')</entry>
	    <entry>(967,969)</entry>
	    <entry>conserved positions as array</entry>
	  </row>
	  <row>
	    <entry>HSP</entry>
	    <entry>seq_inds('hit','identical')</entry>
	    <entry>(197,202,203,204,205 ...)</entry>
	    <entry>identical positions as array</entry>
	  </row>
	  <row>
	    <entry>HSP</entry>
	    <entry>seq_inds('hit','conserved')</entry>
	    <entry>(198,200)</entry>
	    <entry>conserved positions as array</entry>
	  </row>
	  <row>
	    <entry>HSP</entry>
	    <entry>score</entry>
	    <entry>227</entry>
	    <entry>score</entry>
	  </row>
	  <row>
	    <entry>HSP</entry>
	    <entry>bits</entry>
	    <entry>92.0</entry>
	    <entry>bits</entry>
	  </row>
	  <row>
	    <entry>HSP</entry>
	    <entry>range('query')</entry>
	    <entry>(2896,3051)</entry>
	    <entry>start and end as array</entry>
	  </row>
	  <row>
	    <entry>HSP</entry>
	    <entry>range('hit')</entry>
	    <entry>(197,246)</entry>
	    <entry>start and end as array</entry>
	  </row>
	  <row>
	    <entry>HSP</entry>
	    <entry>percent_identity</entry>
	    <entry>88.4615384615385</entry>
	    <entry>% identical</entry>
	  </row>
	  <row>
	    <entry>HSP</entry>
	    <entry>strand('hit')</entry>
	    <entry>1</entry>
	    <entry>strand of the hit</entry>
	  </row>
	  <row>
	    <entry>HSP</entry>
	    <entry>strand('query')</entry>
	    <entry>1</entry>
	    <entry>strand of the query</entry>
	  </row>
	  <row>
	    <entry>HSP</entry>
	    <entry>start('query')</entry>
	    <entry>2896</entry>
	    <entry>start position from alignment</entry>
	  </row>
	  <row>
	    <entry>HSP</entry>
	    <entry>end('query')</entry>
	    <entry>3051</entry>
	    <entry>end position from alignment</entry>
	  </row>
	  <row>
	    <entry>HSP</entry>
	    <entry>start('hit')</entry>
	    <entry>197</entry>
	    <entry>start position from alignment</entry>
	  </row>
	  <row>
	    <entry>HSP</entry>
	    <entry>end('hit')</entry>
	    <entry>246</entry>
	    <entry>end position from alignment</entry>
	  </row>
	  <row>
	    <entry>HSP</entry>
	    <entry>matches('hit')</entry>
	    <entry>(46,48)</entry>
	    <entry>number of identical and conserved as array</entry>
	  </row>
	  <row>
	    <entry>HSP</entry>
	    <entry>matches('query')</entry>
	    <entry>(46,48)</entry>
	    <entry>number of identical and conserved as array</entry>
	  </row>
	  <row>
	    <entry>HSP</entry>
	    <entry>alignment</entry>
	    <entry> </entry>
	    <entry>Bio::SimpleAlign object</entry>
	  </row>
	</tbody>
      </tgroup>
    </table>
    
    <para>
      Table 1 shows that a method can return a string, an array,
      or an object. When an object is returned some
      additional code will probably be needed to get the data of
      interest. For example, if you wanted a printable alignment after
      you'd parsed BLAST output you could use the get_aln() method,
      retrieve a Bio::SimpleAlign object and use it like this:
      <programlisting>
	use Bio::AlignIO;
	# $aln will be a Bio::SimpleAlign object
	my $aln = $hsp->get_aln;
	my $alnIO = Bio::AlignIO->new(-format=>"msf");
	my $alignment_as_string = $alnIO->write_aln($aln);
   </programlisting>
      On one hand it appears to be a complication, but by entering the
      worlds of the AlignIO and SimpleAlign objects you now have access to their
      functionality and flexibility. This is the beauty of Bioperl!
    </para>
    <para>
      Some of these methods deserve a bit more explanation since they do
      more than simply extract data directly from the output. For example, the
      ambiguous_aln() method is designed to tell us whether two or more
      HSPs from a given hit overlap, and whether the overlap refers to the
      queries or the hits, or both. One situation is where overlaps would be
      found in one but not the other arises where there are repeats in the
      query or hit. The ambiguous_aln() method will return one of
      these 4 values:
      <simplelist type="horiz" columns="2">
	<member>q</member>
	<member>query sequence contains overlapping sub-sequences
             while hit sequence does not</member>
	<member>s</member>
	<member>hit sequence contains overlapping sub-sequences
             while query does not</member>
	<member>qw</member>
	<member>query and hit sequences contain overlapping sub-sequences
             relative to each other</member>
	<member>-</member>
	<member>query and hit sequence do not contain multiple domains
             relative to each other OR both contain the same distribution
             of similar domains</member>
      </simplelist>
    </para>
    <para>
      Another method that's useful in dissecting an HSP is the
      seq_inds() method of the HSP object. What this method does is tell
      us what the positions are of all the identical or conserved
      residues in an alignment, query or hit. It could be used like this:
      <programlisting>
	# put all the conserved matches in query strand into an array
	my @str_array = split "",$hsp->query_strand;
	foreach ( $hsp->seq_inds('query','conserved') ){
	   push @conserved,$str_array[$_ - 1];
	}
    </programlisting>
    </para>
    <para>
      In most cases the SearchIO methods extract data directly from
      output but there's one important exception, the frame() method of the HSP
      object. Instead of using the values in the BLAST report it converts
      them to values according to the GFF specification, which is a format used by
      many Bioperl modules involved in gene annotation (for more on GFF see
      <ulink
	url="http://www.sanger.ac.uk/Software/formats/GFF/GFF_Spec.shtml">http://www.sanger.ac.uk/Software/formats/GFF/GFF_Spec.shtml</ulink>).
    </para>
    <para>
      Specifically, the frame() method returns 0, 1, or 2 instead of
      the expected -3, -2, -1, +1, +2, or +3 in BLAST. GFF frame values are
      meaningful relative to the strand of the
      hit or query sequence so in order to reconstruct the BLAST
      frame you need to both the strand, 1 or -1, and the GFF frame value:
      <programlisting>
	my $blast_frame = ($hsp->query->frame + 1) * $hsp->query->strand;
      </programlisting>
    </para>
    <para>
      Our simple table of methods does not show all available
      arguments or returned values for all the SearchIO methods. The
      best place to explore any method in detail is 
      <ulink url="http://doc.bioperl.org">http://doc.bioperl.org</ulink>
      which provides the HTML versions of the Perl POD (Plain Old
      Documentation) that is embedded in every well-written Perl
      module. Other sources of code are the examples/searchio/
      and scripts/searchio directories in the Bioperl package.
    </para>
  </section>
  
  <section id="creating_reports">
    <title>Creating Reports for SearchIO</title>
    <para>
      One note on creating reports that can be parsed by SearchIO: the
      developers haven't attempted to parse all the possible reports 
      that could be created by programs with many command-line options,
      like blastall. Certainly you should be able to parse reports
      created using the default settings, but if you're running
      blastall, say, using some special set of options and you've
      encountered a parsing problem this may be the explanation.
    </para>
    <para>
      For example, one can currently parse BLAST output created with the default
      settings as well as the reports created when using the "-m 8"
      or "-m 9" options (use format "blasttable") or the XML-formatted
      reports but it's still possible to find sets of options that
      SearchIO can't parse.
    </para>
  </section>

  <section id="implementation">
    <title>Implementation</title>
    <para>
      This section is going to describe how the SearchIO system was
      implemented, it is probably not necessary to understand all of
      this unless you are curious or want to implement your own
      Bio::SearchIO parser.  We have utilized an event-based system to
      process these reports.  This is analagous to the SAX 
      (<emphasis>S</emphasis>imple <emphasis>A</emphasis>PI
      for <emphasis>X</emphasis>ML) system used to process XML documents.  Event based
      parsing can be simply thought of as simple start and end events.
      When you hit the beginning of a report a start event is thrown,
      when you hit the end of the report an end event is thrown.  So
      the report events are paired, and everything else that is thrown
      in between the paired start and end events is related to that
      report.  Another way to think of it is as if you pick a number
      and color for a card in a standard deck.  Let's say you pick red
      and 2.  The you start dealing cards from our deck and pile them
      one on top of each other.  When you see your first red 2 you
      start a new pile, and start dealing cards onto that pile until
      you see the next red 2.  Everything in your pile that happened
      between when you saw the beginning red 2 and ending red 2 is data
      you'll want to keep and process.  In the same way all the events
      you see between a pair of start and end events (like 'report' or
      'hsp') are data associated with object or child object in its
      hierarchy.  A listener object processes all of these events, in
      our example the listener is the table where the stack of cards is
      sitting, and later it is the hand which moves the pile of cards
      when a new stack is started.  The listener will take the events
      and process them.  We've neglected to tell you of a third event
      that is thrown and caught.  This is the characters event in SAX
      terminology, which is simply data.  So one sends a start event,
      then some data, then an end event.  This process is analagous to
      a finite state machine in computer science (and I'm sure the
      computer scientists reading this right are already yawning) where
      what we do with data received is dependent on the state we're in.
      The state that the listener is in is affected by the events that
      are processed.
    </para>
    <para>
      A small caveat: in an ideal situation a processor would throw
      events and not need to maintain any state information, it would
      just be processing data and the listener would manage the
      information and state.  However, a lot of the parsing of these
      human-readable reports requires contextual information to apply
      the correct regular expressions.  So in fact the event thrower
      has to know what state it is in and apply different methods based
      on this.  In contrast the XML parsers simply keep track of what
      state they are in, but can process all the data with the same
      system of reading the tag and sending the data that is in between
      the XML start and end tags.
    </para>
    <para>
      All of this framework has been built up so to implement a new
      parser one only needs to write a module that produces the
      appropriate start and end events and the existing framework will
      do the work of creating the objects for you.  Here's how we've
      implemented event-based parsing for Bio::SearchIO.  The
      Bio::SearchIO is just the front-end to this process, in fact the
      processing of these reports is done by different modules in the
      Bio/SearchIO directory.  So if you look at your bioperl
      distribution and the modules in Bio/SearchIO you'll see modules in
      there like blast.pm, fasta.pm, blastxml.pm,
      SearchResultEventBuilder.pm, EventHandlerI.pm (depending on what
      version of the toolkit there may be more modules in there).
      There is also a SearchWriterI.pm and Writer directory in there
      but we'll save that for later.  If you don't have the
      distribution handy you can navigate this at the bioperl <ulink
	url="http://cvs.open-bio.org/cgi-bin/viewcvs/viewcvs.cgi/bioperl-live/Bio/SearchIO/?cvsroot=bioperl">CVS web page</ulink>.  
    </para>
    <para>
      Let's use the blast.pm module as an example to describe the
      relationship of the modules in this dir (could have substituted
      any of the other format parsers like fasta.pm or blastxml.pm -
      these are always lowercase for historical reasons).  The module
      has some features you should look for - the first is the hash in
      the BEGIN block called %MAPPING.  This key/value pairs here are
      the shorthand for how we map events from this module to general
      event names.  This is only necessary because if we have an XML
      processor (see the blastxml.pm module) the event names will be
      the same as the XML tag names (like &lt;Hsp_bit-score&gt; in the
      NCBI BLAST XML DTD).  So to make this general we'll make sure all
      of the events inside our parser map to the values in the %MAPPING
      hash - we can call them whatever we want inside this module.
      Some of the events map to hash references (like
      Statistics_db-len) and this is so we can map multiple values to the
      same top-level attribute field but we know they will be stored as
      a hash value in the subsequent object (in this example, keyed by
      the name 'dbentries').  The capital "RESULT", "HSP", or "HIT" in
      the value name allow us to encode the event state in the event so
      we don't have to pass in two values.  It also easy for someone to
      quickly read the list of events and know which ones are related
      to Hits and which ones are related to HSPs.  The listener in our
      architecture is the Bio::SearchIO::SearchResultEventBuilder.
      This object is attached as a listener through the Bio::SearchIO
      method add_EventListener.  In fact you could have multiple event
      listeners and they could do different things.  In our case we
      want to create Bio::Search objects, but an event listener could
      just as easily be propagating data directly into a database based
      on the events.  The SearchResultEventBuilder takes the events
      thrown by the SearchIO classes and builds the appropriate
      Bio::Search::HSP object from it.
    </para>
    <para>
      Sometimes special objects are needed that are extensions beyond
      what the GenericHSP or GenericHit objects are meant to
      represent.  For this case we have implemented
      Bio::SearchIO::SearchResultEventBuilder so that it can use
      factories for creating its resulting Bio::Search objects - see
      the Bio::SearchIO::hmmer::_initialize method for an example of how
      this can be set.
    </para>
  </section>
  <section id="outputting">
    <title>Writing and formatting output</title>
    <para>
      Often people want to write back out a BLAST report for users who are
      most comfortable with that output or if you want to visualize the
      context of a weakly aligned region to use human intuition to score
      the confidence of a putative homologue.  Bio::SearchIO is for parsing
      in the data but Bio::SearchIO::Writer is for outputting the
      information. The simplest way to output data as a pseudo-BLAST
      HTML format is as follows.
    </para>
    <para>
      <programlisting>
   my $writerhtml = new Bio::SearchIO::Writer::HTMLResultWriter();
   my $outhtml = new Bio::SearchIO(-writer => $writerhtml,
				   -file   => ">searchio.html");
   # get a result from Bio::SearchIO parsing or build it up in memory
   $outhtml->write_result($result);
      </programlisting>
    </para>
    <para>
      If you wanted to get the output as a string rather than write it out to
      a file, simply use the following.
      <programlisting>
   $writerhtml->to_string($result);
      </programlisting>
    </para>
    <para>
      The HTMLResultWriter supports setting your own remote database url for
      the sequence links in the event you'd like to point to your own 
      SRS or local HTTP-based connection to the sequence data. Simply use
      the remote_database_url method which accepts a sequence type as input
      (protein or nucleotide).
    </para>
    <para>
      You can also override the id_parser() method to
      define what the unique IDs are from these sequence ids in the event
      you would like to use something other than the accession number that is
      gleaned from the sequence string.
    </para>
    <para>
      If your data is instead stored in a database you could build the
      Bio::Search objects up in memory directly from your database and
      then use the Writer object to output the data.  Currently there is
      also a Bio::SearchIO::Writer::TextResultWriter which supports
      writing BLAST text file output.
    </para>
  </section>
  <section id="extending">
    <title>Extending SearchIO</title>
    <para>
      The framework for Bio::SearchIO is just a starting point for parsing
      these reports and creating objects which represent the information.
      If you would like to create your own set of objects which extend the
      current functionality we have built the system so that it will
      support this.  For example, you may have built your own HSP object
      which supports a special operation like realign_with_sw(), which might
      realign the HSP via a Smith-Waterman algorithm, pulling extra bases
      from the flanking sequence.  You might call your module
      Bio::Search::HSP::RealignHSP and put it in a file called
      Bio/Search/HSP/RealignHSP.pm.  Note that you don't have to put this
      file directly in the bioperl source directory - you can create your own
      local directory structure that is in parallel to the bioperl release
      source code as long as you have updated your PERL5LIB to contain your
      local directory or you are using the 'use lib' directive in your script.
      Also, you don't have to use the namespace Bio::Search::HSP as
      namespaces don't mean anything to perl with respect to object inheritance, but
      do we recommend you name things in a logical manner so that others
      might read and understand your code (and if you feel encouraged
      to donate your code to the project it might easily integrated
      with existing modules).
    </para>
    <para>
      So, you're going to write your new special module, you do need to
      make sure it inherits from the base Bio::Search::HSP::HSPI object.
      Additionally unless you want to reimplement all the initialization
      state in the current Bio::Search::HSP::GenericHSP you should just
      plan to extend that object.  You need to follow the chained
      constructor system that we have set up so that the arguments are
      properly processed.  Here is a sample of what your code might look
      like (don't forget to write your own POD so that it will be
      documented, we've left it off here to keep things simple).
    </para>
    <para>
      <programlisting>
   package Bio::Search::HSP::RealignHSP;
   use strict;
   use Bio::Search::HSP::GenericHSP;
   use vars qw(@ISA); # for inheritance
   @ISA = qw(Bio::Search::HSP::GenericHSP); # RealignHSP inherits from GenericHSP

   sub new { 
      my ($class,@args) = @_;
      my $self = $class->SUPER::new(@args); # chained contructor

      # process the 1 additional argument this object supports
      my ($ownarg1) = $self->_rearrange([OWNARG1],@args); 

      return $self; # remember to pass the object reference back out   
   }

   sub realign_hsp { 
      my ($self) = @_;
      # implement my special realign method here
   }
      </programlisting>
    </para>
    <para>
      The above code gives you a skeleton of how to start to implement
      your object. To register it so that it is used when the SearchIO
      system makes HSPs you just need to call a couple of functions.  The
      code below outlines them.   
    </para>
    <para>
      <programlisting>
use Bio::SearchIO;
use Bio::Search::HSP::HSPFactory;
use Bio::Search::Hit::HitFactory;

# setup the blast parser, you can do this with and SearchIO parser however
my $searchio = new Bio::SearchIO(-file => $blastfile, 
                                 -format =>'blast'); 
# build HSP factory with a certain type of HSPs to make
# the default is Bio::Search::HSP::GenericHSP
my $hspfact = new Bio::Search::HSP::HSPFactory(-type => 
                 'Bio::Search::HSP::RealignHSP');
# if you wanted to replace the Hit factory you can do this as well
# additionally there is an analagous
# Bio::Search::Result::ResultFactory for setting custom Result objects
my $hitfact = new Bio::Search::Hit::HitFactory(-type =>
                  'Bio::Search::Hit::SUPERDUPER_Hit');
$searchio->_eventHandler->register_factory('hsp', $hspfact);
$searchio->_eventHandler->register_factory('hit', $hitfact);
      </programlisting>
    </para>
    <para>
      We have to register the HSPFactory, which is the object which will
      create HSPI objects, by allowing this to be built by a factory
      rather than a hard-coded Bio::Search::HSP::GenericHSP->new(...) call.
      We are allowing the user to take advantage of the whole parsing
      structure and the ability to slot their own object into the process
      rather than re-implementing very much.  We think this is very
      powerful and worth the system overhead, but it may not permit this
      to be as efficient in parsing as we would like.  Future work will
      hopefully address speed and memory issues with this parser.
      Volunteers and improvement code are always welcome.
    </para>
  </section>
</article>
