comparison variant_effect_predictor/Bio/Tools/Blast.pm @ 0:1f6dce3d34e0

Uploaded

0:1f6dce3d34e0

#----------------------------------------------------------------------------
# PACKAGE : Bio::Tools::Blast
# PURPOSE : To encapsulate code for running, parsing, and analyzing
#           BLAST reports.
# AUTHOR  : Steve Chervitz ([e-mail])
# CREATED : March 1996
# REVISION: $Id: Blast.pm,v 1.30 2002/11/04 09:12:50 heikki Exp $
# STATUS  : Alpha
#
# For the latest version and documentation, visit:
#    http://bio.perl.org/Projects/Blast
#
# To generate documentation, run this module through pod2html
# (preferably from Perl v5.004 or better).
#
# Copyright (c) 1996-2000 Steve Chervitz. All Rights Reserved.
#           This module is free software; you can redistribute it and/or
#           modify it under the same terms as Perl itself.
#----------------------------------------------------------------------------

package Bio::Tools::Blast;
use strict;
use Exporter;

use Bio::Tools::SeqAnal;
use Bio::Root::Global     qw(:std);
use Bio::Root::Utilities  qw(:obj);

require 5.002;

use vars qw(@ISA @EXPORT @EXPORT_OK %EXPORT_TAGS
            $ID $VERSION $Blast @Blast_programs $Revision $Newline);

@ISA        = qw( Bio::Tools::SeqAnal Exporter);
@EXPORT     = qw();
@EXPORT_OK  = qw($VERSION $Blast);
%EXPORT_TAGS = ( obj => [qw($Blast)],
		 std => [qw($Blast)]);

$ID = 'Bio::Tools::Blast';
$VERSION  = 0.09; 
$Revision = '$Id: Blast.pm,v 1.30 2002/11/04 09:12:50 heikki Exp $';  #'

## Static Blast object.
$Blast = {};
bless $Blast, $ID;
$Blast->{'_name'} = "Static Blast object";

@Blast_programs  = qw(blastp blastn blastx tblastn tblastx);

use vars qw($DEFAULT_MATRIX $DEFAULT_SIGNIF);
my $DEFAULT_MATRIX   = 'BLOSUM62';
my $DEFAULT_SIGNIF   = 999;# Value used as significance cutoff if none supplied.
my $MAX_HSP_OVERLAP  = 2;  # Used when tiling multiple HSPs.

## POD Documentation:

=head1 NAME

Bio::Tools::Blast - Bioperl BLAST sequence analysis object

=head1 SYNOPSIS

=head2 Parsing Blast reports

Parse an existing Blast report from file:

    use Bio::Tools::Blast;

    $blastObj = Bio::Tools::Blast->new( -file   => '/tmp/blast.out',
					-parse  => 1,
					-signif => '1e-10',
					);

Parse an existing Blast report from STDIN:

    $blastObj = Bio::Tools::Blast->new( -parse  => 1,
					-signif => '1e-10',
					);

Then send a Blast report to your script via STDIN.

Full parameters for parsing Blast reports.

 %blastParam = (
		-run             => \%runParam,
		-file            => '',
		-parse           => 1,
		-signif          => 1e-5,
		-filt_func       => \&my_filter,
		-min_len         => 15,
		-check_all_hits  => 0,
		-strict          => 0,
		-stats           => 1,
		-best            => 0,
		-share           => 0,
		-exec_func       => \&process_blast,
		-save_array      => \@blast_objs,  # not used if -exce_func defined.
	       );

See L<parse()|parse> for a description of parameters and see L<USAGE | USAGE> for
more examples including how to parse streams containing multiple Blast
reports L<Using the Static $Blast Object>.

See L<Memory Usage Issues> for information about how to make Blast
parsing be more memory efficient.

=head2 Running Blast reports

Run a new Blast2 at NCBI and then parse it:

    %runParam = (
		  -method   => 'remote',
		  -prog     => 'blastp',
		  -database => 'swissprot',
		  -seqs     => [ $seq ],  # Bio::Seq.pm objects.
		  );

    $blastObj = Bio::Tools::Blast->new( -run     => \%runParam,
					-parse   => 1,
					-signif  => '1e-10',
					-strict  => 1,
					);

Full parameters for running Blasts at NCBI using Webblast.pm:

 %runParam = (
	      -method   => 'remote',
	      -prog     => 'blastp',
	      -version  => 2,      # BLAST2
	      -database =>'swissprot',
	      -html     => 0,
	      -seqs     => [ $seqObject ],  # Bio::Seq.pm object(s)
	      -descr    => 250,
	      -align    => 250,
	      -expect   => 10,
	      -gap      => 'on',
	      -matrix   => 'PAM250',
	      -email    => undef,  # don't send report via e-mail if parsing.
	      -filter   => undef,  # use default
	      -gap_c    => undef,  # use default
	      -gap_e    => undef,  # use default
	      -word     => undef,  # use default
	      -min_len  => undef,  # use default
	      );

See L<run()|run> and L<USAGE | USAGE> for more information about running Blasts.

=head2 HTML-formatting Blast reports

Print an HTML-formatted version of a Blast report:

    use Bio::Tools::Blast qw(:obj);

    $Blast->to_html($filename);
    $Blast->to_html(-file   => $filename,
		    -header => "<H1>Blast Results</H1>");
    $Blast->to_html(-file   => $filename,
		    -out    => \@array);  # store output
    $Blast->to_html();  # use STDIN

Results are sent directly to STDOUT unless an C<-out =E<gt> array_ref>
parameter is supplied. See L<to_html()|to_html> for details.

=head1 INSTALLATION

This module is included with the central Bioperl distribution:

   http://bio.perl.org/Core/Latest
   ftp://bio.perl.org/pub/DIST

Follow the installation instructions included in the README file.

=head1 DESCRIPTION

The Bio::Tools::Blast.pm module encapsulates data and methods for
running, parsing, and analyzing pre-existing BLAST reports. This
module defines an application programming interface (API) for working
with Blast reports. A Blast object is constructed from raw Blast
output and encapsulates the Blast results which can then be accessed
via the interface defined by the Blast object.

The ways in which researchers use Blast data are many and varied. This
module attempts to be general and flexible enough to accommodate
different uses. The Blast module API is still at an early stage of
evolution and I expect it to continue to evolve as new uses for Blast
data are developed. Your L<FEEDBACK | FEEDBACK> is welcome.

B<FEATURES:>

=over 2

=item * Supports NCBI Blast1.x, Blast2.x, and WashU-Blast2.x, gapped
and ungapped.

Can parse HTML-formatted as well as non-HTML-formatted reports.

=item * Launch new Blast analyses remotely or locally.

Blast objects can be constructed directly from the results of the
run. See L<run()|run>.

=item * Construct Blast objects from pre-existing files or from a new run.

Build a Blast object from a single file or build multiple Blast
objects from an input stream containing multiple reports. See
L<parse()|parse>.

=item * Add hypertext links from a BLAST report.

See L<to_html()|to_html>.

=item * Generate sequence and sequence alignment objects from HSP
sequences.

If you have Bio::Seq.pm and Bio::UnivAln.pm installed on your system,
they can be used for working with high-scoring segment pair (HSP)
sequences in the Blast alignment.  (A new version of Bio::Seq.pm is
included in the distribution, see L<INSTALLATION | INSTALLATION>).  For more
information about them, see:

    http://bio.perl.org/Projects/Sequence/
    http://bio.perl.org/Projects/SeqAlign/

=back

A variety of different data can be extracted from the Blast report by
querying the Blast.pm object. Some basic examples are given in the
L<USAGE | USAGE> section. For some working scripts, see the links provided in
the L<the DEMO SCRIPTS section | DEMO> section.

As a part of the incipient Bioperl framework, the Bio::Tools::Blast.pm
module inherits from B<Bio::Tools::SeqAnal.pm>, which provides some
generic functionality for biological sequence analysis. See the
documentation for that module for details 
(L<Links to related modules>).

=head2 The BLAST Program

BLAST (Basic Local Alignment Search Tool) is a widely used algorithm
for performing rapid sequence similarity searches between a single DNA
or protein sequence and a large dataset of sequences.  BLAST analyses
are typically performed by dedicated remote servers, such as the ones
at the NCBI. Individual groups may also run the program on local
machines.

The Blast family includes 5 different programs:

              Query Seq        Database
             ------------     ----------
 blastp  --  protein          protein
 blastn  --  nucleotide       nucleotide
 blastx  --  nucleotide*      protein
 tblastn --  protein          nucleotide*
 tblastx --  nucleotide*      nucleotide*

            * = dynamically translated in all reading frames, both strands

See L<References & Information about the BLAST program>.

=head2 Versions Supported

BLAST reports generated by different application front ends are similar
but not exactly the same. Blast reports are not intended to be exchange formats,
making parsing software susceptible to obsolescence. This module aims to
support BLAST reports generated by different implementations:

  Implementation    Latest version tested
  --------------    --------------------
  NCBI Blast1       1.4.11   [24-Nov-97]
  NCBI Blast2       2.0.8    [Jan-5-1999]
  WashU-BLAST2      2.0a19MP [05-Feb-1998]
  GCG               1.4.8    [1-Feb-95]

Support for both gapped and ungapped versions is included. Currently, there
is only rudimentary support for PSI-BLAST in that these reports can be parsed but
there is no special treatment of separate iteration rounds (they are all
merged together).

=head2 References & Information about the BLAST program

B<WEBSITES:>

   http://www.ncbi.nlm.nih.gov/BLAST/                 - Homepage at NCBI
   http://www.ncbi.nlm.nih.gov/BLAST/blast_help.html  - Help manual
   http://blast.wustl.edu/                            - WashU-Blast2

B<PUBLICATIONS:> (with PubMed links)

     Altschul S.F., Gish W., Miller W., Myers E.W., Lipman D.J. (1990).
     "Basic local alignment search tool", J Mol Biol 215: 403-410.

http://www.ncbi.nlm.nih.gov/htbin-post/Entrez/query?uid=2231712&form=6&db=m&Dopt=r

     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.

http://www.ncbi.nlm.nih.gov/htbin-post/Entrez/query?uid=9254694&form=6&db=m&Dopt=r

     Karlin, Samuel and Stephen F. Altschul (1990).  Methods  for
     assessing the statistical significance of molecular sequence
     features by using general scoring schemes. Proc. Natl. Acad.
     Sci. USA 87:2264-68.

http://www.ncbi.nlm.nih.gov/htbin-post/Entrez/query?uid=2315319&form=6&db=m&Dopt=b

     Karlin, Samuel and Stephen F. Altschul (1993).  Applications
     and statistics for multiple high-scoring segments in molecu-
     lar sequences. Proc. Natl. Acad. Sci. USA 90:5873-7.

http://www.ncbi.nlm.nih.gov/htbin-post/Entrez/query?uid=8390686&form=6&db=m&Dopt=b

=head1 USAGE

=head2 Creating Blast objects

A Blast object can be constructed from the contents of a Blast report
using a set of named parameters that specify significance criteria for
parsing.  The report data can be read in from an existing file
specified with the C<-file =E<gt> 'filename'> parameter or from a
STDIN stream containing potentially multiple Blast reports. If the
C<-file> parameter does not contain a valid filename, STDIN will be
used. Separate Blast objects will be created for each report in the
stream.

To parse the report, you must include a C<-parse =E<gt> 1> parameter
in addition to any other parsing parameters
See L<parse()|parse> for a full description of parsing parameters.
To run a new report and then parse it, include a C<-run =E<gt> \%runParams>
parameter containing a reference to a hash
that hold the parameters required by the L<run()|run> method.

The constructor for Blast objects is inherited from Bio::Tools::SeqAnal.pm.
See the B<_initialize>() method of that package for general information
relevant to creating Blast objects. (The B<new>() method, inherited from
B<Bio::Root::Object.pm>, calls B<_initialize>(). See L<Links to related modules>).

The Blast object can read compressed (gzipped) Blast report
files. Compression/decompression uses the gzip or compress programs
that are standard on Unix systems and should not require special
configuration. If you can't or don't want to use gzip as the file
compression tool, either pre-uncompress your files before parsing with
this module or modify B<Bio::Root::Utilities.pm> to your liking.

Blast objects can be generated either by direct instantiation as in:

 use Bio::Tools::Blast;		
 $blast = new Bio::Tools::Blast (%parameters);

=head2 Using the Static $Blast Object

 use Bio::Tools::Blast qw(:obj);		

This exports the static $Blast object into your namespace. "Static"
refers to the fact that it has class scope and there is one of these
created when you use this module. The static $Blast object is
basically an empty object that is provided for convenience and is also
used for various internal chores.

It is exported by this module and can be used for
parsing and running reports as well as HTML-formatting without having
to first create an empty Blast object.

Using the static $Blast object for parsing a STDIN stream of Blast reports:

    use Bio::Tools::Blast qw(:obj);

    sub process_blast {
	my $blastObj = shift;
	print $blastObj->table();
	$blastObj->destroy;
    }

    $Blast->parse( -parse     => 1,
		   -signif    => '1e-10',
		   -exec_func => \&process_blast,
		   );

Then pipe a stream of Blast reports into your script via STDIN.  For
each Blast report extracted from the input stream, the parser will
generate a new Blast object and pass it to the function specified by
C<-exec_func>.  The destroy() call tells Perl to free the memory
associated with the object, important if you are crunching through
many reports. This method is inherited from B<Bio::Root::Object.pm>
(see L<Links to related modules>). See L<parse()|parse> for a full
description of parameters and L<the DEMO SCRIPTS section | DEMO> for
additional examples.

=head2 Running Blasts

To run a Blast, create a new Blast object with a C<-run =E<gt>
\%runParams> parameter.  Remote Blasts are performed by including a
C<-method =E<gt> 'remote'> parameter; local Blasts are performed by
including a C<-method =E<gt> 'local'> parameter.  See 
L<Running Blast reports> as well as the
L<the DEMO SCRIPTS section | DEMO> for examples.  
Note that running local Blasts is not yet supported, see below.

Note that the C<-seqs =E<gt> [ $seqs ]> run parameter must contain a
reference to an array of B<Bio::Seq.pm> objects 
(L<Links to related modules>). Encapsulating the sequence in an 
object makes sequence information much easier to handle as it can 
be supplied in a variety of formats. Bio::Seq.pm is included with 
this distribution (L<INSTALLATION | INSTALLATION>).

Remote Blasts are implemented using the
B<Bio::Tools::Blast::Run::Webblast.pm> module.  Local Blasts require
that you customize the B<Bio::Tools::Blast::Run::LocalBlast.pm>
module.  The version of LocalBlast.pm included with this distribution
provides the basic framework for running local Blasts.
See L<Links to related modules>.

=head2 Significance screening

A C<-signif> parameter can be used to screen out all hits with
P-values (or Expect values) above a certain cutoff. For example, to
exclude all hits with Expect values above 1.0e-10: C<-signif =E<gt>
1e-10>. Providing a C<-signif> cutoff can speed up processing
tremendously, since only a small fraction of the report need be
parsed. This is because the C<-signif> value is used to screen hits
based on the data in the "Description" section of the Blast report:

For NCBI BLAST2 reports:

                                                                     Score     E
  Sequences producing significant alignments:                        (bits)  Value

  sp|P31376|YAB1_YEAST  HYPOTHETICAL 74.1 KD PROTEIN IN CYS3-MDM10...   957  0.0

For BLAST1 or WashU-BLAST2 reports:

                                                                       Smallest
                                                                         Sum
                                                                High  Probability
  Sequences producing High-scoring Segment Pairs:              Score  P(N)      N

  PDB:3PRK_E Proteinase K complexed with inhibitor ...........   504  1.8e-50   1

Thus, the C<-signif> parameter will screen based on Expect values for
BLAST2 reports and based on P-values for BLAST1/WashU-BLAST2 reports.

To screen based on other criteria, you can supply a C<-filt_func>
parameter containing a function reference that takes a
B<Bio::Tools::Sbjct.pm> object as an argument and returns a boolean,
true if the hit is to be screened out. See example below for
L<Screening hits using arbitrary criteria>.

=head2 Get the best hit.

     $hit = $blastObj->hit;

A "hit" is contained by a B<Bio::Tools::Blast::Sbjct.pm> object.

=head2 Get the P-value or Expect value of the most significant hit.

     $p = $blastObj->lowest_p;
     $e = $blastObj->lowest_expect;

Alternatively:

     $p = $blastObj->hit->p;
     $e = $blastObj->hit->expect;

Note that P-values are not reported in NCBI Blast2 reports.

=head2 Iterate through all the hits

     foreach $hit ($blastObj->hits) {
	 printf "%s\t %.1e\t %d\t %.2f\t %d\n",
	                  $hit->name, $hit->expect, $hit->num_hsps,
                          $hit->frac_identical, $hit->gaps;
     }

Refer to the documentation for B<Bio::Tools::Blast::Sbjct.pm>
for other ways to work with hit objects (L<Links to related modules>).

=head2 Screening hits using arbitrary criteria

    sub filter { $hit=shift;
		 return ($hit->gaps == 0 and
			 $hit->frac_conserved > 0.5); }

     $blastObj = Bio::Tools::Blast->new( -file      => '/tmp/blast.out',
					 -parse     => 1,
					 -filt_func => \&filter );

While the Blast object is parsing the report, each hit checked by calling
&filter($hit). All hits that generate false return values from &filter
are screened out and will not be added to the Blast object.
Note that the Blast object will normally stop parsing the report after
the first non-significant hit or the first hit that does not pass the
filter function. To force the Blast object to check all hits,
include a C<-check_all_hits =E<gt> 1>  parameter.
Refer to the documentation for B<Bio::Tools::Blast::Sbjct.pm>
for other ways to work with hit objects.

=over 4

=item Hit start, end coordinates.

      print $sbjct->start('query');
      print $sbjct->end('sbjct');

In array context, you can get information for both query and sbjct with one call:

      ($qstart, $sstart) = $sbjct->start();
      ($qend, $send)     = $sbjct->end();

For important information regarding coordinate information, see
the L<HSP start, end, and strand> section below.
Also check out documentation for the start and end methods in B<Bio::Tools::Blast::Sbjct.pm>,
which explains what happens if there is more than one HSP.

=back

=head2 Working with HSPs

=over 4

