Bio/SearchDist.pm

   1 # $Id$
   2
   3 #
   4 # BioPerl module for Bio::SearchDist
   5 #
   6 # Please direct questions and support issues to <bioperl-l@bioperl.org>
   7 #
   8 # Cared for by Ewan Birney <birney@ebi.ac.uk>
   9 #
  10 # Copyright Ewan Birney
  11 #
  12 # You may distribute this module under the same terms as perl itself
  13
  14 # POD documentation - main docs before the code
  15
  16 =head1 NAME
  17
  18 Bio::SearchDist - A perl wrapper around Sean Eddy's histogram object
  19
  20 =head1 SYNOPSIS
  21
  22   $dis = Bio::SearchDist->new();
  23   foreach $score ( @scores ) {
  24      $dis->add_score($score);
  25   }
  26
  27   if( $dis->fit_evd() ) {
  28     foreach $score ( @scores ) {
  29       $evalue = $dis->evalue($score);
  30       print "Score $score had an evalue of $evalue\n";
  31     }
  32   } else {
  33     warn("Could not fit histogram to an EVD!");
  34   }
  35
  36 =head1 DESCRIPTION
  37
  38 The Bio::SearchDist object is a wrapper around Sean Eddy's excellent
  39 histogram object. The histogram object can bascially take in a number
  40 of scores which are sensibly distributed somewhere around 0 that come
  41 from a supposed Extreme Value Distribution. Having add all the scores
  42 from a database search via the add_score method you can then fit a
  43 extreme value distribution using fit_evd(). Once fitted you can then
  44 get out the evalue for each score (or a new score) using
  45 evalue($score).
  46
  47 The fitting procedure is better described in Sean Eddy's own code
  48 (available from http://hmmer.janelia.org/, or in the histogram.h header
  49 file in Compile/SW). Bascially it fits a EVD via a maximum likelhood
  50 method with pruning of the top end of the distribution so that real
  51 positives are discarded in the fitting procedure. This comes from
  52 an orginally idea of Richard Mott's and the likelhood fitting
  53 is from a book by Lawless [should ref here].
  54
  55
  56 The object relies on the fact that the scores are sensibly distributed
  57 around about 0 and that integer bins are sensible for the
  58 histogram. Scores based on bits are often ideal for this (bits based
  59 scoring mechanisms is what this histogram object was originally
  60 designed for).
  61
  62
  63 =head1 CONTACT
  64
  65 The original code this was based on comes from the histogram module as
  66 part of the HMMer2 package. Look at http://hmmer.janelia.org/
  67
  68 Its use in Bioperl is via the Compiled XS extension which is cared for
  69 by Ewan Birney (birney@ebi.ac.uk). Please contact Ewan first about
  70 the use of this module
  71
  72 =head1 FEEDBACK
  73
  74 =head2 Mailing Lists
  75
  76 User feedback is an integral part of the evolution of this and other
  77 Bioperl modules. Send your comments and suggestions preferably to one
  78 of the Bioperl mailing lists.  Your participation is much appreciated.
  79
  80   bioperl-l@bioperl.org                  - General discussion
  81   http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
  82
  83 =head2 Support
  84
  85 Please direct usage questions or support issues to the mailing list:
  86
  87 I<bioperl-l@bioperl.org>
  88
  89 rather than to the module maintainer directly. Many experienced and
  90 reponsive experts will be able look at the problem and quickly
  91 address it. Please include a thorough description of the problem
  92 with code and data examples if at all possible.
  93
  94 =head2 Reporting Bugs
  95
  96 Report bugs to the Bioperl bug tracking system to help us keep track
  97 the bugs and their resolution.  Bug reports can be submitted via the
  98 web:
  99
 100   http://bugzilla.open-bio.org/
 101
 102 =head1 APPENDIX
 103
 104 The rest of the documentation details each of the object
 105 methods. Internal methods are usually preceded with a _
 106
 107 =cut
 108
 109
 110 # Let the code begin...
 111
 112
 113 package Bio::SearchDist;
 114 use strict;
 115
 116
 117 BEGIN {
 118     eval {
 119         require Bio::Ext::Align;
 120     };
 121     if ( $@ ) {
 122 print $@;
 123         print STDERR ("\nThe C-compiled engine for histogram object (Bio::Ext::Align) has not been installed.\n Please install the bioperl-ext package\n\n");
 124         exit(1);
 125     }
 126 }
 127
 128
 129 use base qw(Bio::Root::Root);
 130
 131 sub new {
 132   my($class,@args) = @_;
 133   my $self = $class->SUPER::new(@args);
 134   my($min, $max, $lump) =
 135         $self->_rearrange([qw(MIN MAX LUMP)], @args);
 136
 137   if( !  $min ) {
 138     $min = -100;
 139   }
 140
 141   if( !  $max ) {
 142     $max = +100;
 143   }
 144
 145   if( ! $lump ) {
 146     $lump = 50;
 147   }
 148
 149   $self->_engine(&Bio::Ext::Align::new_Histogram($min,$max,$lump));
 150
 151   return $self;
 152 }
 153
 154 =head2 add_score
 155
 156  Title   : add_score
 157  Usage   : $dis->add_score(300);
 158  Function: Adds a single score to the distribution
 159  Returns : nothing
 160  Args    :
 161
 162
 163 =cut
 164
 165 sub add_score{
 166    my ($self,$score) = @_;
 167    my ($eng);
 168    $eng = $self->_engine();
 169    #$eng->AddToHistogram($score);
 170    $eng->add($score);
 171 }
 172
 173 =head2 fit_evd
 174
 175  Title   : fit_evd
 176  Usage   : $dis->fit_evd();
 177  Function: fits an evd to the current distribution
 178  Returns : 1 if it fits successfully, 0 if not
 179  Args    :
 180
 181
 182 =cut
 183
 184 sub fit_evd{
 185    my ($self,@args) = @_;
 186
 187    return $self->_engine()->fit_EVD(10000,1);
 188 }
 189
 190 =head2 fit_Gaussian
 191
 192  Title   : fit_Gaussian
 193  Usage   :
 194  Function:
 195  Example :
 196  Returns :
 197  Args    :
 198
 199
 200 =cut
 201
 202 sub fit_Gaussian{
 203    my ($self,$high) = @_;
 204
 205    if( ! defined $high ) {
 206        $high = 10000;
 207    }
 208
 209    return $self->_engine()->fit_Gaussian($high);
 210 }
 211
 212
 213 =head2 evalue
 214
 215  Title   : evalue
 216  Usage   : $eval = $dis->evalue($score)
 217  Function: Returns the evalue of this score
 218  Returns : float
 219  Args    :
 220
 221
 222 =cut
 223
 224 sub evalue{
 225    my ($self,$score) = @_;
 226
 227    return $self->_engine()->evalue($score);
 228
 229 }
 230
 231
 232
 233 =head2 _engine
 234
 235  Title   : _engine
 236  Usage   : $obj->_engine($newval)
 237  Function: underlyine bp_sw:: histogram engine
 238  Returns : value of _engine
 239  Args    : newvalue (optional)
 240
 241
 242 =cut
 243
 244 sub _engine{
 245     my ($self,$value) = @_;
 246     if( defined $value) {
 247         $self->{'_engine'} = $value;
 248     }
 249     return $self->{'_engine'};
 250 }
 251
 252
 253 ## End of Package
 254
 255 1;
 256 __END__
 257