=item Iterate through all the HSPs of every hit

     foreach $hit ($blastObj->hits) {
	 foreach $hsp ($hit->hsps) {
	 printf "%.1e\t %d\t %.1f\t %.2f\t %.2f\t %d\t %d\n",
	                  $hsp->expect, $hsp->score, $hsp->bits,
                          $hsp->frac_identical, $hsp->frac_conserved,
	                  $hsp->gaps('query'), $hsp->gaps('sbjct');
     }

Refer to the documentation for B<Bio::Tools::Blast::HSP.pm>
for other ways to work with hit objects (L<Links to related modules>).

=back

=over 4

=item Extract HSP sequence data as strings or sequence objects

Get the first HSP of the first hit and the sequences
of the query and sbjct as strings.

      $hsp = $blast_obj->hit->hsp;
      $query_seq = $hsp->seq_str('query');
      $hsp_seq = $hsp->seq_str('sbjct');

Get the indices of identical and conserved positions in the HSP query seq.

      @query_iden_indices = $hsp->seq_inds('query', 'identical');
      @query_cons_indices = $hsp->seq_inds('query', 'conserved');

Similarly for the sbjct sequence.

      @sbjct_iden_indices = $hsp->seq_inds('sbjct', 'identical');
      @sbjct_cons_indices = $hsp->seq_inds('sbjct', 'conserved');

      print "Query in Fasta format:\n", $hsp->seq('query')->layout('fasta');
      print "Sbjct in Fasta format:\n", $hsp->seq('sbjct')->layout('fasta');

See the B<Bio::Seq.pm> package for more information about using these sequence objects
(L<Links to related modules>).

=back

=over 4

=item Create sequence alignment objects using HSP sequences

      $aln = $hsp->get_aln;
      print " consensus:\n", $aln->consensus();
      print $hsp->get_aln->layout('fasta');

      $ENV{READSEQ_DIR} = '/home/users/sac/bin/solaris';
      $ENV{READSEQ} = 'readseq';
      print $hsp->get_aln->layout('msf');

MSF formated layout requires Don Gilbert's ReadSeq program (not included).
See the B<Bio::UnivAln.pm> for more information about using these alignment objects
(L<Links to related modules>)'.

=back

=over 4

=item HSP start, end, and strand

To facilitate HSP processing, endpoint data for each HSP sequence are
normalized so that B<start is always less than end>. This affects TBLASTN
and TBLASTX HSPs on the reverse complement or "Minus" strand.

Some examples of obtaining start, end coordinates for HSP objects:

      print $hsp->start('query');
      print $hsp->end('sbjct');
      ($qstart, $sstart) = $hsp->start();
      ($qend, $send) = $hsp->end();

Strandedness of the HSP can be assessed using the strand() method
on the HSP object:

      print $hsp->strand('query');
      print $hsp->strand('sbjct');

These will return 'Minus' or 'Plus'.
Or, to get strand information for both query and sbjct with a single call:

      ($qstrand, $sstrand) = $hsp->strand();

=back

=head2 Report Generation

=over 4

=item Generate a tab-delimited table of all results.

     print $blastObj->table;
     print $blastObj->table(0);   # don't include hit descriptions.
     print $blastObj->table_tiled;

The L<table()|table> method returns data for each B<HSP> of each hit listed one per
line. The L<table_tiled()|table_tiled> method returns data for each B<hit, i.e., Sbjct>
listed one per line; data from multiple HSPs are combined after tiling to
reduce overlaps. See B<Bio::Tools::Blast::Sbjct.pm> for more information about
HSP tiling.  These methods generate stereotypical, tab-delimited data for each
hit of the Blast report. The output is suitable for importation into
spreadsheets or database tables. Feel free to roll your own table function if
you need a custom table.

For either table method, descriptions of each hit can be included if a
single, true argument is supplied (e.g., $blastObj-E<gt>table(1)). The description
will be added as the last field. This will significantly increase the size of
the table. Labels for the table columns can be obtained with L<table_labels()|table_labels>
and L<table_labels_tiled()|table_labels_tiled>.

=back

=over 4

=item Print a summary of the Blast report

     $blastObj->display();
     $blastObj->display(-show=>'hits');

L<display()|display> prints various statistics extracted from the Blast report
such as database name, database size, matrix used, etc. The
C<display(-show=E<gt>'hits')> call prints a non-tab-delimited table
attempting to line the data up into more readable columns. The output
generated is similar to L<table_tiled()|table_tiled>.

=back

=over 4

=item HTML-format an existing report

     use Bio::Tools::Blast qw(:obj);

     # Going straight from a non HTML report file to HTML output using
     # the static $Blast object exported by Bio::Tools::Blast.pm
     $Blast->to_html(-file   => '/usr/people/me/blast.output.txt',
		     -header => qq|<H1>BLASTP Results</H1><A HREF="home.html">Home</A>|
		     );

     # You can also use a specific Blast object created previously.
     $blastObj->to_html;

L<to_html()|to_html> will send HTML output, line-by-line, directly to STDOUT
unless an C<-out =E<gt> array_ref> parameter is supplied (e.g., C<-out
=E<gt> \@array>), in which case the HTML will be stored in @array, one
line per array element.  The direct outputting permits faster response
time since Blast reports can be huge. The -header tag can contain a
string containing any HTML that you want to appear at the top of the
Blast report.

=back

=head1 DEMO SCRIPTS

Sample Scripts are included in the central bioperl distribution in the
'examples/blast/' directory (see L<INSTALLATION | INSTALLATION>):

=head2 Handy library for working with Bio::Tools::Blast.pm

   examples/blast/blast_config.pl

=head2 Parsing Blast reports one at a time.

   examples/blast/parse_blast.pl
   examples/blast/parse_blast2.pl
   examples/blast/parse_positions.pl

=head2 Parsing sets of Blast reports.

   examples/blast/parse_blast.pl
   examples/blast/parse_multi.pl

   B<Warning:> See note about L<Memory Usage Issues>.

=head2 Running Blast analyses one at a time.

   examples/blast/run_blast_remote.pl

=head2 Running Blast analyses given a set of sequences.

   examples/blast/blast_seq.pl

=head2 HTML-formatting Blast reports.

   examples/blast/html.pl

=head1 TECHNICAL DETAILS

=head2 Blast Modes

A BLAST object may be created using one of three different modes as
defined by the B<Bio::Tools::SeqAnal.pm> package 
(See L<Links to related modules>):

 -- parse - Load a BLAST report and parse it, storing parsed data in
    Blast.pm object.
 -- run      - Run the BLAST program to generate a new report.
 -- read     - Load a BLAST report into the Blast object without parsing.

B<Run mode support has recently been added>.  The module
B<Bio::Tools::Blast::Run::Webblast.pm> is an modularized adaptation of
the webblast script by Alex Dong Li:

   http://www.genet.sickkids.on.ca/bioinfo_resources/software.html#webblast

for running remote Blast analyses and saving the results locally.  Run
mode can be combined with a parse mode to generate a Blast report and
then build the Blast object from the parsed results of this report
(see L<run()|run> and L<SYNOPSIS | SYNOPSIS>).

In read mode, the BLAST report is read in by the Blast object but is
not parsed.  This could be used to internalize a Blast report but not
parse it for results (e.g., generating HTML formatted output).

=head2 Significant Hits

This module permits the screening of hits on the basis of
user-specified criteria for significance. Currently, Blast reports can
be screened based on:

   CRITERIA                            PARAMETER       VALUE
   ----------------------------------  ---------      ----------------
  1) the best Expect (or P) value      -signif        float or sci-notation
  2) the length of the query sequence  -min_length    integer
  3) arbitrary criteria                -filt_func     function reference

The parameters are used for construction of the BLAST object or when
running the L<parse()|parse> method on the static $Blast object.  The
-SIGNIF value represents the number listed in the description section
at the top of the Blast report. For Blast2, this is an Expect value,
for Blast1 and WashU-Blast2, this is a P-value.  The idea behind the
C<-filt_func> parameter is that the hit has to pass through a filter
to be considered significant. Refer to the documentation for
B<Bio::Tools::Blast::Sbjct.pm> for ways to work with hit objects.

Using a C<-signif> parameter allows for the following:

=over 2

=item Faster parsing.

Each hit can be screened by examination of the description line alone
without fully parsing the HSP alignment section.

=item Flexibility.

The C<-signif> tag provides a more semantic-free way to specify the
value to be used as a basis for screening hits. Thus, C<-signif> can
be used for screening Blast1 or Blast2 reports. It is up to the user
to understand whether C<-signif> represents a P-value or an Expect
value.

=back

Any hit not meeting the significance criteria will not be added to the
"hit list" of the BLAST object. Also, a BLAST object without any hits
meeting the significance criteria will throw an exception during
object construction (a fatal event).

=head2 Statistical Parameters

There are numerous parameters which define the behavior of the BLAST
program and which are useful for interpreting the search
results. These parameters are extracted from the Blast report:

  filter  --  for masking out low-complexity sequences or short repeats
  matrix  --  name of the substitution scoring matrix (e.g., BLOSUM62)
  E       --  Expect filter (screens out frequent scores)
  S       --  Cutoff score for segment pairs
  W       --  Word length
  T       --  Threshold score for word pairs
  Lambda, --  Karlin-Altschul "sum" statistical parameters dependent on
   K, H        sequence composition.
  G       --  Gap creation penalty.
  E       --  Gap extension penalty.

These parameters are not always needed. Extraction may be turned off
explicitly by including a C<-stats =E<gt> 0> parameter during object
construction.  Support for all statistical parameters is not complete.

For more about the meaning of parameters, check out the NCBI URLs given above.

=head2 Module Organization

The modules that comprise this Bioperl Blast distribution are location in the
Bio:: hierarchy as shown in the diagram below.

                            Bio/
                             |
               +--------------------------+
               |                          |
          Bio::Tools                  Bio::Root
               |                          |
    +----------------------+           Object.pm
    |          |           |
 SeqAnal.pm  Blast.pm    Blast/
                           |
            +---------+---------+------------+
            |         |         |            |
          Sbjct.pm   HSP.pm   HTML.pm       Run/
                                             |
                                       +------------+
                                       |            |
                                  Webblast.pm   LocalBlast.pm

Bio::Tools::Blast.pm is a concrete class that inherits from
B<Bio::Tools::SeqAnal.pm> and relies on other modules for parsing and
managing BLAST data.  Worth mentioning about this hierarchy is the
lack of a "Parse.pm" module.  Since parsing is considered central to
the purpose of the Bioperl Blast module (and Bioperl in general), it
seems somewhat unnatural to segregate out all parsing code. This
segregation could also lead to inefficiencies and harder to maintain
code. I consider this issue still open for debate.

Bio::Tools::Blast.pm, B<Bio::Tools::Blast::Sbjct.pm>, and
B<Bio::Tools::Blast::HSP.pm> are mostly dedicated to parsing and all
can be used to instantiate objects.  Blast.pm is the main "command and
control" module, inheriting some basic behaviors from SeqAnal.pm
(things that are not specific to Blast I<per se>).

B<Bio::Tools::Blast::HTML.pm> contains functions dedicated to
generating HTML-formatted Blast reports and does not generate objects.

=head2 Running Blasts: Details

B<Bio::Tools::Blast::Run::Webblast.pm> contains a set of functions for
running Blast analyses at a remote server and also does not
instantiate objects.  It uses a helper script called postclient.pl,
located in the Run directory.  The proposed LocalBlast.pm module would
be used for running Blast reports on local machines and thus would be
customizable for different sites. It would operate in a parallel
fashion to Webblast.pm (i.e., being a collection of functions, taking
in sequence objects or files, returning result files).

The Run modules are considered experimental. In particular,
Webblast.pm catures an HTML-formatted version of the Blast report from
the NCBI server and strips out the HTML in preparation for parsing. A
more direct approach would be to capture the Blast results directly
from the server using an interface to the NCBI toolkit.  This approach
was recently proposed on the Bioperl mailing list:
http://www.uni-bielefeld.de/mailinglists/BCD/vsns-bcd-perl/9805/0000.html

=head2 Memory Usage Issues

Parsing large numbers of Blast reports (a few thousand or so) with
Bio::Tools::Blast.pm may lead to unacceptable memory usage situations.
This is somewhat dependent of the size and complexity of the reports.

While this problem is under investigation, here are some workarounds
that fix the memory usage problem:

=over 4

=item 1 Don't specify a -signif criterion when calling L<parse()|parse>.

The C<-signif> value is used for imposing a upper limit to the expect- or
P-value for Blast hits to be parsed. For reasons that are still under
investigation, specifying a value for C<-signif> in the L<parse()|parse>
method prevents Blast objects from being fully
garbage collected. When using the B<parse_blast.pl> or B<parse_multi.pl>
scripts in C<examples/blast/> of the bioperl distribution), don't supply
a C<-signif> command-line parameter.

=item 2 If you want to impose a -signif criterion, put it inside a
-filt_func.

For the L<parse()|parse> method, a -signif =E<gt> 1e-5 parameter is equivalent
to using a filter function parameter of

 -filt_func => sub { my $hit = shift; return $hit->signif <= 1e-5; }

Using the B<examples/blast/parse_multi.pl> script, you can supply a
command-line argument of

 -filt_func '$hit->signif <= 1e-5'

For more information, see L<parse()|parse> and the section 
L<Screening hits using arbitrary criteria>.

=back

=head1 TODO

=over 4

=item * Develop a functional, prototype Bio::Tools::Blast::Run::LocalBlast.pm module.

=item * Add support for PSI-BLAST and PHI-BLAST

=item * Parse histogram of expectations and retrieve gif image in
Blast report (if present).

=item * Further investigate memory leak that occurs when parsing Blast
streams whe supplying a -signif parameter to L<parse()|parse>.

=item * Access Blast results directly from the NCBI server using a
Perl interface to the NCBI toolkit or XML formated Blast reports (when
available).

=item * Further exploit Bio::UnivAln.pm and multiple-sequence
alignment programs using HSP sequence data. Some of this may best go
into a separate, dedicated module or script as opposed to burdening
Blast.pm, Sbjct.pm, and HSP.pm with additional functionality that is
not always required.

=item * Add an example script for parsing Blast reports containing
HTML formatting.

=back

=head1 VERSION

Bio::Tools::Blast.pm, 0.09

=head1 FEEDBACK

=head2 Mailing Lists

User feedback is an integral part of the evolution of this and other
Bioperl modules.  Send your comments and suggestions preferably to one
of the Bioperl mailing lists.  Your participation is much appreciated.

    [e-mail]             - General discussion
    http://bio.perl.org/MailList.html - About the mailing lists

=head2 Reporting Bugs

Report bugs to the Bioperl bug tracking system to help us keep track
the bugs and their resolution. Bug reports can be submitted via email
or the web:

    [e-mail]
    http://bugzilla.bioperl.org/

=head1 AUTHOR

Steve Chervitz, [e-mail]

See the L<FEEDBACK | FEEDBACK> section for where to send bug reports and comments.

=head1 ACKNOWLEDGEMENTS

This module was developed under the auspices of the Saccharomyces Genome
Database:
    http://genome-www.stanford.edu/Saccharomyces

Other contributors include: Alex Dong Li (webblast), Chris Dagdigian
(Seq.pm), Steve Brenner (Seq.pm), Georg Fuellen (Seq.pm, UnivAln.pm),
and untold others who have offered comments (noted in the
Bio/Tools/Blast/CHANGES file of the distribution).

=head1 COPYRIGHT

Copyright (c) 1996-98 Steve Chervitz. All Rights Reserved.  This
module is free software; you can redistribute it and/or modify it
under the same terms as Perl itself.

=head1 SEE ALSO

 Bio::Tools::SeqAnal.pm                  - Sequence analysis object base class.
 Bio::Tools::Blast::Sbjct.pm             - Blast hit object.
 Bio::Tools::Blast::HSP.pm               - Blast HSP object.
 Bio::Tools::Blast::HTML.pm              - Blast HTML-formating utility class.
 Bio::Tools::Blast::Run::Webblast.pm     - Utility module for running Blasts remotely.
 Bio::Tools::Blast::Run::LocalBlast.pm   - Utility module for running Blasts locally.
 Bio::Seq.pm                             - Biosequence object
 Bio::UnivAln.pm                         - Biosequence alignment object.
 Bio::Root::Object.pm                    - Proposed base class for all Bioperl objects.

=head2 Links to related modules

 Bio::Tools::SeqAnal.pm
      http://bio.perl.org/Core/POD/Bio/Tools/SeqAnal.html

 Bio::Tools::Blast::Sbjct.pm
      http://bio.perl.org/Core/POD/Bio/Tools/Blast/Sbjct.html

 Bio::Tools::Blast::HSP.pm
      http://bio.perl.org/Core/POD/Bio/Tools/Blast/HSP.html

 Bio::Tools::Blast::HTML.pm
      http://bio.perl.org/Core/POD/Bio/Tools/Blast/HTML.html

 Bio::Tools::Blast::Run::Webblast.pm
      http://bio.perl.org/Core/POD/Bio/Tools/Blast/Run/Webblast.html

 Bio::Tools::Blast::Run::LocalBlast.pm
      http://bio.perl.org/Core/POD/Bio/Tools/Blast/Run/LocalBlast.html

 Bio::Seq.pm
      http://bio.perl.org/Core/POD/Seq.html

 Bio::UnivAln.pm
      http://bio.perl.org/Projects/SeqAlign/
      Europe:  http://www.techfak.uni-bielefeld.de/bcd/Perl/Bio/#univaln

 Bio::Root::Object.pm
      http://bio.perl.org/Core/POD/Root/Object.html

 http://bio.perl.org/Projects/modules.html  - Online module documentation
 http://bio.perl.org/Projects/Blast/        - Bioperl Blast Project
 http://bio.perl.org/                       - Bioperl Project Homepage

L<References & Information about the BLAST program>.

=head1 KNOWN BUGS

There is a memory leak that occurs when parsing parsing streams
containing large numbers of Blast reports (a few thousand or so) and
specifying a -signif parameter to the L<parse()|parse> method. For a
workaround, see L<Memory Usage Issues>.

Not sharing statistical parameters between different Blast objects
when parsing a multi-report stream has not been completely tested and
may be a little buggy.

Documentation inconsistencies or inaccuracies may exist since this
module underwend a fair bit of re-working going from 0.75 to 0.80
(corresponds to versions 0.04.4 to 0.05 of the bioperl distribution).

=cut

#
##
###
#### END of main POD documentation.
###
##
#

=head1 APPENDIX

Methods beginning with a leading underscore are considered private and
are intended for internal use by this module. They are B<not>
considered part of the public interface and are described here for
documentation purposes only.

=cut

##############################################################################
##                          CONSTRUCTOR                                     ##
##############################################################################


sub new { 
	my ($class,@args) = @_;
	$class->warn("Bio::Tools::BLAST is deprecated, use Bio::SearchIO system or Bio::Tools::BPlite");
return	$class->SUPER::new(@args);
 }

## The Blast.pm object relies on the the superclass constructor:
## Bio::Tools::SeqAnal::_initialize(). See that module for details.

#-------------
sub destroy {
#-------------
    my $self=shift;
    $DEBUG==2 && print STDERR "DESTROYING $self ${\$self->name}";
    if($self->{'_hits'}) {
	foreach($self->hits) {
	    $_->destroy;
	    undef $_;
	}
	undef $self->{'_hits'};
        #$self->{'_hits'}->remove_all;  ## When and if this member becomes a vector.
    }

    $self->SUPER::destroy;
}

#####################################################################################
##                                  ACCESSORS                                      ##
#####################################################################################

=head2 run

 Usage     : $object->run( %named_parameters )
 Purpose   : Run a local or remote Blast analysis on one or more sequences.
 Returns   : String containing name of Blast output file if a single Blast
           : is run.
           :  -- OR --
           : List of Blast objects if multiple Blasts are being run as a group.
 Argument  : Named parameters:  (PARAMETER TAGS CAN BE UPPER OR LOWER CASE).
           :    -METHOD  => 'local' or 'remote' (default = remote),
           :    -PARSE   => boolean, (true if the results are to be parsed after the run)
           :    -STRICT  => boolean, the strict mode to use for the resulting Blast objects.
           :  ADDITIONAL PARAMETERS:
           :      See methods _run_remote() and _run_local() for required
           :      parameters necessary for running the blast report.
 Throws    : Exception if no Blast output file was obtained.
 Comments  : This method is called automatically during construction of a
           : Blast.pm object when run parameters are sent to the constructor:
           :  $blastObj = new Bio::Tools::Blast (-RUN =>\%runParam,
	   :					 %parseParam );
           :
           : The specific run methods (local or remote) called by run()
           : must return a list containing  the file name(s) with the Blast output.
           :
           : The run() method can perform single or multiple Blast runs
           : (analogous to the way parse() works) depending on how many
           : sequences are submitted. However, the running of multiple
           : Blasts is probably better handled at the script level. See notes in
           : the "TODO" section below.
           :
           : As for what to do with the Blast result file, that decision is
           : left for the user who can direct the Blast object to delete, compress,
           : or leave it alone.
           :
           : This method does not worry about load balancing, which
           : is probably best handled at the server level.
           :
 TODO:     : Support for running+parsing multiple Blast analyses with a
           : single run() call is incomplete. One can generate multiple
           : reports by placing more than one sequence object in the -seqs
           : reference parameter. This saves some overhead in the code
           : that executes the Blasts since all options are configured once.
           : (This is analogous to parsing using the static $Blast object
           : see parse() and _parse_stream()).
           :
           : The trouble is that Blast objects for all runs are constructed,
           : parsed (if necessary), and then returned as a group
           : This can require lots of memory when run+parsing many Blasts
           : but should be fine if you just want to run a bunch Blasts.
           :
           : For now, when running+parsing Blasts, stick to running one
           : Blast at a time, building the Blast object with the results
           : of that report, and processing as necessary.
           :
           : Support for running PSI-Blast is not complete.

See Also:  L<_run_remote()|_run_remote>, L<_run_local()|_run_local>, L<parse()|parse>

=cut

#---------
sub run {
#---------
    my ($self, %param) = @_;
    my($method, $parse, $strict) =
	$self->_rearrange([qw(METHOD PARSE STRICT)], %param);

    $strict = $self->strict($strict) if $strict;

    my (@files);
    if($method =~ /loc/i) {
	@files = $self->_run_local(%param);
	
    } else {
	@files = $self->_run_remote(%param);
    }

    $self->throw("Run Blast failed: no Blast output created.") if !@files;

    if(scalar(@files) == 1) {
	# If there was just one Blast output file, prepare to incorporate it
	# into the current Blast object. run() is called before parse() in the
	# SeqAnal.pm constructor.
	if($files[0] ne 'email') {
	    $self->file($files[0]);
	} else {
	    # Can't do anything with the report.
	    $self->throw("Blast report to be sent via e-mail.");
	}

    } else {
	# If there are multiple report files, build individual Blast objects foreach.
	# In this situation, the static $Blast object is being used to run
	# a set of related Blasts, similar to the way parse() can be used.
	# This strategy is not optimal since all reports are generated first
	# before any are parsed.
	# Untested.

	my(@objs);
	foreach(@files) {
	    push @objs, new Bio::Tools::Blast(-FILE   => $_,
					      -PARSE  => $parse || 0,
					      -STRICT => $strict,
					      );
	}
	return @objs;
    }
}

=head2 _run_remote

 Usage     : n/a; internal method called by run()
           : $object->_run_remote( %named_parameters )
 Purpose   : Run Blast on a remote server.
 Argument  : Named parameters:
           :   See documentation for function &blast_remote in
           :   Bio::Tools::Blast::Run::Webblast.pm for description
           :   of parameters.
 Comments  : This method requires the Bio::Tools::Blast::Run::Webblast.pm
           : which conforms to this minimal API:
           :    * export a method called &blast_remote that accepts a
           :      Bio::Tools::Blast.pm object + named parameters
           :      (specified in the Webblast.pm module).
           :    * return a list of names of files containing the raw Blast reports.
           :      (When building a Blast object, this list would contain a
           :       single file from which the Blast object is to be constructed).

See Also   : L<run()|run>, L<_run_local()|_run_local>, B<Bio::Tools::Blast::Run::Webblast.pm::blast_remote()>, L<Links to related modules>

=cut

#----------------
sub _run_remote {
#----------------
    my ($self, %param) = @_;

    require Bio::Tools::Blast::Run::Webblast;
    Bio::Tools::Blast::Run::Webblast->import(qw(&blast_remote));

    &blast_remote($self, %param);
}

=head2 _run_local

 Usage     : n/a; internal method called by run()
           : $object->_run_local(%named_parameters)
 Purpose   : Run Blast on a local machine.
 Argument  : Named parameters:
           :   See documentation for function &blast_local in
           :   Bio::Tools::Blast::Run::LocalBlast.pm for description
           :   of parameters.
 Comments  : This method requires the Bio::Tools::Blast::Run::LocalBlast.pm
           : module which should be customized for your site. This module would
           : contain all the commands, paths, environment variables, and other
           : data necessary to run Blast commands on a local machine, but should
           : not contain any semantics for specific query sequences.
           :
           : LocalBlast.pm should also conform to this minimal API:
           :    * export a method called &blast_local that accepts a
           :       Bio::Tools::Blast.pm object + named parameters
           :      (specified in the LocalBlast.pm module).
           :    * return a list of names of files containing the raw Blast reports.
           :      (When building a Blast object, this list would contain a
           :       single file from which the Blast object is to be constructed).

See Also   : L<run()|run>, L<_run_remote()|_run_remote>, B<Bio::Tools::Blast::Run::LocalBlast::blast_local()>, L<Links to related modules>

=cut

#--------------
sub _run_local {
#--------------
    my ($self, %param) = @_;

    require Bio::Tools::Blast::Run::Webblast;
    Bio::Tools::Blast::Run::Webblast->import(qw(&blast_local));

    &blast_local($self, %param);
}

=head2 db_remote

 Usage     : @dbs = $Blast->db_remote( [seq_type] );
 Purpose   : Get a list of available sequence databases for remote Blast analysis.
 Returns   : Array of strings
 Argument  : seq_type = 'p' or 'n'
           :  'p' = Gets databases for peptide searches  (default)
           :  'n' = Gets databases for nucleotide searches
 Throws    : n/a
 Comments  : Peptide databases are a subset of the nucleotide databases.
           : It is convenient to call this method on the static $Blast object
           : as shown in Usage.

See Also   : L<db_local()|db_local>

=cut

#----------------
sub db_remote {
#----------------
    my ($self, $type) = @_;
    $type ||= 'p';

    require Bio::Tools::Blast::Run::Webblast;
    Bio::Tools::Blast::Run::Webblast->import(qw(@Blast_dbp_remote
						 @Blast_dbn_remote));

    # We shouldn't have to fully qualify the Blast_dbX_remote arrays. Hm.

    my(@dbs);
    if( $type =~ /^p|amino/i) {
	@dbs = @Bio::Tools::Blast::Run::Webblast::Blast_dbp_remote;
    } else {
	@dbs = @Bio::Tools::Blast::Run::Webblast::Blast_dbn_remote;
    }
    @dbs;
}

=head2 db_local

 Usage     : @dbs = $Blast->db_local( [seq_type] );
 Purpose   : Get a list of available sequence databases for local Blast analysis.
 Returns   : Array of strings
 Argument  : seq_type = 'p' or 'n'
           :  'p' = Gets databases for peptide searches  (default)
           :  'n' = Gets databases for nucleotide searches
 Throws    : n/a
 Comments  : Peptide databases are a subset of the nucleotide databases.
           : It is convenient to call this method on the static $Blast object.
             as shown in Usage.

See Also   : L<db_remote()|db_remote>

=cut

#----------------
sub db_local {
#----------------
    my ($self, $type) = @_;
    $type ||= 'p';

    require Bio::Tools::Blast::Run::LocalBlast;
    Bio::Tools::Blast::Run::LocalBlast->import(qw(@Blast_dbp_local
						  @Blast_dbn_local));

    # We shouldn't have to fully qualify the Blast_dbX_local arrays. Hm.

    my(@dbs);
    if( $type =~ /^p|amino/i) {
	@dbs = @Bio::Tools::Blast::Run::LocalBlast::Blast_dbp_local;
    } else {
	@dbs = @Bio::Tools::Blast::Run::LocalBlast::Blast_dbn_local;
    }
    @dbs;
}

=head2 parse

 Usage     : $blast_object->parse( %named_parameters )
 Purpose   : Parse a Blast report from a file or STDIN.
           :   * Parses a raw BLAST data, populating Blast object with report data.
           :   * Sets the significance cutoff.
           :   * Extracts statistical parameters about the BLAST run.
           :   * Handles both single files and streams containing multiple reports.
 Returns   : integer (number of Blast reports parsed)
 Argument  : <named parameters>:  (PARAMETER TAGS CAN BE UPPER OR LOWER CASE).
	   : -FILE       => string (name of file containing raw Blast output.
           :                         Optional. If a valid file is not supplied,
	   :		             STDIN will be used).
	   : -SIGNIF     => number (float or scientific notation number to be used
	   :                         as a P- or Expect value cutoff;
	   :			     default =  $DEFAULT_SIGNIF (999)).
	   : -FILT_FUNC  => func_ref (reference to a function to be used for
           :                          filtering out hits based on arbitrary criteria.
           :                          This function should take a
           :                          Bio::Tools::Blast::Sbjct.pm object as its first
           :                          argument and return a boolean value,
	   :                          true if the hit should be filtered out).
           :                          Sample filter function:
           :                          -FILT_FUNC => sub { $hit = shift;
	   :				                  $hit->gaps == 0; },
           : -CHECK_ALL_HITS => boolean (check all hits for significance against
           :                             significance criteria.  Default = false.
	   :			         If false, stops processing hits after the first
           :                             non-significant hit or the first hit that fails
           :                             the filt_func call. This speeds parsing,
           :                             taking advantage of the fact that the hits
           :                             are processed in the order they are ranked.)
           : -MIN_LEN     => integer (to be used as a minimum query sequence length
           :                          sequences below this length will not be processed).
           :                          default = no minimum length).
	   : -STATS       => boolean (collect stats for report: matrix, filters, etc.
           :                          default = false).
	   : -BEST        => boolean (only process the best hit of each report;
           :                          default = false).
           : -OVERLAP     => integer (the amount of overlap to permit between
           :                          adjacent HSPs when tiling HSPs,
           :                          Default = $MAX_HSP_OVERLAP (2))
           :
           : PARAMETERS USED WHEN PARSING MULTI-REPORT STREAMS:
           : --------------------------------------------------
	   : -SHARE       => boolean (set this to true if all reports in stream
	   :			      share the same stats. Default = true)
           :                          Must be set to false when parsing both Blast1 and
           :                          Blast2 reports in the same run or if you need
           :                          statistical params for each report, Lambda, K, H).
	   : -STRICT      => boolean (use strict mode for all Blast objects created.
           :                          Increases sensitivity to errors. For single
           :                          Blasts, this is parameter is sent to new().)
           : -EXEC_FUNC   => func_ref (reference to a function for processing each
           :                           Blast object after it is parsed. Should accept a
           :                           Blast object as its sole argument. Return value
           :                           is ignored. If an -EXEC_FUNC parameter is supplied,
           :                           the -SAVE_ARRAY parameter will be ignored.)
           : -SAVE_ARRAY  =>array_ref, (reference to an array for storing all
           :                            Blast objects as they are created.
           :                            Experimental. Not recommended.)
           : -SIGNIF_FMT  => boolean   String of 'exp' or 'parts'. Sets the format
           :                           for reporting P/Expect values. 'exp' reports
           :                           only the exponent portion. 'parts' reports
           :                           them as a 2 element list. See signif_fmt()..
           :
 Throws    : Exception if BLAST report contains a FATAL: error.
           : Propagates any exception thrown by read().
           : Propagates any exception thrown by called parsing methods.
 Comments  : This method can be called either directly using the static $Blast object
           : or indirectly (by Bio::Tools::SeqAnal.pm) during constuction of an
           : individual Blast object.
           :
           : HTML-formatted reports can be parsed as well. No special flag is required
           : since it is detected automatically. The presence of HTML-formatting
           : will result in slower performace, however, since it must be removed
           : prior to parsing. Parsing HTML-formatted reports is highly
           : error prone and is generally not recommended.
           :
           : If one has an HTML report, do NOT remove the HTML from it by using the
           : "Save As" option of a web browser to save it as text. This renders the
           : report unparsable.
           : HTML-formatted reports can be parsed after running through the strip_html
           : function of Blast::HTML.pm as in:
           :    require Bio::Tools::Blast::HTML;
           :    Bio::Tools::Blast::HTML->import(&strip_html);
           :    &strip_html(\$data);
           :    # where data contains full contents of an HTML-formatted report.
           : TODO: write a demo script that does this.

See Also   : L<_init_parse_params()|_init_parse_params>, L<_parse_blast_stream()|_parse_blast_stream>, L<overlap()|overlap>, L<signif_fmt()|signif_fmt>, B<Bio::Root::Object::read()>, B<Bio::Tools::Blast::HTML.pm::strip_html()>, L<Links to related modules>

=cut

#---------
sub parse {
#---------
# $self might be the static $Blast object.
    my ($self, @param) = @_;

    my($signif, $filt_func, $min_len, $check_all, $overlap, $stats,
       $share, $strict, $best, $signif_fmt, $no_aligns) =
	$self->_rearrange([qw(SIGNIF FILT_FUNC MIN_LEN CHECK_ALL_HITS
			      OVERLAP STATS SHARE STRICT
			      BEST EXPONENT NO_ALIGNS )], @param);

    ## Initialize the static Blast object with parameters that
    ## apply to all Blast objects within a parsing session.

    &_init_parse_params($share, $filt_func, $check_all,
			$signif, $min_len, $strict,
			$best, $signif_fmt, $stats, $no_aligns
		       );

    my $count = $self->_parse_blast_stream(@param);
	
#    print STDERR "\nDONE PARSING STREAM.\n";

    if($Blast->{'_blast_errs'}) {
      my @errs = @{$Blast->{'_blast_errs'}};
      printf STDERR "\n*** %d BLAST REPORTS HAD FATAL ERRORS:\n", scalar(@errs);
      foreach(@errs) { print STDERR "$_\n"; }
      @{$Blast->{'_blast_errs'}} = ();
    }

    return $count;
}

=head2 _init_parse_params

 Title   : _init_parse_params
 Usage   : n/a; called automatically by parse()
 Purpose : Initializes parameters used during parsing of Blast reports.
         : This is a static method used by the $Blast object.
         : Calls _set_signif().
 Example :
 Returns : n/a
 Args    : Args extracted by parse().

See Also: L<parse()|parse>, L<_set_signif()|_set_signif>

=cut

#----------------------
sub _init_parse_params {
#----------------------
    my ($share, $filt_func, $check_all,
	$signif, $min_len, $strict,
	$best, $signif_fmt, $stats, $no_aligns) = @_;

    ## Default is to share stats.
    $Blast->{'_share'}  = defined($share) ? $share : 1;
    $Blast->{'_filt_func'} = $filt_func || 0;
    $Blast->{'_check_all'} = $check_all || 0;
    $Blast->{'_signif_fmt'} ||= $signif_fmt || '';
    $Blast->{'_no_aligns'} = $no_aligns || 0;

    &_set_signif($signif, $min_len, $filt_func);
    $Blast->strict($strict) if defined $strict;
    $Blast->best($best) if $best;
    $Blast->{'_blast_count'} = 0;

    ## If $stats is false, miscellaneous statistical and other parameters
    ## are NOT extracted from the Blast report (e.g., matrix name, filter used, etc.).
    ## This can speed processing when crunching tons of Blast reports.
    ## Default is to NOT get stats.
    $Blast->{'_get_stats'} = defined($stats) ? $stats : 0;

    # Clear any errors from previous parse.
    undef $Blast->{'_blast_errs'};
}

=head2 _set_signif

 Usage     : n/a; called automatically by _init_parse_params()
           : This is now a "static" method used only by $Blast.
           : _set_signif($signif, $min_len, $filt_func);
 Purpose   : Sets significance criteria for the BLAST object.
 Argument  : Obligatory three arguments:
           :   $signif = float or sci-notation number or undef
           :   $min_len = integer or undef
           :   $filt_func = function reference or undef
           :
           :   If $signif is undefined, a default value is set
           :   (see $DEFAULT_SIGNIF; min_length = not set).
 Throws    : Exception if significance value is defined but appears
           :   out of range or invalid.
           : Exception if $filt_func if defined and is not a func ref.
 Comments  : The significance of a BLAST report can be based on
           : the P (or Expect) value and/or the length of the query sequence.
           : P (or Expect) values GREATER than '_significance' are not significant.
           : Query sequence lengths LESS than '_min_length' are not significant.
           :
           : Hits can also be screened using arbitrary significance criteria
           : as discussed in the parse() method.
           :
           : If no $signif is defined, the '_significance' level is set to
           : $Bio::Tools::Blast::DEFAULT_SIGNIF (999).

See Also   : L<signif()|signif>, L<min_length()|min_length>, L<_init_parse_params()|_init_parse_params>, L<parse()|parse>

=cut

#-----------------
sub _set_signif {
#-----------------
    my( $sig, $len, $func ) = @_;

    if(defined $sig) {
	$Blast->{'_confirm_significance'} = 1;
	if( $sig =~ /[^\d.e-]/ or $sig <= 0) {
	    $Blast->throw("Invalid significance value: $sig",
			 "Must be greater than zero.");
	}
	$Blast->{'_significance'} = $sig;
    } else {
	$Blast->{'_significance'}   = $DEFAULT_SIGNIF;
	$Blast->{'_check_all'}      = 1 if not $Blast->{'_filt_func'};
    }

    if(defined $len) {
	if($len =~ /\D/ or $len <= 0) {
	    $Blast->warn("Invalid minimum length value: $len",
			"Value must be an integer > 0. Value not set.");
	} else {
	    $Blast->{'_min_length'} = $len;
	}
    }

    if(defined $func) {
	$Blast->{'_confirm_significance'} = 1;
	if($func and not ref $func eq 'CODE') {
	    $Blast->throw("Not a function reference: $func",
			  "The -filt_func parameter must be function reference.");
	  }
      }
  }

=head2 _parse_blast_stream

 Usage     : n/a. Internal method called by parse()
 Purpose   : Obtains the function to be used during parsing and calls read().
 Returns   : Integer (the number of blast reports read)
 Argument  : Named parameters  (forwarded from parse())
 Throws    : Propagates any exception thrown by _get_parse_blast_func() and read().

See Also   : L<_get_parse_blast_func()|_get_parse_blast_func>, B<Bio::Root::Object::read()>

=cut

#----------------------
sub _parse_blast_stream {
#----------------------
    my ($self, %param) = @_;

    my $func = $self->_get_parse_blast_func(%param);
#    my $func = sub { my $data = shift;
#		      printf STDERR "Chunk length = %d\n", length($data);
#		      sleep(3);
#		    };

    # Only setting the newline character once per session.
    $Newline ||= $Util->get_newline(-client => $self, %param);

    $self->read(-REC_SEP  =>"$Newline>",
		-FUNC     => $func,
		%param);

    return $Blast->{'_blast_count'};
}

=head2 _get_parse_blast_func

 Usage     : n/a; internal method used by _parse_blast_stream()
           : $func_ref = $blast_object->_get_parse_blast_func()
 Purpose   : Generates a function ref to be used as a closure for parsing
           : raw data as it is being loaded by Bio::Root::IOManager::read().
 Returns   : Function reference (closure).
 Comments  : The the function reference contains a fair bit of logic
           : at present. It could perhaps be split up into separate
           : functions to make it more 'digestible'.

See Also   : L<_parse_blast_stream()|_parse_blast_stream>

=cut

#--------------------------
sub _get_parse_blast_func {
#--------------------------
    my ($self, @param) = @_;

    my ($save_a, $exec_func) =
	$self->_rearrange([qw(SAVE_ARRAY EXEC_FUNC)], @param);

#    $MONITOR && print STDERR "\nParsing Blast stream (5/dot, 250/line)\n";
    my $count = 0;
    my $strict = $self->strict();

    # Some parameter validation.
    # Remember, all Blast parsing will use this function now.
    # You won't need a exec-func or save_array when just creating a Blast object
    #  as in: $blast = new Bio::Tools::Blast();
    if($exec_func and not ref($exec_func) eq 'CODE') {
	$self->throw("The -EXEC_FUNC parameter must be function reference.",
		    "exec_func = $exec_func");

    } elsif($save_a and not ref($save_a) eq 'ARRAY') {
	$self->throw("The -SAVE_ARRAY parameter must supply an array reference".
		     "when not using an -EXEC_FUNC parameter.");
    }

    ## Might consider breaking this closure up if possible.

     return sub {
	my ($data) = @_;
	## $data should contain one of three possible fragment types
        ## from a Blast report:
        ##   1. Header with description section,
        ##   2. An alignment section for a single hit, or
        ##   3. The final alignment section plus the footer section.
	## (record separator = "Newline>").

#	print STDERR "\n(BLAST) DATA CHUNK: $data\n";

	my ($current_blast, $current_prog, $current_vers, $current_db);
	my $prev_blast;
        my $contains_translation = 0;

### steve --- Wed Mar 15 02:48:07 2000
### In the process of addressing bug PR#95. Tricky.
### Using the $contains_translation to do so. Not complete
### and possibly won't fix. We'll see.

	# Check for header section. Start a new Blast object and
	# parse the description section.
#        if ($data =~ /\sQuery\s?=/s || ($contains_translation && $data =~ /Database:/s)) {
        if ($data =~ /\sQuery\s?=/s) {
	    $Blast->{'_blast_count'}++;
	    print STDERR ".", $Blast->{'_blast_count'} % 50 ? '' : "\n" if $MONITOR;

            if($data =~ /$Newline\s+Translating/so) {
                print STDERR "\nCONTAINS TRANSLATION\n";
                $contains_translation = 1;
            }

	    # If we're parsing a stream containing multiple reports,
	    # all subsequent header sections will contain the last hit of
	    # the previous report which needs to be parsed and added to that
	    # report if signifcant. It also contains the run parameters
	    # at the bottom of the Blast report.
#	    if($Blast->{'_blast_count'} > 1 || $contains_translation) {
            if($Blast->{'_blast_count'} > 1) {
#	      print STDERR "\nMULTI-BLAST STREAM.\n";
	      $Blast->{'_multi_stream'} = 1;

	      if($data =~ /(.+?)$Newline(<\w+>)?(T?BLAST[NPX])\s+(.+?)$Newline(.+)/so) {
		($current_prog, $current_vers, $data) = ($3, $4, $5);
		# Final chunk containing last hit and last footer.
		$Blast->{'_current_blast'}->_parse_alignment($1);
		$prev_blast = $Blast->{'_current_blast'}; # finalized.
#              }	elsif($contains_translation) {
#                  $data =~ /(T?BLAST[NPX])\s+(.+?)$Newline(.+)/so;
#                  ($current_prog, $current_vers, $data) = ($1, $2, $3);
	      }	else {
		$Blast->throw("Can't determine program type from BLAST report.",
			      "Checked for: @Blast_programs.");
		# This has important implications for how to handle interval
		# information for HSPs. TBLASTN uses nucleotides in query HSP
		# but amino acids in the sbjct HSP sequence.
	      }

	      if($data =~ m/Database:\s+(.+?)$Newline/so ) {
		$current_db = $1;
	      } else {
		# In some reports, the Database is only listed at end.
		#$Blast->warn("Can't determine database name from BLAST report.");
	      }

	      # Incyte_Fix:   Nasty Invisible Bug.
	      # Records in blast report are delimited by '>', but... when
	      #  there are no hits for a query, there won't be a '>'.  That
	      #  causes several blast reports to run together in the data
	      #  passed to this routine.  Need to get rid of non-hits in data
	      if ($data =~ /.+(No hits? found.+Sequences.+)/so) {
		  $data = $1;
	      }
	      # End Incyte_Fix

	    }

	    # Determine if we need to create a new Blast object
	    # or use the $self object for this method.

	    if($Blast->{'_multi_stream'} or $self->name eq 'Static Blast object') {
	      # Strict mode is not object-specific but may be someday.
#	      print STDERR "\nCreating new Blast object.\n";
	      $current_blast = new Bio::Tools::Blast(-STRICT => $strict);
	    } else {
	      $current_blast = $self;
	    }
	    $Blast->{'_current_blast'} = $current_blast;

	    # If we're not sharing stats, set data on current blast  object.
	    if(defined $current_prog and not $Blast->{'_share'}) {
	      $current_blast->program($current_prog);
	      $current_blast->program_version($current_vers);
	      $current_blast->database($current_db);
	    }

#	    print STDERR "CURRENT BLAST = ", $current_blast->name, "\n";
	    $current_blast->_parse_header($data);
	
	    # If there were any descriptions in the header,
	    # we know if there are any significant hits.
	    # No longer throwing exception if there were no significant hits
	    # and a -signif parameter was specified. Doing so prevents the
	    # construction of a Blast object, which could still be useful.
#	    if($current_blast->{'_has_descriptions'} and $Blast->{'_confirm_significance'} and not $current_blast->is_signif) {
#	      $current_blast->throw("No significant BLAST hits for ${\$current_blast->name}");

#	    }

	} # Done parsing header/description section

### For use with $contains_translation - not right - breaks regular report parsing.
# 	elsif(ref $Blast->{'_current_blast'} && $data !~ /\s*\w*\s*/s) {
 	elsif(ref $Blast->{'_current_blast'} ) {
	    # Process an alignment section.
	    $current_blast = $Blast->{'_current_blast'};
#	    print STDERR "\nCONTINUING PROCESSING ALN WITH ", $current_blast->name, "\n";
#	    print STDERR "DATA: $data\n";
            eval {
                $current_blast->_parse_alignment($data);
            };
            if($@) {
     #           push @{$self->{'_blast_errs'}}, $@;
            }
	}
	
	# If the current Blast object has been completely parsed
	# (occurs with a single Blast stream), or if there is a previous
	# Blast object (occurs with a multi Blast stream),
	# execute a supplied function on it or store it in a supplied array.

	if( defined $prev_blast or $current_blast->{'_found_params'}) {
	  my $finished_blast = defined($prev_blast) ? $prev_blast : $current_blast;
	
	  $finished_blast->_report_errors();
#	  print STDERR "\nNEW BLAST OBJECT: ${\$finished_blast->name}\n";

	  if($exec_func) {
#	    print STDERR "  RUNNING EXEC_FUNC...\n";
	    &$exec_func($finished_blast);  # ignoring any return value.
	    # Report processed, no longer need object.
	    $finished_blast->destroy;
	    undef $finished_blast;
	  } elsif($save_a) {
#	    print STDERR "  SAVING IN ARRAY...\n";
	    # We've already verified that if there is no exec_func
	    # then there must be a $save_array
	    push @$save_a, $finished_blast;
	  }
	}
	1;
      }
  }

=head2 _report_errors

 Title   : _report_errors
 Usage   : n/a; Internal method called by _get_parse_blast_func().
 Purpose : Throw or warn about any errors encountered.
 Returns : n/a
 Args    : n/a
 Throws  : If all hits generated exceptions, raise exception
         :   (a fatal event for the Blast object.)
         : If some hits were okay but some were bad, generate a warning
         :   (a few bad applies should not spoil the bunch).
         :   This usually indicates a limiting B-value.
         : When the parsing code fails, it is either all or nothing.

=cut

#-------------------
sub _report_errors {
#-------------------
  my $self = shift;

  return unless ref($self->{'_blast_errs'});
#  ref($self->{'_blast_errs'}) || (print STDERR "\nNO ERRORS\n", return );

  my @errs = @{$self->{'_blast_errs'}};

  if(scalar @errs) {
    my ($str);
    @{$self->{'_blast_errs'}} = (); # clear the errs on the object.
    # When there are many errors, in most of the cases, they are
    # caused by the same problem. Only need to see full data for
    # the first one.
    if(scalar @errs > 2) {
      $str = "SHOWING FIRST EXCEPTION ONLY:\n$errs[0]";
      $self->clear_err();  # clearing the existing set of errors (conserve memory).
      # Not necessary, unless the -RECORD_ERR =>1
      # constructor option was used for Blast object.
    } else {
      $str = join("\n",@errs);
    }

    if(not $self->{'_num_hits_significant'}) {
      $self->throw(sprintf("Failed to parse any hit data (n=%d).", scalar(@errs)),
		   "\n\nTRAPPED EXCEPTION(S):\n$str\nEND TRAPPED EXCEPTION(S)\n"
		  );
    } else {
      $self->warn(sprintf("Some potential hits were not parsed (n=%d).", scalar(@errs)),
		  @errs > 2 ? "This may be due to a limiting B value (max alignment listings)." : "",
		  "\n\nTRAPPED EXCEPTION(S):\n$str\nEND TRAPPED EXCEPTION(S)\n"
		 );
    }
  }
}

=head2 _parse_header

 Usage     : n/a; called automatically by the _get_parse_blast_func().
 Purpose   : Parses the header section of a BLAST report.
 Argument  : String containing the header+description section of a BLAST report.
 Throws    : Exception if description data cannot be parsed properly.
           : Exception if there is a 'FATAL' error in the Blast report.
           : Warning if there is a 'WARNING' in the Blast report.
           : Warning if there are no significant hits.
 Comments  : Description section contains a single line for each hit listing
           : the seq id, description, score, Expect or P-value, etc.

See Also   : L<_get_parse_blast_func()|_get_parse_blast_func>

=cut

#----------------------
sub _parse_header {
#----------------------
    my( $self, $data ) = @_;

#    print STDERR "\n$ID: PARSING HEADER\n"; #$data\n";

    $data =~ s/^\s+|\s+>?$//sg;

    if($data =~ /<HTML/i) {
      $self->throw("Can't parse HTML-formatted BLAST reports.",
#		    "Such reports can be parsed with a special parsing \n".
#		    "script included in the examples/blast directory \n".
#		    "of the Bioperl distribution. (TODO)"
		  );
      # This was the old strategy, can't do it with new strategy
      # since we don't have the whole report in one chunk.
      # This could be the basis for the "special parsing script".
#	 require Bio::Tools::Blast::HTML;
#	 Bio::Tools::Blast::HTML->import(&strip_html);
#	 &strip_html(\$data);
    }
	
    $data =~ /WARNING: (.+?)$Newline$Newline/so and $self->warn("$1") if $self->strict;
    $data =~ /FATAL: (.+?)$Newline$Newline/so and $self->throw("FATAL BLAST ERROR = $1");
    # No longer throwing exception when no hits were found. Still reporting it.
    $data =~ /No hits? found/i and $self->warn("No hits were found.") if $self->strict;

    # If this is the first Blast, the program, version, and database info
    # pertain to it. Otherwise, they are for the previous report and have
    # already been parsed out.
    # Data is stored in the static Blast object. Data for subsequent reports
    # will be stored in separate objects if the -share parameter is not set.
    # See _get_parse_blast_func().

    if($Blast->{'_blast_count'} == 1) {
      if($data =~ /(<\w+>)?(T?BLAST[NPX])\s+(.+?)$Newline/so) {
	$Blast->program($2);
	$Blast->program_version($3);
      }	else {
	$self->throw("Can't determine program type from BLAST report.",
		     "Checked for: @Blast_programs.");
	# This has important implications for how to handle interval
	# information for HSPs. TBLASTN uses nucleotides in query HSP
	# but amino acids in the sbjct HSP sequence.
      }

      if($data =~ m/Database:\s+(.+?)$Newline/so ) {
	$Blast->database($1);
      } else {
	# In some reports, the Database is only listed at end.
	#$self->warn("Can't determine database name from BLAST report (_parse_header)\n$data\n.");
      }
    }

    my ($header, $descriptions);

    ## For efficiency reasons, we want to to avoid using $' and $`.
    ## Therefore using single-line mode pattern matching.

    if($data =~ /(.+?)\nSequences producing.+?\n(.+)/s ) {
        ($header, $descriptions) = ($1, $2);
	$self->{'_has_descriptions'} = 1;
    } else {
        $header = $data;
	$self->{'_has_descriptions'} = 0;
	# Blast reports can legally lack description section. No need to warn.
	#push @{$self->{'_blast_errs'}}, "Can't parse description data.";
    }

    $self->_set_query($header);  # The name of the sequence will appear in error report.
#    print STDERR "\nQUERY = ", $Blast->{'_current_blast'}->query, "\n";

    $self->_set_date($header) if $Blast->{'_get_stats'};
    $self->_set_length($header);

#    not $Blast->{'_confirm_significance'} and print STDERR "\nNOT PARSING DESCRIPTIONS.\n";

    # Setting the absolute max and min significance levels.
    $self->{'_highestSignif'} = 0;
    $self->{'_lowestSignif'} = $DEFAULT_SIGNIF;

    if ($Blast->{'_confirm_significance'} || $Blast->{'_no_aligns'}) {
      $self->_parse_descriptions($descriptions) if $descriptions;
    } else {
      $self->{'_is_significant'} = 1;
    }
  }

#-----------------------
sub _parse_descriptions {
#-----------------------
  my ($self, $desc) = @_;

    # NOTE: This method will not be called if the report lacks
    #       a description section.

#    print STDERR "\nPARSING DESCRIPTION DATA\n";

    my @descriptions = split( $Newline, $desc);
    my($line);

    # NOW step through each line parsing out the P/Expect value
    # All we really need to do is check the first one, if it doesn't
    # meet the significance requirement, we can skip the report.
    # BUT: we want to collect data for all hits anyway to get min/max signif.

    my $my_signif = $self->signif;
    my $layout_set = $Blast->{'_layout'} || 0;
    my $layout;
    my $count = 0;
    my $sig;

    desc_loop:
  foreach $line (@descriptions) {
      $count++;
      last desc_loop if $line =~ / NONE |End of List/;
      next desc_loop if $line =~ /^\s*$/;
      next desc_loop if $line =~ /^\.\./;

	## Checking the significance value (P- or Expect value) of the hit
	## in the description line.

      # These regexps need testing on a variety of reports.
      if ( $line =~ /\d+\s{1,5}[\de.-]+\s*$/) {
	$layout = 2;
      } elsif( $line =~ /\d+\s{1,5}[\de.-]+\s{1,}\d+\s*$/) {
	$layout = 1;
      } else {
	$self->warn("Can't parse significance data in description line $line");
	next desc_loop;
      }
      not $layout_set and ($self->_layout($layout), $layout_set = 1);

      $sig = &_parse_signif( $line, $layout );

#      print STDERR "  Parsed signif ($layout) = $sig\n";

      last desc_loop if ($sig > $my_signif and not $Blast->{'_check_all'});
      $self->_process_significance($sig, $my_signif);
    }

#  printf "\n%d SIGNIFICANT HITS.\nDONE PARSING DESCRIPTIONS.\n", $self->{'_num_hits_significant'};
}

sub _process_significance {
    my($self, $sig, $my_signif) = @_;

    $self->{'_highestSignif'} = ($sig > $self->{'_highestSignif'})
   	                        ? $sig : $self->{'_highestSignif'};

    $self->{'_lowestSignif'} = ($sig < $self->{'_lowestSignif'})
                                 ? $sig : $self->{'_lowestSignif'};

    # Significance value assessment.
    $sig <= $my_signif and $self->{'_num_hits_significant'}++;
    $self->{'_num_hits'}++;

    $self->{'_is_significant'} = 1 if $self->{'_num_hits_significant'};
}

=head2 _parse_alignment

 Usage     : n/a; called automatically by the _get_parse_blast_func().
 Purpose   : Parses a single alignment section of a BLAST report.
 Argument  : String containing the alignment section.
 Throws    : n/a; All errors are trapped while parsing the hit data
           : and are processed as a group when the report is
           : completely processed (See _report_errors()).
           :
 Comments  : Alignment section contains all HSPs for a hit.
           : Requires Bio::Tools::Blast::Sbjct.pm.
           : Optionally calls a filter function to screen the hit on arbitrary
           : criteria. If the filter function returns true for a given hit,
           : that hit will be skipped.
           :
           : If the Blast object was created with -check_all_hits set to true,
           : all hits will be checked for significance and processed if necessary.
           : If this field is false, the parsing will stop after the first
           : non-significant hit.
           : See parse() for description of parsing parameters.

See Also   : L<parse()|parse>, L<_get_parse_blast_func()|_get_parse_blast_func>, L<_report_errors()|_report_errors>, B<Bio::Tools::Blast::Sbjct()>, L<Links to related modules>

=cut

#----------------------
sub _parse_alignment {
#----------------------
# This method always needs to check detect if the $data argument
# contains the footer of a Blast report, indicating the last chunk
# of a single Blast stream.

    my( $self, $data ) = @_;

#    printf STDERR "\nPARSING ALIGNMENT DATA for %s $self.\n", $self->name;

    # NOTE: $self->{'_current_hit'} is an instance variable
    #       The $Blast object will not have this member.

    # If all of the significant hits have been parsed,
    # return if we're not checking all or if we don't need to get
    # the Blast stats (parameters at footer of report).
    if(defined $self->{'_current_hit'} and
      defined $self->{'_num_hits_significant'}) {
      return if $self->{'_current_hit'} >= $self->{'_num_hits_significant'} and
	not ($Blast->{'_check_all'} or $Blast->{'_get_stats'});
    }

    # Check for the presence of the Blast footer section.
    # _parse_footer returns the alignment section.
    $data = $self->_parse_footer($data);

    # Return if we're only interested in the best hit.
    # This has to occur after checking for the parameters section
    # in the footer (since we may still be interested in them).
    return if $Blast->best and ( defined $self->{'_current_hit'} and $self->{'_current_hit'} >=1);

#    print "RETURNED FROM _parse_footer (", $self->to_string, ")";
#    print "\n  --> FOUND PARAMS.\n" if $self->{'_found_params'};
#    print "\n  --> DID NOT FIND PARAMS.\n" unless $self->{'_found_params'};

    require Bio::Tools::Blast::Sbjct;

    $data =~ s/^\s+|\s+>?$//sg;
    $data =~ s/$Newline$Newline/$Newline/sog;  # remove blank lines.
    my @data = split($Newline, $data);
    push @data, 'end';

#    print STDERR "\nALIGNMENT DATA:\n$data\n";

    my $prog       = $self->program;
    my $check_all  = $Blast->{'_check_all'};
    my $filt_func  = $Blast->{'_filt_func'} || 0;
    my $signif_fmt = $Blast->{'_signif_fmt'};
    my $my_signif  = $self->signif;
    my $err;

    # Now construct the Sbjct objects from the alignment section

#	debug(1);

    $self->{'_current_hit'}++;

    # If not confirming significance, _parse_descriptions will not have been run,
    # so we need to count the total number of hits here.
    if( not $Blast->{'_confirm_significance'}) {
      $self->{'_num_hits'}++;
    }

    if($Blast->{'_no_aligns'}) {
#        printf STDERR "\nNOT PARSING ALIGNMENT DATA\n";
        return;
    }

    my $hit;  # Must be my'ed within hit_loop.
    eval {
      $hit = new Bio::Tools::Blast::Sbjct (-DATA      =>\@data,
					   -PARENT    =>$self,
					   -NAME      =>$self->{'_current_hit'},
					   -RANK      =>$self->{'_current_hit'},
					   -RANK_BY   =>'order',
					   -PROGRAM   =>$prog,
					   -SIGNIF_FMT=>$signif_fmt,
					   -OVERLAP   =>$Blast->{'_overlap'} || $MAX_HSP_OVERLAP,
					  );
#      printf STDERR "NEW HIT: %s, SIGNIFICANCE = %g\n", $hit->name, $hit->expect;  <STDIN>;
      # The BLAST report may have not had a description section.
      if(not $self->{'_has_descriptions'}) {
	  $self->_process_significance($hit->signif, $my_signif);
      }
    };

    if($@) {
      # Throwing lots of errors can slow down the code substantially.
      # Error handling code is not that efficient.
      #print STDERR "\nERROR _parse_alignment: $@\n";
      push @{$self->{'_blast_errs'}}, $@;
      $hit->destroy if ref $hit;
      undef $hit;
    } else {
      # Collect overall signif data if we don't already have it,
      # (as occurs if no -signif parameter is supplied).
      my $hit_signif = $hit->signif;

      if (not $Blast->{'_confirm_significance'} ) {
	$self->{'_highestSignif'} = ($hit_signif > $self->{'_highestSignif'})
                                    ? $hit_signif : $self->{'_highestSignif'};

	$self->{'_lowestSignif'} = ($hit_signif < $self->{'_lowestSignif'})
                                    ? $hit_signif : $self->{'_lowestSignif'};
      }

      # Test significance using custom function (if supplied)
      if($filt_func) {
	if(&$filt_func($hit)) {
	  push @{$self->{'_hits'}}, $hit;
	} else {
	  $hit->destroy; undef $hit;
	}
      } elsif($hit_signif <= $my_signif) {
	push @{$self->{'_hits'}}, $hit;
      }
    }

  }

=head2 _parse_footer

 Usage     : n/a; internal function. called by _parse_alignment()
 Purpose   : Extracts statistical and other parameters from the BLAST report.
           : Sets various key elements such as the program and version,
           : gapping, and the layout for the report (blast1 or blast2).
 Argument  : Data to be parsed.
 Returns   : String containing an alignment section for processing by
           : _parse_alignment().
 Throws    : Exception if cannot find the parameters section of report.
           : Warning if cannot determine if gapping was used.
           : Warning if cannot determine the scoring matrix used.
 Comments  : This method must always get called, even if the -STATS
           : parse() parameter is false. The reason is that the layout
           : of the report  and the presence of gapping must always be set.
           : The determination whether to set additional stats is made
           : by methods called by _parse_footer().

See Also   : L<parse()|parse>, L<_parse_alignment()|_parse_alignment>, L<_set_database()|_set_database>

=cut

#---------------------
sub _parse_footer {
#---------------------
# Basic strategy:
# 1. figure out if we're supposed to get the stats,
# 2. figure out if the stats are to be shared. some, not all can be shared
#    (eg., db info and matrix can be shared, karlin altschul params cannot.
#    However, this method assumes they are all sharable.)
# 3. Parse the stats.
# 4. return the block before the parameters section if the supplied data
#    contains a footer parameters section.

    my ($self, $data) = @_;
    my ($client, $last_align, $params);

#    printf STDERR "\nPARSING PARAMETERS for %s $self.\n", $self->name;

    # Should the parameters be shared?
    # If so, set $self to be the static $Blast object and return if
    # the parameters were already set.
    # Before returning, we need to extract the last alignment section
    # from the parameter section, if any.

    if ($Blast->{'_share'}) {
      $client = $self;
      $self = $Blast if $Blast->{'_share'};
    }

    my $get_stats = $Blast->{'_get_stats'};

    if( $data =~ /(.+?)${Newline}CPU time: (.*)/so) {
	# NCBI-Blast2 format (v2.04).
	($last_align, $params) = ($1, $2);
	return $last_align if $client->{'_found_params'};
	$self->_set_blast2_stats($params);

    } elsif( $data =~ /(.+?)${Newline}Parameters:(.*)/so) {
    # NCBI-Blast1 or WashU-Blast2 format.
    ($last_align, $params) = ($1, $2);
    return $last_align if $client->{'_found_params'};
    $self->_set_blast1_stats($params);

    } elsif( $data =~ /(.+?)$Newline\s+Database:(.*)/so) {
        # Gotta watch out for confusion with the Database: line in the header
        # which will be present in the last hit of an internal Blast report
        # in a multi-report stream.

	# NCBI-Blast2 format (v2.05).
	($last_align, $params) = ($1, $2);
	return $last_align if $client->{'_found_params'};
        $self->_set_blast2_stats($params);

    } elsif( $data =~ /(.+?)$Newline\s*Searching/so) {
	# trying to detect a Searching at the end of a PSI-blast round.
        # Gotta watch out for confusion with the Searching line in the header
        # which will be present in the last hit of an internal Blast report
        # in a multi-report, non-PSI-blast stream.

	# PSI-Blast format (v2.08).
	($last_align) = ($1);
	return $last_align; # if $client->{'_found_params'};
    }

    # If parameter section was found, set a boolean,
    # otherwise return original data.

    if (defined($params)) {
      $client->{'_found_params'} = 1;
    } else {
      return $data;
    }

    $self->_set_database($params) if $get_stats;

    # The {'_gapped'} member should be set in the _set_blast?_stats() call.
    # This is a last minute attempt to deduce it.

    if(!defined($self->{'_gapped'})) {
	if($self->program_version() =~ /^1/) {
	    $self->{'_gapped'} = 0;
	} else {
	    if($self->strict > 0) {
		$self->warn("Can't determine if gapping was used. Assuming gapped.");
	    }
	    $self->{'_gapped'} = 1;
	}
    }

    return $last_align;
}

=head2 _set_blast2_stats

 Usage     : n/a; internal function called by _parse_footer()
 Purpose   : Extracts statistical and other parameters from BLAST2 report footer.
           : Stats collected: database release, gapping,
           : posted date, matrix used, filter used, Karlin-Altschul parameters,
           : E, S, T, X, W.
 Throws    : Exception if cannot get "Parameters" section of Blast report.

See Also   : L<parse()|parse>, L<_parse_footer()|_parse_footer>, L<_set_database()|_set_database>, B<Bio::Tools::SeqAnal::set_date()>,L<Links to related modules>

=cut

#---------------------'
sub _set_blast2_stats {
#---------------------
    my ($self, $data) = (@_);

    if($data =~ /$Newline\s*Gapped/so) {
	$self->{'_gapped'} = 1;
    } else {
 	$self->{'_gapped'} = 0;
    }

    # Other stats are not always essential.
    return unless $Blast->{'_get_stats'};

    # Blast2 Doesn't report what filter was used in the parameters section.
    # It just gives a warning that *some* filter was used in the header.
    # You just have to know the defaults (currently: protein = SEG, nucl = DUST).
    if($data =~ /\bfiltered\b/si) {
	$self->{'_filter'} = 'DEFAULT FILTER';
    } else {
	$self->{'_filter'} = 'NONE';
    }

    if($data =~ /Gapped$Newline\s*Lambda +K +H$Newline +(.+?)$Newline/so) {
	my ($l, $k, $h) = split(/\s+/, $1);
	$self->{'_lambda'} = $l || 'UNKNOWN';
	$self->{'_k'} = $k || 'UNKNOWN';
	$self->{'_h'} = $h || 'UNKNOWN';
    } elsif($data =~ /Lambda +K +H$Newline +(.+?)$Newline/so) {
	my ($l, $k, $h) = split(/\s+/, $1);
	$self->{'_lambda'} = $l || 'UNKNOWN';
	$self->{'_k'} = $k || 'UNKNOWN';
	$self->{'_h'} = $h || 'UNKNOWN';
    }

    if($data =~ /$Newline\s*Matrix: (.+?)$Newline/so) {
	$self->{'_matrix'} = $1;
    } else {
	$self->{'_matrix'} = $DEFAULT_MATRIX.'?';
	if($self->strict > 0) {
	    $self->warn("Can't determine scoring matrix. Assuming $DEFAULT_MATRIX.");
	}
    }

    if($data =~ /$Newline\s*Gap Penalties: Existence: +(\d+), +Extension: (\d+)$Newline/so) {
	$self->{'_gapCreation'} = $1;
	$self->{'_gapExtension'} = $2;
    }
    if($data =~ /sequences better than (\d+):/s) {
	$self->{'_expect'} = $1;
    }

    if($data =~ /$Newline\s*T: (\d+)/o) { $self->{'_word_size'} = $1; }
    if($data =~ /$Newline\s*A: (\d+)/o) { $self->{'_a'} = $1; }
    if($data =~ /$Newline\s*S1: (\d+)/o) { $self->{'_s'} = $1; }
    if($data =~ /$Newline\s*S2: (\d+)/o) { $self->{'_s'} .= ", $1"; }
    if($data =~ /$Newline\s*X1: (\d+)/o) { $self->{'_x1'} = $1; }
    if($data =~ /$Newline\s*X2: (\d+)/o) { $self->{'_x2'} = $1; }
}

=head2 _set_blast1_stats

 Usage     : n/a; internal function called by _parse_footer()
 Purpose   : Extracts statistical and other parameters from BLAST 1.x style eports.
           : Handles NCBI Blast1 and WashU-Blast2 formats.
           : Stats collected: database release, gapping,
           : posted date, matrix used, filter used, Karlin-Altschul parameters,
           : E, S, T, X, W.

See Also   : L<parse()|parse>, L<_parse_footer()|_parse_footer>, L<_set_database()|_set_database>, B<Bio::Tools::SeqAnal::set_date()>,L<Links to related modules>

=cut

#----------------------
sub _set_blast1_stats {
#----------------------
    my ($self, $data) = (@_);

    if(!$self->{'_gapped'} and $self->program_version() =~ /^2[\w\-\.]+WashU/) {
	$self->_set_gapping_wu($data);
    } else {
	$self->{'_gapped'} = 0;
    }

    # Other stats are not always essential.
    return unless $Blast->{'_get_stats'};

    if($data =~ /filter=(.+?)$Newline/so) {
	$self->{'_filter'} = $1;
    } elsif($data =~ /filter$Newline +(.+?)$Newline/so) {
	$self->{'_filter'} = $1;
    } else {
	$self->{'_filter'} = 'NONE';
    }

    if($data =~ /$Newline\s*E=(\d+)$Newline/so) {  $self->{'_expect'} = $1; }

    if($data =~ /$Newline\s*M=(\w+)$Newline/so) {  $self->{'_matrix'} = $1; }

    if($data =~ /\s*Frame  MatID Matrix name .+?$Newline +(.+?)$Newline/so) {
	## WU-Blast2.
	my ($fr, $mid, $mat, $lu, $ku, $hu, $lc, $kc, $hc) = split(/\s+/,$1);
	$self->{'_matrix'} = $mat || 'UNKNOWN';
	$self->{'_lambda'} = $lu || 'UNKNOWN';
	$self->{'_k'} = $ku || 'UNKNOWN';
	$self->{'_h'} = $hu || 'UNKNOWN';
	
    } elsif($data =~ /Lambda +K +H$Newline +(.+?)$Newline/so) {
	## NCBI-Blast1.
	my ($l, $k, $h) = split(/\s+/, $1);
	$self->{'_lambda'} = $l || 'UNKNOWN';
	$self->{'_k'} = $k || 'UNKNOWN';
	$self->{'_h'} = $h || 'UNKNOWN';
    }

    if($data =~ /E +S +W +T +X.+?$Newline +(.+?)$Newline/so) {
	# WashU-Blast2
	my ($fr, $mid, $len, $elen, $e, $s, $w, $t, $x, $e2, $s2) = split(/\s+/,$1);
	$self->{'_expect'} ||= $e || 'UNKNOWN';
	$self->{'_s'} = $s || 'UNKNOWN';
	$self->{'_word_size'} = $w || 'UNKNOWN';
	$self->{'_t'} = $t || 'UNKNOWN';
	$self->{'_x'} = $x || 'UNKNOWN';

    } elsif($data =~ /E +S +T1 +T2 +X1 +X2 +W +Gap$Newline +(.+?)$Newline/so) {
	## NCBI-Blast1.
	my ($e, $s, $t1, $t2, $x1, $x2, $w, $gap) = split(/\s+/,$1);
	$self->{'_expect'} ||= $e || 'UNKNOWN';
	$self->{'_s'} = $s || 'UNKNOWN';
	$self->{'_word_size'} = $w || 'UNKNOWN';
	$self->{'_t1'} = $t1 || 'UNKNOWN';
	$self->{'_t2'} = $t2 || 'UNKNOWN';
	$self->{'_x1'} = $x1 || 'UNKNOWN';
	$self->{'_x2'} = $x2 || 'UNKNOWN';
	$self->{'_gap'} = $gap || 'UNKNOWN';
    }

    if(!$self->{'_matrix'}) {
	$self->{'_matrix'} = $DEFAULT_MATRIX.'?';
	if($self->strict > 0) {
	    $self->warn("Can't determine scoring matrix. Assuming $DEFAULT_MATRIX.");
	}
    }
}

=head2 _set_gapping_wu

 Usage     : n/a; internal function called by _set_blast1_stats()
 Purpose   : Determine if gapping_wu was on for WashU Blast reports.
 Comments  : In earlier versions, gapping was always specified
           : but in the current version (2.0a19MP), gapping is on by default
           : and there is no positive "gapping" indicator in the Parameters
           : section.

See Also   : L<_set_blast1_stats()|_set_blast1_stats>

=cut

#--------------------
sub _set_gapping_wu {
#--------------------
    my ($self, $data) = @_;

    if($data =~ /gaps?$Newline/so) {
	$self->{'_gapped'} = ($data =~ /nogaps?$Newline/so) ? 0 : 1;
    } else {
	$self->{'_gapped'} = 1;
    }
}

=head2 _set_date

 Usage     : n/a; internal function called by _parse_footer()
 Purpose   : Determine the date on which the Blast analysis was performed.
 Comments  : Date information is not consistently added to Blast output.
           : Uses superclass method set_date() to set date from the file,
           : (if any).

See Also   : L<_parse_footer()|_parse_footer>, B<Bio::Tools::SeqAnal::set_date()>,L<Links to related modules>

=cut

#--------------
sub _set_date {
#--------------
    my $self = shift;
    my $data = shift;

    ### Network BLAST reports from NCBI are time stamped as follows:
    #Fri Apr 18 15:55:41 EDT 1997, Up 1 day, 19 mins, 1 user, load: 19.54, 19.13, 17.77
    if($data =~ /Start:\s+(.+?)\s+End:/s) {
	## Calling superclass method to set the date.
	## If we can't get date from the report, file date is obtained.
	$self->set_date($1);
    } elsif($data =~ /Date:\s+(.*?)$Newline/so) {
	## E-mailed reports have a Date: field
	$self->set_date($1);
    } elsif( $data =~ /done\s+at (.+?)$Newline/so ) {
	$self->set_date($1);
    } elsif( $data =~ /$Newline([\w:, ]+), Up \d+/so ) {
	$self->set_date($1);
    } else {
	## Otherwise, let superclass attempt to get the file creation date.
	$self->set_date() if $self->file;
    }
}

=head2 _set_length

 Usage     : n/a; called automatically during Blast report parsing.
 Purpose   : Sets the length of the query sequence (extracted from report).
 Returns   : integer (length of the query sequence)
 Throws    : Exception if cannot determine the query sequence length from
           :           the BLAST report.
           : Exception if the length is below the min_length cutoff (if any).
 Comments  : The logic here is a bit different from the other _set_XXX()
           : methods since the significance of the BLAST report is assessed
           : if MIN_LENGTH is set.

See Also   : B<Bio::Tools::SeqAnal::length()>, L<Links to related modules>

=cut

#---------------
sub _set_length {
#---------------
    my ($self, $data) = @_;

    my ($length);
    if( $data =~ m/$Newline\s+\(([\d|,]+) letters[\);]/so ) {
	$length = $1;
	$length =~ s/,//g;
#	printf "Length = $length in BLAST for %s$Newline",$self->name; <STDIN>;
    } else {
	$self->throw("Can't determine sequence length from BLAST report.");
    }

    my($sig_len);
    if(defined($Blast->{'_min_length'})) {
      local $^W = 0;
      if($length < $Blast->{'_min_len'}) {
	$self->throw("Query sequence too short for ${\$self->name} ($length)",
		     "Minimum  length is $Blast->{'_min_len'}");
      }
    }

    $self->length($length);  # defined in superclass.
}

=head2 _set_database

 Usage     : n/a; called automatically during Blast report parsing.
 Purpose   : Sets the name of the database used by the BLAST analysis.
           : Extracted from raw BLAST report.
 Throws    : Exception if the name of the database cannot be determined.
 Comments  : The database name is used by methods or related objects
           : for database-specific parsing.

See Also   : L<parse()|parse>, B<Bio::Tools::SeqAnal::database()>,B<Bio::Tools::SeqAnal::_set_db_stats()>,L<Links to related modules>

=cut

#------------------
sub _set_database {
#------------------
# This now only sets data base information extracted from the report footer.

    my ($self, $data) = @_;

    my ($name, $date, $lets, $seqs);

    my $strict = $self->strict > 0;

    # This is fail-safe since DB name usually gets set in _parse_header()
    # In some reports, the database is only listed at bottom (NCBI 2.0.8).
    if($data =~ m/Database: +(.+?)$Newline/so ) {
      $name = $1;
    } elsif(not $self->database) {
      $self->warn("Can't determine database name from BLAST report.");
    }

    if($data =~ m/Posted date: +(.+?)$Newline/so ) {
	$date = $1;
    } elsif($data =~ m/Release date: +(.+?)$Newline/so ) {
	$date = $1;
    } elsif($strict) {
	$self->warn("Can't determine database release date.");
    }

    if($data =~ m/letters in database: +([\d,]+)/si ||
       $data =~ m/length of database: +([\d,]+)/si ) {
	$lets = $1;
    } elsif($strict) {
	$self->warn("Can't determine number of letters in database.\n$data\n");
    }

    if($data =~ m/sequences in database: +([\d,]+)/si ||
      $data =~ m/number of sequences: +([\d,]+)/si ) {
	$seqs = $1;
    } elsif($strict) {
	$self->warn("Can't determine number of sequences in database.\n$data\n");
    }

    $self->_set_db_stats( -NAME    => $name,
			  -RELEASE => $date || '',
			  -LETTERS => $lets || '',
			  -SEQS    => $seqs || ''
			  );
}

=head2 _set_query

 Usage     : n/a; called automatically during Blast report parsing.
 Purpose   : Set the name of the query and the query description.
           : Extracted from the raw BLAST report.
 Returns   : String containing name of query extracted from report.
 Throws    : Warning if the name of the query cannont be obtained.

See Also   : B<Bio::Tools::SeqAnal::query_desc()>,L<Links to related modules>

=cut

#---------------
sub _set_query {
#---------------
    my $self = shift;
    my $data = shift;

    if($data =~ m/${Newline}Query= *(.+?)$Newline/so ) {
	my $info = $1;
	$info =~ s/TITLE //;
	# Split the query line into two parts.
	# Using \s instead of ' '
	$info =~ /(\S+?)\s(.*)/;
	$self->query_desc($2 || '');
	# set name of Blast object and return.
	$self->name($1 || 'UNKNOWN');
    } else {
	$self->warn("Can't determine query sequence name from BLAST report.");
    }
#    print STDERR "$Newline  NAME = ${\$self->name}$Newline";
}

=head2 _parse_signif

 Usage     : &_parse_signif(string, layout, gapped);
           : This is a class function.
 Purpose   : Extracts the P- or Expect value from a single line of a BLAST description section.
 Example   : &_parse_signif("PDB_UNIQUEP:3HSC_  heat-shock cognate ...   799  4.0e-206  2", 1);
           : &_parse_signif("gi|758803  (U23828) peritrophin-95 precurs   38  0.19", 2);
 Argument  : string = line from BLAST description section
           : layout = integer (1 or 2)
           : gapped = boolean (true if gapped Blast).
 Returns   : Float (0.001 or 1e-03)
 Status    : Static

=cut

#------------------
sub _parse_signif {
#------------------
    my ($line, $layout, $gapped) = @_;

    local $_ = $line;
    my @linedat = split();

    # When processing both Blast1 and Blast2 reports
    # in the same run, offset needs to be configured each time.

    my $offset  = 0;
    $offset  = 1 if $layout == 1 or not $gapped;

    my $signif = $linedat[ $#linedat - $offset ];

    # fail-safe check
    if(not $signif =~ /[.-]/) {
	$offset = ($offset == 0 ? 1 : 0);
	$signif = $linedat[ $#linedat - $offset ];
    }

    $signif = "1$signif" if $signif =~ /^e/i;
    return $signif;
}

##
## BEGIN ACCESSOR METHODS THAT INCORPORATE THE STATIC $Blast OBJECT.
##

sub program {
## Overridden method to incorporate the BLAST object.
    my $self = shift;
    return $self->SUPER::program(@_) if @_;          # set
    $self->SUPER::program || $Blast->SUPER::program; # get
}

sub program_version {
## Overridden method to incorporate the BLAST object.
    my $self = shift;
    return $self->SUPER::program_version(@_) if @_;                  # set
    $self->SUPER::program_version || $Blast->SUPER::program_version; # get
}

sub database {
## Overridden method to incorporate the BLAST object.
    my $self = shift;
    return $self->SUPER::database(@_) if @_;           # set
    $self->SUPER::database || $Blast->SUPER::database; # get
}

sub database_letters {
## Overridden method to incorporate the BLAST object.
    my $self = shift;
    return $self->SUPER::database_letters(@_) if @_;           # set
    $self->SUPER::database_letters || $Blast->SUPER::database_letters; # get
}

sub database_release {
## Overridden method to incorporate the BLAST object.
    my $self = shift;
    return $self->SUPER::database_release(@_) if @_;           # set
    $self->SUPER::database_release || $Blast->SUPER::database_release; # get
}

sub database_seqs {
## Overridden method to incorporate the BLAST object.
    my $self = shift;
    return $self->SUPER::database_seqs(@_) if @_;           # set
    $self->SUPER::database_seqs || $Blast->SUPER::database_seqs; # get
}

sub date {
## Overridden method to incorporate the BLAST object.
    my $self = shift;
    return $self->SUPER::date(@_) if @_;       # set
    $self->SUPER::date || $Blast->SUPER::date; # get
}

sub best {
## Overridden method to incorporate the BLAST object.
    my $self = shift;
    return $Blast->SUPER::best(@_) if @_;       # set
    $Blast->SUPER::best; # get
}

=head2 signif

 Usage     : $blast->signif();
 Purpose   : Gets the P or Expect value used as significance screening cutoff.
 Returns   : Scientific notation number with this format: 1.0e-05.
 Argument  : n/a
 Comments  : Screening of significant hits uses the data provided on the
           : description line. For Blast1 and WU-Blast2, this data is P-value.
           : for Blast2 it is an Expect value.
           :
           : Obtains info from the static $Blast object if it has not been set
           : for the current object.

See Also   : L<_set_signif()|_set_signif>

=cut

#-----------
sub signif {
#-----------
    my $self = shift;
    my $sig = $self->{'_significance'} || $Blast->{'_significance'};
    sprintf "%.1e", $sig;
}

=head2 is_signif

 Usage     : $blast->is_signif();
 Purpose   : Determine if the BLAST report contains significant hits.
 Returns   : Boolean
 Argument  : n/a
 Comments  : BLAST reports without significant hits but with defined
           : significance criteria will throw exceptions during construction.
           : This obviates the need to check significant() for
           : such objects.

See Also   : L<_set_signif()|_set_signif>

=cut

#------------
sub is_signif { my $self = shift; return $self->{'_is_significant'}; }
#------------

# is_signif() doesn't incorporate the static $Blast object but is included
# here to be with the other 'signif' methods.

=head2 signif_fmt

 Usage     : $blast->signif_fmt( [FMT] );
 Purpose   : Allows retrieval of the P/Expect exponent values only
           : or as a two-element list (mantissa, exponent).
 Usage     : $blast_obj->signif_fmt('exp');
           : $blast_obj->signif_fmt('parts');
 Returns   : String or '' if not set.
 Argument  : String, FMT = 'exp' (return the exponent only)
           :             = 'parts'(return exponent + mantissa in 2-elem list)
           :              = undefined (return the raw value)
 Comments  : P/Expect values are still stored internally as the full,
           : scientific notation value.
           : This method uses the static $Blast object since this issue
           : will pertain to all Blast reports within a given set.
           : This setting is propagated to Bio::Tools::Blast::Sbjct.pm.

=cut

#-------------
sub signif_fmt {
#-------------
    my $self = shift;
    if(@_) { $Blast->{'_signif_fmt'} = shift; }
    $Blast->{'_signif_fmt'} || '';
}

=head2 min_length

 Usage     : $blast->min_length();
 Purpose   : Gets the query sequence length used as significance screening criteria.
 Returns   : Integer
 Argument  : n/a
 Comments  : Obtains info from the static $Blast object if it has not been set
           : for the current object.

See Also   : L<_set_signif()|_set_signif>, L<signif()|signif>

=cut

#--------------
sub min_length {
#--------------
    my $self = shift;
    $self->{'_min_length'} || $Blast->{'_min_length'};
}

=head2 gapped

 Usage     : $blast->gapped();
 Purpose   : Set/Get boolean indicator for gapped BLAST.
 Returns   : Boolean
 Argument  : n/a
 Comments  : Obtains info from the static $Blast object if it has not been set
           : for the current object.

=cut

#-----------
sub gapped {
#-----------
    my $self = shift;
    if(@_) { $self->{'_gapped'} = shift; }
    $self->{'_gapped'} || $Blast->{'_gapped'};
}

=head2 _get_stats

 Usage     : n/a; internal method.
 Purpose   : Set/Get indicator for collecting full statistics from report.
 Returns   : Boolean (0 | 1)
 Comments  : Obtains info from the static $Blast object which gets set
           : by _init_parse_params().

=cut

#---------------
sub _get_stats {
#---------------
    my $self = shift;
    $Blast->{'_get_stats'};
}

=head2 _layout

 Usage     : n/a; internal method.
 Purpose   : Set/Get indicator for the layout of the report.
 Returns   : Integer (1 | 2)
           : Defaults to 2 if not set.
 Comments  : Blast1 and WashU-Blast2 have a layout = 1.
           : This is intended for internal use by this and closely
           : allied modules like Sbjct.pm and HSP.pm.
           :
           : Obtains info from the static $Blast object if it has not been set
           : for the current object.

=cut

#------------
sub _layout {
#------------
    my $self = shift;
    if(@_) {
      # Optimization if we know all reports share the same stats.
      if($Blast->{'_share'}) {
	$Blast->{'_layout'} = shift;
      } else {
	$self->{'_layout'} = shift;
      }
    }
    $self->{'_layout'} || $Blast->{'_layout'} || 2;
}

##
## END ACCESSOR METHODS THAT INCORPORATE THE STATIC $Blast OBJECT.
##

=head2 hits

 Usage     : $blast->hits();
 Purpose   : Get a list containing all BLAST hit (Sbjct) objects.
           : Get the numbers of significant hits.
 Examples  : @hits       = $blast->hits();
           : $num_signif = $blast->hits();
 Returns   : List context : list of Bio::Tools::Blast::Sbjct.pm objects
           :                or an empty list if there are no hits.
           : Scalar context: integer (number of significant hits)
           :                 or zero if there are no hits.
           :                 (Equivalent to num_hits()).
 Argument  : n/a. Relies on wantarray.
 Throws    : n/a.
           : Not throwing exception because the absence of hits may have
           : resulted from stringent significance criteria, not a failure
           : set the hits.

See Also   : L<hit()|hit>, L<num_hits()|num_hits>, L<is_signif()|is_signif>, L<_set_signif()|_set_signif>

=cut

#----------
sub hits {
#----------
    my $self = shift;

    if(wantarray) {
        my @ary = ref($self->{'_hits'}) ? @{$self->{'_hits'}} : ();
        return @ary;
    } else {
        return $self->num_hits();
    }

#    my $num = ref($self->{'_hits'}) ? scalar(@{$self->{'_hits'}}) : 0;
#    my @ary = ref($self->{'_hits'}) ? @{$self->{'_hits'}} : ();
#
#    return wantarray
#        #  returning list containing all hits or empty list.
#        ?  $self->{'_is_significant'} ? @ary : ()
#        #  returning number of hits or 0.
#        :  $self->{'_is_significant'} ? $num : 0;
}

=head2 hit

 Example   : $blast_obj->hit( [class] )
 Purpose   : Get a specific hit object.
           : Provides some syntactic sugar for the hits() method.
 Usage     : $hitObj = $blast->hit();
           : $hitObj = $blast->hit('best');
           : $hitObj = $blast->hit('worst');
           : $hitObj = $blast->hit( $name );
 Returns   : Object reference for a Bio::Tools::Blast::Sbjct.pm object.
           : undef if there are no hit (Sbjct) objects defined.
 Argument  : Class (or no argument).
           :   No argument (default) = highest scoring hit (same as 'best').
           :   'best' or 'first' = highest scoring hit.
           :   'worst' or 'last' = lowest scoring hit.
           :   $name = retrieve a hit by seq id (case-insensitive).
 Throws    : Exception if the Blast object has no significant hits.
           : Exception if a hit cannot be found when supplying a specific
           : hit sequence identifier as an argument.
 Comments  : 'best'  = lowest significance value (P or Expect) among significant hits.
           : 'worst' = highest sigificance value (P or Expect) among significant hits.

See Also   : L<hits()|hits>, L<num_hits()|num_hits>, L<is_signif()|is_signif>

=cut

#---------
sub hit {
#---------
    my( $self, $option) = @_;
    $option ||= 'best';

    if($Blast->{'_no_aligns'} || ! ref($self->{'_hits'})) {
        return undef;
    }

    $self->{'_is_significant'} or
	$self->throw("There were no significant hits.",
		     "Use num_hits(), hits(), is_signif() to check.");

    my @hits = @{$self->{'_hits'}};

    return $hits[0]      if $option =~ /^(best|first|1)$/i;
    return $hits[$#hits] if $option =~ /^(worst|last)$/i;

    # Get hit by name.	
    foreach ( @hits ) {
	return $_ if $_->name() =~ /$option/i;
    }

    $self->throw("Can't get hit for: $option");
}

=head2 num_hits

 Usage     : $blast->num_hits( ['total'] );
 Purpose   : Get number of significant hits or number of total hits.
 Examples  : $num_signif = $blast-num_hits;
           : $num_total  = $blast->num_hits('total');
 Returns   : Integer
 Argument  : String = 'total' (or no argument).
           :   No argument (Default) = return number of significant hits.
           :   'total' = number of total hits.
 Throws    : n/a.
           : Not throwing exception because the absence of hits may have
           : resulted from stringent significance criteria, not a failure
           : set the hits.
 Comments  : A significant hit is defined as a hit with an expect value
           : (or P value for WU-Blast) at or below the -signif parameter
           : used when parsing the report. Additionally, if a filter function
           : was supplied, the significant hit must also pass that
           : criteria.

See Also   : L<hits()|hits>, L<hit()|hit>, L<is_signif()|is_signif>, L<_set_signif()|_set_signif>, L<parse()|parse>

=cut

#-------------
sub num_hits {
#-------------
    my( $self, $option) = @_;
    $option ||= '';

   $option =~ /total/i and return $self->{'_num_hits'} || 0;

    # Default: returning number of significant hits.
#    return $self->{'_num_hits_significant'} || 0;
#    return 0 if not ref $self->{'_hits'};

    if(ref $self->{'_hits'}) {
        return scalar(@{$self->{'_hits'}});
    } else {
        return $self->{'_num_hits_significant'} || 0;
    }
}

=head2 lowest_p

 Usage     : $blast->lowest_p()
 Purpose   : Get the lowest P-value among all hits in a BLAST report.
           : Syntactic sugar for $blast->hit('best')->p().
 Returns   : Float or scientific notation number.
           : Returns -1.0 if lowest_p has not been set.
 Argument  : n/a.
 Throws    : Exception if the Blast report does not report P-values
           : (as is the case for NCBI Blast2).
 Comments  : A value is returned regardless of whether or not there were
           : significant hits ($DEFAULT_SIGNIF, currently  999).

See Also   : L<lowest_expect()|lowest_expect>, L<lowest_signif()|lowest_signif>, L<highest_p()|highest_p>, L<signif_fmt()|signif_fmt>

=cut

#------------
sub lowest_p {
#------------
    my $self = shift;

    # Layout 2 = NCBI Blast 2.x does not report P-values.
    $self->_layout == 2 and
	$self->throw("Can't get P-value with BLAST2.",
		     "Use lowest_signif() or lowest_expect()");

    return $self->{'_lowestSignif'} || -1.0;
}

=head2 lowest_expect

 Usage     : $blast->lowest_expect()
 Purpose   : Get the lowest Expect value among all hits in a BLAST report.
           : Syntactic sugar for $blast->hit('best')->expect()
 Returns   : Float or scientific notation number.
           : Returns -1.0 if lowest_expect has not been set.
 Argument  : n/a.
 Throws    : Exception if there were no significant hits and the report
           : does not have Expect values on the description lines
           : (i.e., Blast1, WashU-Blast2).

See Also   : L<lowest_p()|lowest_p>, L<lowest_signif()|lowest_signif>, L<highest_expect()|highest_expect>, L<signif_fmt()|signif_fmt>

=cut

#------------------
sub lowest_expect {
#------------------
    my $self = shift;

    if ($self->_layout == 2) {
      return $self->{'_lowestSignif'} || -1.0;
    }

    if($self->{'_is_significant'}) {
	my $bestHit = $self->{'_hits'}->[0];
	return $bestHit->expect();
    } else {
	$self->throw("Can't get lowest expect value: no significant hits ",
		     "The format of this report requires expect values to be extracted$Newline".
		     "from the hits themselves.");
    }
}

=head2 highest_p

 Example   : $blast->highest_p( ['overall'])
 Purpose   : Get the highest P-value among all hits in a BLAST report.
           : Syntactic sugar for $blast->hit('worst')->p()
           : Can also get the highest P-value overall (not just among signif hits).
 Usage     : $p_signif = $blast->highest_p();
           : $p_all    = $blast->highest_p('overall');
 Returns   : Float or scientific notation number.
           : Returns -1.0 if highest_p has not been set.
 Argument  : String 'overall' or no argument.
           : No argument = get highest P-value among significant hits.
 Throws    : Exception if object is created from a Blast2 report
           : (which does not report P-values).

See Also   : L<highest_signif()|highest_signif>, L<lowest_p()|lowest_p>, L<_set_signif()|_set_signif>, L<signif_fmt()|signif_fmt>

=cut

#---------------
sub highest_p {
#---------------
    my ($self, $overall) = @_;

    # Layout 2 = NCBI Blast 2.x does not report P-values.
    $self->_layout == 2 and
	$self->throw("Can't get P-value with BLAST2.",
		     "Use highest_signif() or highest_expect()");

    $overall and  return $self->{'_highestSignif'} || -1.0;
    $self->hit('worst')->p();
}

=head2 highest_expect

 Usage     : $blast_object->highest_expect( ['overall'])
 Purpose   : Get the highest Expect value among all significant hits in a BLAST report.
           : Syntactic sugar for $blast->hit('worst')->expect()
 Examples  : $e_sig = $blast->highest_expect();
           : $e_all = $blast->highest_expect('overall');
 Returns   : Float or scientific notation number.
           : Returns -1.0 if highest_exoect has not been set.
 Argument  : String 'overall' or no argument.
           : No argument = get highest Expect-value among significant hits.
 Throws    : Exception if there were no significant hits and the report
           : does not have Expect values on the description lines
           : (i.e., Blast1, WashU-Blast2).

See Also   : L<lowest_expect()|lowest_expect>, L<highest_signif()|highest_signif>, L<signif_fmt()|signif_fmt>

=cut

#-------------------
sub highest_expect {
#-------------------
    my ($self, $overall) = @_;

    if ( $overall and $self->_layout == 2) {
      return $self->{'_highestSignif'} || -1.0;
    }

    if($self->{'_is_significant'}) {
	return $self->hit('worst')->expect;
    } else {
	$self->throw("Can't get highest expect value: no significant hits ",
		     "The format of this report requires expect values to be extracted$Newline".
		     "from the hits themselves.");
    }
}

=head2 lowest_signif

 Usage     : $blast_obj->lowest_signif();
           : Syntactic sugar for $blast->hit('best')->signif()
 Purpose   : Get the lowest P or Expect value among all hits
           : in a BLAST report.
           : This method is syntactic sugar for $blast->hit('best')->signif()
           : The value returned is the one which is reported in the decription
           : section of the Blast report.
           : For Blast1 and WU-Blast2, this is a P-value,
           : for NCBI Blast2, it is an Expect value.
 Example   : $blast->lowest_signif();
 Returns   : Float or scientific notation number.
           : Returns -1.0 if lowest_signif has not been set.
 Argument  : n/a.
 Throws    : n/a.
 Status    : Deprecated. Use lowest_expect() or lowest_p().
 Comments  : The signif() method provides a way to deal with the fact that
           : Blast1 and Blast2 formats differ in what is reported in the
           : description lines of each hit in the Blast report. The signif()
           : method frees any client code from having to know if this is a P-value
           : or an Expect value, making it easier to write code that can process
           : both Blast1 and Blast2 reports. This is not necessarily a good thing, since
           : one should always know when one is working with P-values or
           : Expect values (hence the deprecated status).
           : Use of lowest_expect() is recommended since all hits will have an Expect value.

See Also   : L<lowest_p()|lowest_p>, L<lowest_expect()|lowest_expect>, L<signif()|signif>, L<signif_fmt()|signif_fmt>, L<_set_signif()|_set_signif>

=cut

#------------------
sub lowest_signif {
#------------------
    my ($self) = @_;

    return $self->{'_lowestSignif'} || -1.0;
}

=head2 highest_signif

 Usage     : $blast_obj->highest_signif('overall');
           : Syntactic sugar for $blast->hit('worst')->signif()
 Purpose   : Get the highest P or Expect value among all hits
           : in a BLAST report.
           : The value returned is the one which is reported in the decription
           : section of the Blast report.
           : For Blast1 and WU-Blast2, this is a P-value,
           : for NCBI Blast2, it is an Expect value.
 Example   : $blast->highest_signif();
 Returns   : Float or scientific notation number.
           : Returns -1.0 if highest_signif has not been set.
 Argument  : Optional  string 'overall' to get the highest overall significance value.
 Throws    : n/a.
 Status    : Deprecated. Use highest_expect() or highest_p().
 Comments  : Analogous to lowest_signif(), q.v.

See Also   : L<lowest_signif()|lowest_signif>, L<lowest_p()|lowest_p>, L<lowest_expect()|lowest_expect>, L<signif()|signif>, L<signif_fmt()|signif_fmt>, L<_set_signif()|_set_signif>

=cut

#---------------------
sub highest_signif {
#---------------------
    my ($self, $overall) = @_;

    $overall and return $self->{'_highestSignif'} || -1.0;

    if($self->{'_is_significant'}) {
        my $worst_hit = $self->hit('worst');
        if(defined $worst_hit) {
            return $worst_hit->signif;
        } else {
            return $self->{'_highestSignif'};
        }
    }
}

=head2 matrix

 Usage     : $blast_object->matrix();
 Purpose   : Get the name of the scoring matrix used.
           : This is extracted from the report.
 Argument  : n/a
 Returns   : string or undef if not defined

=cut

#------------
sub matrix { my $self = shift; $self->{'_matrix'} || $Blast->{'_matrix'}; }
#------------

=head2 filter

 Usage     : $blast_object->filter();
 Purpose   : Get the name of the low-complexity sequence filter used.
           : (SEG, SEG+XNU, DUST, NONE).
           : This is extracted from the report.
 Argument  : n/a
 Returns   : string or undef if not defined

=cut

#----------
sub filter { my $self = shift; $self->{'_filter'}  || $Blast->{'_filter'}; }
#----------

=head2 expect

 Usage     : $blast_object->expect();
 Purpose   : Get the expect parameter (E) used for the Blast analysis.
           : This is extracted from the report.
 Argument  : n/a
 Returns   : string or undef if not defined.

=cut

#-----------
sub expect { my $self = shift; $self->{'_expect'} || $Blast->{'_expect'}; }
#-----------

=head2 karlin_altschul

 Usage     : $blast_object->karlin_altschul();
 Purpose   : Get the Karlin_Altschul sum statistics (Lambda, K, H)
           : These are extracted from the report.
 Argument  : n/a
 Returns   : list of three floats (Lambda, K, H)
           : If not defined, returns list of three zeros)

=cut

#---------------------
sub karlin_altschul {
#---------------------
    my $self = shift;
    if(defined($self->{'_lambda'})) {
	($self->{'_lambda'}, $self->{'_k'}, $self->{'_h'});
    } elsif(defined($Blast->{'_lambda'})) {
	($Blast->{'_lambda'}, $Blast->{'_k'}, $Blast->{'_h'});
    } else {
	(0, 0, 0);
    }
}

=head2 word_size

 Usage     : $blast_object->word_size();
 Purpose   : Get the word_size used during the Blast analysis.
           : This is extracted from the report.
 Argument  : n/a
 Returns   : integer or undef if not defined.

=cut

#--------------
sub word_size {
#--------------
    my $self = shift;
    $self->{'_word_size'} || $Blast->{'_word_size'};
}

=head2 s

 Usage     : $blast_object->s();
 Purpose   : Get the s statistic for the Blast analysis.
           : This is extracted from the report.
 Argument  : n/a
 Returns   : integer or undef if not defined.

=cut

#------
sub s { my $self = shift; $self->{'_s'} || $Blast->{'_s'}; }
#------

=head2 gap_creation

 Usage     : $blast_object->gap_creation();
 Purpose   : Get the gap creation penalty used for a gapped Blast analysis.
           : This is extracted from the report.
 Argument  : n/a
 Returns   : integer or undef if not defined.

See Also   : L<gap_extension()|gap_extension>

=cut

#-----------------
sub gap_creation {
#-----------------
    my $self = shift;
    $self->{'_gapCreation'} || $Blast->{'_gapCreation'};
}

=head2 gap_extension

 Usage     : $blast_object->gap_extension();
 Purpose   : Get the gap extension penalty used for a gapped Blast analysis.
           : This is extracted from the report.
 Argument  : n/a
 Returns   : integer or undef if not defined.

See Also   : L<gap_extension()|gap_extension>

=cut

#-------------------
sub gap_extension {
#-------------------
    my $self = shift;
    $self->{'_gapExtension'} || $Blast->{'_gapExtension'};
}

=head2 ambiguous_aln

 Usage     : $blast_object->ambiguous_aln();
 Purpose   : Test all hits and determine if any have an ambiguous alignment.
 Example   : print "ambiguous" if $blast->ambiguous_aln();
 Returns   : Boolean (true if ANY significant hit has an ambiguous alignment)
 Argument  : n/a
 Throws    : n/a
 Status    : Experimental
 Comments  : An ambiguous BLAST alignment is defined as one where two or more
           : different HSPs have significantly overlapping sequences such
           : that it is not possible to create a unique alignment
           : by simply concatenating HSPs. This may indicate the presence
           : of multiple domains in one sequence relative to another.
           : This method only indicates the presence of ambiguity in at
           : least one significant hit. To determine the nature of the
           : ambiguity, each hit must be examined.

See Also   : B<Bio::Tools::Blast::Sbjct::ambiguous_aln()>,L<Links to related modules>

=cut

#----------------
sub ambiguous_aln {
#----------------
    my $self = shift;
    foreach($self->hits()) {
	return 1 if ($_->ambiguous_aln() ne '-');
    }
    0;
}

=head2 overlap

 Usage     : $blast_object->overlap([integer]);
 Purpose   : Set/Get the number of overlapping residues allowed when tiling multiple HSPs.
           : Delegates to Bio::Tools::Blast::Sbjct::overlap().
 Throws    : Exception if there are no significant hits.
 Status    : Experimental

See Also   : B<Bio::Tools::Blast::Sbjct::overlap()>,L<Links to related modules>

=cut

#------------
sub overlap {
#------------
    my $self = shift;
    if(not $self->hits) {
	$self->throw("Can't get overlap data without significant hits.");
    }
    $self->hit->overlap();
}

=head2 homol_data

 Usage     : @data = $blast_object->homo_data( %named_params );
 Purpose   : Gets specific similarity data about each significant hit.
 Returns   : Array of strings:
           : "Homology data" for each HSP is in the format:
           :  "<integer> <start> <stop>"
           : Data for different HSPs are tab-delimited.
 Argument  : named parameters passed along to the hit objects.
 Throws    : n/a
 Status    : Experimental
 Comments  : This is a very experimental method used for obtaining an
           : indication of:
           :   1) how many HSPs are in a Blast alignment
           :   2) how strong the similarity is between sequences in the HSP
           :   3) the endpoints of the alignment (sequence monomer numbers)

See Also   : B<Bio::Tools::Blast::Sbjct::homol_data()>,L<Links to related modules>

=cut

#----------------
sub homol_data {
#----------------

    my ($self, %param) = @_;
    my @hits = $self->hits();
    my @data = ();

    ## Note: Homology data can be either for the query sequence or the hit
    ##       (Sbjct) sequence. Default is for sbjct. This is specifyable via
    ##       $param{-SEQ}='sbjct' || 'query'.

    foreach ( @hits ) {
	push @data, $_->homol_data(%param);
    }
    @data;
}

=head1 REPORT GENERATING METHODS

=head2 table

 Usage     : $blast_obj->table( [get_desc]);
 Purpose   : Output data for each HSP of each hit in tab-delimited format.
 Example   : print $blast->table;
           : print $blast->table(0);
           : # Call table_labels() to print labels.
 Argument  : get_desc = boolean, if false the description of each hit is not included.
           :            Default: true (if not defined, include description column).
 Returns   : String containing tab-delimited set of data for each HSP
           : of each significant hit. Different HSPs are separated by newlines.
           : Left-to-Right order of fields:
           : 1 QUERY_NAME             # Sequence identifier of the query.
           : 2 QUERY_LENGTH           # Full length of the query sequence.
           : 3 SBJCT_NAME             # Sequence identifier of the sbjct ("hit".
           : 4 SBJCT_LENGTH           # Full length of the sbjct sequence.
           : 5 EXPECT                 # Expect value for the alignment.
           : 6 SCORE                  # Blast score for the alignment.
           : 7 BITS                   # Bit score for the alignment.
           : 8 NUM_HSPS               # Number of HSPs (not the "N" value).
           : 9 HSP_FRAC_IDENTICAL     # fraction of identical substitutions.
           : 10 HSP_FRAC_CONSERVED    # fraction of conserved ("positive") substitutions.
           : 11 HSP_QUERY_ALN_LENGTH  # Length of the aligned portion of the query sequence.
           : 12 HSP_SBJCT_ALN_LENGTH  # Length of the aligned portion of the sbjct sequence.
           : 13 HSP_QUERY_GAPS        # Number of gaps in the aligned query sequence.
           : 14 HSP_SBJCT_GAPS        # Number of gaps in the aligned sbjct sequence.
           : 15 HSP_QUERY_START       # Starting coordinate of the query sequence.
           : 16 HSP_QUERY_END         # Ending coordinate of the query sequence.
           : 17 HSP_SBJCT_START       # Starting coordinate of the sbjct sequence.
           : 18 HSP_SBJCT_END         # Ending coordinate of the sbjct sequence.
           : 19 HSP_QUERY_STRAND      # Strand of the query sequence (TBLASTN/X only)
           : 20 HSP_SBJCT_STRAND      # Strand of the sbjct sequence (TBLASTN/X only)
           : 21 HSP_FRAME             # Frame for the sbjct translation (TBLASTN/X only)
           : 22 SBJCT_DESCRIPTION  (optional)  # Full description of the sbjct sequence from
           :                                  # the alignment section.
 Throws    : n/a
 Comments  : This method does not collect data based on tiling of the HSPs.
           : The table will contains redundant information since the hit name,
           : id, and other info for the hit are listed for each HSP.
           : If you need more flexibility in the output format than this
           : method provides, design a custom function.

See Also  : L<table_tiled()|table_tiled>, L<table_labels()|table_labels>, L<_display_hits()|_display_hits>

=cut

#-----------
sub table {
#-----------
    my ($self, $get_desc) = @_;
    my $str = '';

    $get_desc = defined($get_desc) ? $get_desc : 1;
#    $str .= $self->_table_labels($get_desc) unless $self->{'_labels'};

    my $sigfmt = $self->signif_fmt();
    $sigfmt eq 'parts' and $sigfmt = 'exp';  # disallow 'parts' format for this table.
    my $sigprint = $sigfmt eq 'exp' ? 'd' : '.1e';

    my ($hit, $hsp);
    foreach $hit($self->hits) {
	foreach $hsp($hit->hsps) {
	    # Note: range() returns a 2-element list.
	    $str .= sprintf "%s\t%d\t%s\t%d\t%$sigprint\t%d\t%d\t%d\t%.2f\t%.2f\t%d\t%d\t%d\t%d\t%d\t%d\t%d\t%d\t%s\t%s\t%s\t%s$Newline",
                   $self->name, $self->length, $hit->name, $hit->length,
	           $hit->expect($sigfmt), $hit->score, $hit->bits,
 	           $hit->num_hsps, $hsp->frac_identical, $hsp->frac_conserved,
	           $hsp->length('query'), $hsp->length('sbjct'),
	           $hsp->gaps('list'),
	           $hsp->range('query'), $hsp->range('sbjct'),
	           $hsp->strand('query'), $hsp->strand('sbjct'), $hsp->frame,
	           ($get_desc ? $hit->desc  : '');
	}
    }
    $str =~ s/\t$Newline/$Newline/gs;
    $str;
}

=head2 table_labels

 Usage     : print $blast_obj->table_labels( [get_desc] );
 Purpose   : Get column labels for table().
 Returns   : String containing column labels. Tab-delimited.
 Argument  : get_desc = boolean, if false the description column is not included.
           : Default: true (if not defined, include description column).
 Throws    : n/a

See Also   : L<table()|table>

=cut

#----------------
sub table_labels {
#----------------
    my ($self, $get_desc) = @_;
    $get_desc = defined($get_desc) ? $get_desc : 1;
    my $descstr = $get_desc ? 'DESC' : '';
    my $descln = $get_desc ? '-----' : '';

    my $str = sprintf "%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s$Newline",
                       'QUERY', 'Q_LEN', 'SBJCT', 'S_LEN', 'EXPCT', 'SCORE', 'BITS', 'HSPS',
                       'IDEN', 'CONSV', 'Q_ALN', 'S_ALN', 'Q_GAP', 'S_GAP',
                       'Q_BEG', 'Q_END', 'S_BEG', 'S_END', 'Q_STR', 'S_STR', 'FRAM', $descstr;
    $str .= sprintf "%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s$Newline",
                       '-----', '-----', '-----', '-----', '-----', '-----', '-----', '-----',
                       '-----', '-----', '-----', '-----', '-----', '-----',
                       '-----', '-----', '-----','-----', '-----', '-----','-----', $descln;

    $self->{'_labels'} = 1;
    $str =~ s/\t$Newline/$Newline/gs;
    $str;
}

=head2 table_tiled

 Purpose   : Get data from tiled HSPs in tab-delimited format.
           : Allows only minimal flexibility in the output format.
           : If you need more flexibility, design a custom function.
 Usage     : $blast_obj->table_tiled( [get_desc]);
 Example   : print $blast->table_tiled;
           : print $blast->table_tiled(0);
           : # Call table_labels_tiled() if you want labels.
 Argument  : get_desc = boolean, if false the description of each hit is not included.
           :            Default: true (include description).
 Returns   : String containing tab-delimited set of data for each HSP
           : of each significant hit. Multiple hits are separated by newlines.
           : Left-to-Right order of fields:
           : 1 QUERY_NAME           # Sequence identifier of the query.
           : 2 QUERY_LENGTH         # Full length of the query sequence.
           : 3 SBJCT_NAME           # Sequence identifier of the sbjct ("hit".
           : 4 SBJCT_LENGTH         # Full length of the sbjct sequence.
           : 5 EXPECT               # Expect value for the alignment.
           : 6 SCORE                # Blast score for the alignment.
           : 7 BITS                 # Bit score for the alignment.
           : 8 NUM_HSPS             # Number of HSPs (not the "N" value).
           : 9 FRAC_IDENTICAL*      # fraction of identical substitutions.
           : 10 FRAC_CONSERVED*     # fraction of conserved ("positive") substitutions .
           : 11 FRAC_ALN_QUERY*     # fraction of the query sequence that is aligned.
           : 12 FRAC_ALN_SBJCT*     # fraction of the sbjct sequence that is aligned.
           : 13 QUERY_ALN_LENGTH*   # Length of the aligned portion of the query sequence.
           : 14 SBJCT_ALN_LENGTH*   # Length of the aligned portion of the sbjct sequence.
           : 15 QUERY_GAPS*         # Number of gaps in the aligned query sequence.
           : 16 SBJCT_GAPS*         # Number of gaps in the aligned sbjct sequence.
           : 17 QUERY_START*        # Starting coordinate of the query sequence.
           : 18 QUERY_END*          # Ending coordinate of the query sequence.
           : 19 SBJCT_START*        # Starting coordinate of the sbjct sequence.
           : 20 SBJCT_END*          # Ending coordinate of the sbjct sequence.
           : 21 AMBIGUOUS_ALN       # Ambiguous alignment indicator ('qs', 'q', 's').
           : 22 SBJCT_DESCRIPTION  (optional)  # Full description of the sbjct sequence from
           :                                  # the alignment section.
           :
           : * Items marked with a "*" report data summed across all HSPs
           :   after tiling them to avoid counting data from overlapping regions
           :   multiple times.
 Throws    : n/a
 Comments  : This function relies on tiling of the HSPs since it calls
           : frac_identical() etc. on the hit as opposed to each HSP individually.

See Also   : L<table()|table>, L<table_labels_tiled()|table_labels_tiled>, B<Bio::Tools::Blast::Sbjct::"HSP Tiling and Ambiguous Alignments">, L<Links to related modules>

=cut

#----------------
sub table_tiled {
#----------------
    my ($self, $get_desc) = @_;
    my $str = '';

    $get_desc = defined($get_desc) ? $get_desc : 1;

    my ($hit);
    my $sigfmt = $self->signif_fmt();
    $sigfmt eq 'parts' and $sigfmt = 'exp';  # disallow 'parts' format for this table.
    my $sigprint = $sigfmt eq 'exp' ? 'd' : '.1e';

    foreach $hit($self->hits) {
	$str .= sprintf "%s\t%d\t%s\t%d\t%$sigprint\t%d\t%d\t%d\t%.2f\t%.2f\t%.2f\t%.2f\t%d\t%d\t%d\t%d\t%d\t%d\t%d\t%d\t%s\t%s$Newline",
	               $self->name, $self->length, $hit->name, $hit->length,
	               $hit->expect($sigfmt), $hit->score, $hit->bits,
	               $hit->num_hsps, $hit->frac_identical, $hit->frac_conserved,
	               $hit->frac_aligned_query, $hit->frac_aligned_hit,
	               $hit->length_aln('query'), $hit->length_aln('sbjct'),
                       $hit->gaps('list'), $hit->range('query'), $hit->range('sbjct'),
	               $hit->ambiguous_aln, ($get_desc ? $hit->desc : '');
    }
    $str =~ s/\t$Newline/$Newline/gs;
    $str;
}

=head2 table_labels_tiled

 Usage     : print $blast_obj->table_labels_tiled( [get_desc] );
 Purpose   : Get column labels for table_tiled().
 Returns   : String containing column labels. Tab-delimited.
 Argument  : get_desc = boolean, if false the description column is not included.
           : Default: true (include description column).
 Throws    : n/a

See Also   : L<table_tiled()|table_tiled>

=cut

#---------------------
sub table_labels_tiled {
#---------------------
    my ($self, $get_desc) = @_;
    my $descstr = $get_desc ? 'DESC' : '';
    my $descln = $get_desc ? '-----' : '';

    my $str = sprintf "%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s$Newline",
                       'QUERY', 'Q_LEN', 'SBJCT', 'S_LEN', 'EXPCT', 'SCORE', 'BITS',
                       'HSPS', 'FR_ID', 'FR_CN', 'FR_ALQ', 'FR_ALS', 'Q_ALN',
                       'S_ALN', 'Q_GAP', 'S_GAP', 'Q_BEG', 'Q_END', 'S_BEG', 'S_END',
                       'AMBIG', $descstr;
    $str =~ s/\t$Newline/$Newline/;
    $str .= sprintf "%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s$Newline",
                       '-----', '-----', '------', '-----', '-----','-----', '-----',
                       '-----', '-----', '-----', '-----', '-----', '-----',
                       '-----', '-----', '-----','-----','-----', '-----',
                       '-----','-----', $descln;

    $self->{'_labels_tiled'} = 1;
    $str =~ s/\t$Newline/$Newline/gs;
    $str;
}

=head2 display

 Usage     : $blast_object->display( %named_parameters );
 Purpose   : Display information about Bio::Tools::Blast.pm data members,
           : E.g., parameters of the report, data for each hit., etc.
           : Overrides Bio::Root::Object::display().
 Example   : $object->display(-SHOW=>'stats');
           : $object->display(-SHOW=>'stats,hits');
 Argument  : Named parameters: (TAGS CAN BE UPPER OR LOWER CASE)
           :     -SHOW  => 'file' | 'hits' | 'homol'
           :     -WHERE => filehandle (default = STDOUT)
 Returns   : n/a (print/printf is called)
 Status    : Experimental
 Comments  : For tab-delimited output, see table().

See Also   : L<_display_homol()|_display_homol>, L<_display_hits()|_display_hits>, L<_display_stats()|_display_stats>, L<table()|table>, B<Bio::Root::Tools::SeqAnal::display()>,L<Links to related modules>,

=cut

#--------------
sub display {
#--------------
    my( $self, %param ) = @_;

    $self->SUPER::display(%param);
    my $OUT = $self->fh();

    $self->show =~ /homol/i and $self->_display_homol($OUT);
    $self->show =~ /hits/i and $self->_display_hits( %param );
    1;
}

=head2 _display_homol

 Usage     : n/a; called automatically by display()
 Purpose   : Print homology data for hits in the BLAST report.
 Example   : n/a
 Argument  : one argument = filehandle object.
 Returns   : printf call.
 Status    : Experimental

See Also   : L<homol_data()|homol_data>, L<display()|display>

=cut

#-------------------
sub _display_homol {
#-------------------
    my( $self, $OUT ) = @_;

    print $OUT "${Newline}BLAST HOMOLOGY DATA FOR: ${\$self->name()}$Newline";
    print $OUT '-'x40,"$Newline";

    foreach ( $self->homol_data()) {
	print $OUT "$_$Newline";
    }
}

=head2 _display_stats

 Usage     : n/a; called automatically by display()
 Purpose   : Display information about the Blast report "meta" data.
           : Overrides Bio::Tools::SeqAnal::_display_stats() calling it first.
 Example   : n/a
 Argument  : one argument = filehandle object.
 Returns   : printf call.
 Status    : Experimental

See Also   : L<display()|display>, B<Bio::Tools::SeqAnal::_display_stats()>,L<Links to related modules>

=cut

#--------------------
sub _display_stats {
#--------------------
    my( $self, $OUT ) = @_;

    $self->SUPER::_display_stats($OUT);
    printf( $OUT "%-15s: %s$Newline", "GAPPED", $self->gapped ? 'YES' : 'NO');
    printf( $OUT "%-15s: %d$Newline", "TOTAL HITS", $self->num_hits('total'));
    printf( $OUT "%-15s: %s$Newline", "CHECKED ALL", $Blast->{'_check_all'} ? 'YES' : 'NO');
    printf( $OUT "%-15s: %s$Newline", "FILT FUNC", $Blast->{'_filt_func'} ? 'YES' : 'NO');
    if($self->min_length) {
	printf( $OUT "%-15s: Length >= %s$Newline", "MIN_LENGTH", $self->min_length);
    }

    my $num_hits =  $self->num_hits;
    my $signif_str = ($self->_layout == 1) ? 'P' : 'EXPECT';

    printf( $OUT "%-15s: %d$Newline", "SIGNIF HITS", $num_hits);
    # Blast1: signif = P-value, Blast2: signif = Expect value.
	
    printf( $OUT "%-15s: %s ($signif_str-VALUE)$Newline", "SIGNIF CUTOFF", $self->signif);
    printf( $OUT "%-15s: %s$Newline", "LOWEST $signif_str", $self->lowest_signif());
    printf( $OUT "%-15s: %s$Newline", "HIGHEST $signif_str", $self->highest_signif());

    printf( $OUT "%-15s: %s (OVERALL)$Newline", "HIGHEST $signif_str", $self->highest_signif('overall'));

    if($Blast->_get_stats) {
	my $warn = ($Blast->{'_share'}) ? '(SHARED STATS)' : '';
	printf( $OUT "%-15s: %s$Newline", "MATRIX", $self->matrix() || 'UNKNOWN');
	printf( $OUT "%-15s: %s$Newline", "FILTER", $self->filter() || 'UNKNOWN');
	printf( $OUT "%-15s: %s$Newline", "EXPECT", $self->expect() || 'UNKNOWN');
	printf( $OUT "%-15s: %s, %s, %s %s$Newline", "LAMBDA, K, H", $self->karlin_altschul(), $warn);
	printf( $OUT "%-15s: %s$Newline", "WORD SIZE", $self->word_size() || 'UNKNOWN');
	printf( $OUT "%-15s: %s %s$Newline", "S", $self->s() || 'UNKNOWN', $warn);
	if($self->gapped) {
	    printf( $OUT "%-15s: %s$Newline", "GAP CREATION", $self->gap_creation() || 'UNKNOWN');
	    printf( $OUT "%-15s: %s$Newline", "GAP EXTENSION", $self->gap_extension() || 'UNKNOWN');
	}
    }
    print $OUT "$Newline";
}

=head2 _display_hits

 Usage     : n/a; called automatically by display()
 Purpose   : Display data for each hit. Not tab-delimited.
 Example   : n/a
 Argument  : one argument = filehandle object.
 Returns   : printf call.
 Status    : Experimental
 Comments  : For tab-delimited output, see table().

See Also   : L<display()|display>, B<Bio::Tools::Blast::Sbjct::display()>, L<table()|table>,  L<Links to related modules>

=cut

sub _display_hits {

    my( $self, %param ) = @_;
    my $OUT = $self->fh();
    my @hits  = $self->hits();

    ## You need a wide screen to see this properly.
    # Header.
    print $OUT "${Newline}BLAST HITS FOR: ${\$self->name()} length = ${\$self->length}$Newline";
    print "(This table requires a wide display.)$Newline";
    print $OUT '-'x80,"$Newline";

    print $self->table_labels_tiled(0);
    print $self->table_tiled(0);

    ## Doing this interactively since there is potentially a lot of data here.
    ## Not quite satisfied with this approach.

    if (not $param{-INTERACTIVE}) {
	return 1;
    } else {
	my ($reply);
	print "${Newline}DISPLAY FULL HSP DATA? (y/n): [n] ";
	chomp( $reply = <STDIN> );
	$reply =~ /^y.*/i;
	
	my $count = 0;
	foreach ( @hits ) {
	    $count++;
	    print $OUT "$Newline$Newline",'-'x80,"$Newline";
	    print $OUT "HSP DATA FOR HIT #$count  (hit <RETURN>)";
		print $OUT "$Newline",'-'x80;<STDIN>;
	    $param{-SHOW} = 'hsp';
	    $_->display( %param );
	}
    }
    1;
}

=head2 to_html

 Usage     : $blast_object->to_html( [%named_parameters] )
 Purpose   : To produce an HTML-formatted version of a BLAST report
           : for efficient navigation of the report using a web browser.
 Example   : # Using the static Blast object:
           : # Can read from STDIN or from a designated file:
           :   $Blast->to_html($file);
           :   $Blast->to_html(-FILE=>$file, -HEADER=>$header);
           :   (if no file is supplied, STDIN will be used).
           : # saving HTML to an array:
           :   $Blast->to_html(-FILE=>$file, -OUT =>\@out);
           : # Using a pre-existing blast object (must have been built from
           : # a file, not STDIN:
           :   $blastObj->to_html();
 Returns   : n/a, either prints report to STDOUT or saves to a supplied array
           : if an '-OUT' parameter is defined (see below).
 Argument  : %named_parameters: (TAGS ARE AND CASE INSENSITIVE).
           :    -FILE   => string containing name of a file to be processed.
           :               If not a valid file or undefined, STDIN will be used.
           :               Can skip the -FILE tag if supplying a filename
           :               as a single argument.
           :    -HEADER => string
           :               This should be an HTML-formatted string to be used
           :               as a header for the page, typically describing query sequence,
           :               database searched, the date of the analysis, and any
           :               additional links.
           :               If not supplied, no special header is used.
           :               Regardless of whether a header is supplied, the
           :               standard info at the top of the report is highlighted.
           :               This should include the <HEADER></HEADER> section
           :               of the page as well.
           :
           :    -IN    => array reference containing a raw Blast report.
           :              each line in a separate element in the array.
           :              If -IN is not supplied, read() is called
           :              and data is then read either from STDIN or a file.
           :
           :    -OUT   => array reference to hold the HTML output.
           :              If not supplied, output is sent to STDOUT.
 Throws    : Exception is propagated from $HTML::get_html_func()
           : and Bio::Root::Object::read().
 Comments  : The code that does the actual work is located in
           :  Bio::Tools::Blast::HTML::get_html_func().
 Bugs      : Some hypertext links to external databases may not be
           : correct. This due in part to the dynamic nature of
           : the web.
           : Hypertext links are not added to hits without database ids.
 TODO      : Possibly create a function to produce fancy default header
           : using data extracted from the report (requires some parsing).
           : For example, it would be nice to always include a date

See Also   : B<Bio::Tools::Blast::HTML::get_html_func()>, B<Bio::Root::Object::read()>, L<Links to related modules>

=cut

#------------
sub to_html {
#------------
    my ($self, @param) = @_;

    # Permits syntax such as: $blast->to_html($filename);
    my ($file, $header_html, $in_aref, $out_aref) =
	$self->_rearrange([qw(FILE HEADER IN OUT)], @param);

    $self->file($file) if $file;

    # Only setting the newline character once for efficiency.
    $Newline ||= $Util->get_newline(-client => $self, @param);

    $header_html ||= '';
    (ref($out_aref) eq 'ARRAY') ? push(@$out_aref, $header_html) : print "$header_html$Newline";

    require Bio::Tools::Blast::HTML;
    Bio::Tools::Blast::HTML->import(qw(&get_html_func));

    my ($func);
    eval{ $func = &get_html_func($out_aref);  };
    if($@) {
	my $err = $@;
	$self->throw($err);
    }

    eval {
	if(!$header_html) {
	    $out_aref ? push(@$out_aref, "<html><body>$Newline") : print "<html><body>$Newline";
	}

	if (ref ($in_aref) =~ /ARRAY/) {
	    # If data is being supplied, process it.
	    foreach(@$in_aref) {
		&$func($_);
	    }
	} else {
	    # Otherwise, read it, processing as we go.
	
  	    $self->read(-FUNC => $func, @param);
	}
	$out_aref ? push(@$out_aref, "$Newline</pre></body></html>") : print "$Newline</pre></body></html>";
    };

    if($@) {
	# Check for trivial error (report already HTML formatted).
	if($@ =~ /HTML formatted/) {
	    print STDERR "\a${Newline}Blast report appears to be HTML formatted already.$Newline$Newline";
	} else {
	    my $err = $@;
	    $self->throw($err);
	}
    }
}

1;
__END__

#####################################################################################
#                                END OF CLASS                                       #
#####################################################################################

=head1 FOR DEVELOPERS ONLY

=head2 Data Members

Information about the various data members of this module is provided for those
wishing to modify or understand the code. Two things to bear in mind:

=over 4

=item 1 Do NOT rely on these in any code outside of this module.

All data members are prefixed with an underscore to signify that they are private.
Always use accessor methods. If the accessor doesn't exist or is inadequate,
create or modify an accessor (and let me know, too!). (An exception to this might
be for Sbjct.pm or HSP.pm which are more tightly coupled to Blast.pm and
may access Blast data members directly for efficiency purposes, but probably
should not).

=item 2 This documentation may be incomplete and out of date.

It is easy for these data member descriptions to become obsolete as
this module is still evolving. Always double check this info and search
for members not described here.

=back

An instance of Bio::Tools::Blast.pm is a blessed reference to a hash containing
all or some of the following fields:

 FIELD           VALUE
 --------------------------------------------------------------
 _significance    P-value or Expect value cutoff (depends on Blast version:
	          Blast1/WU-Blast2 = P-value; Blast2 = Expect value).
   	          Values GREATER than this are deemed not significant.

 _significant     Boolean. True if the query has one or more significant hit.

 _min_length      Integer. Query sequences less than this will be skipped.

 _confirm_significance  Boolean. True if client has supplied significance criteria.

 _gapped          Boolean. True if BLAST analysis has gapping turned on.

 _hits            List of Sbjct.pm objects.

 _num_hits        Number of hits obtained from the BLAST report.

 _num_hits_significant Number of significant based on Significant data members.

 _highestSignif   Highest P or Expect value overall (not just what is stored in _hits).

 _lowestSignif    Lowest P or Expect value overall (not just what is stored in _hits).

The static $Blast object has a special set of members:

  _errs
  _share
  _stream
  _get_stats
  _gapped
  _filt_func

 Miscellaneous statistical parameters:
 -------------------------------------
  _filter, _matrix, _word_size, _expect, _gapCreation, _gapExtension, _s,
  _lambda, _k, _h

 INHERITED DATA MEMBERS
 -----------------------
 (See Bio::Tools::SeqAnal.pm for inherited data members.)

=cut

1;