Bio/SearchDist.pm

   1
   2
   3 #
   4 # BioPerl module for Bio::SearchDist
   5 #
   6 # Cared for by Ewan Birney <birney@sanger.ac.uk>
   7 #
   8 # Copyright Ewan Birney
   9 #
  10 # You may distribute this module under the same terms as perl itself
  11
  12 # POD documentation - main docs before the code
  13
  14 =head1 NAME
  15
  16 Bio::SearchDist - A perl wrapper around Sean Eddy's histogram object
  17
  18 =head1 SYNOPSIS
  19
  20   $dis = Bio::SearchDist->new();
  21   foreach $score ( @scores ) {
  22      $dis->add_score($score);
  23   }
  24
  25   if( $dis->fit_evd() ) {
  26     foreach $score ( @scores ) {
  27       $evalue = $dis->evalue($score);
  28       print "Score $score had an evalue of $evalue\n";
  29     }
  30   } else {
  31     warn("Could not fit histogram to an EVD!");
  32   }
  33
  34
  35 =head1 DESCRIPTION
  36
  37 The Bio::SearchDist object is a wrapper around Sean Eddy's excellent
  38 histogram object. The histogram object can bascially take in a number
  39 of scores which are sensibly distributed somewhere around 0 that come
  40 from a supposed Extreme Value Distribution. Having add all the scores
  41 from a database search via the add_score method you can then fit a
  42 extreme value distribution using ->fit_evd(). Once fitted you can then
  43 get out the evalue for each score (or a new score) using
  44 ->evalue($score).
  45
  46 The fitting procedure is better described in Sean Eddy's own code
  47 (available from http://hmmer.wustl.edu, or in the histogram.h header
  48 file in Compile/SW). Bascially it fits a EVD via a maximum likelhood
  49 method with pruning of the top end of the distribution so that real
  50 positives are discarded in the fitting procedure. This comes from
  51 an orginally idea of Richard Mott's and the likelhood fitting
  52 is from a book by Lawless [should ref here].
  53
  54
  55 The object relies on the fact that the scores are sensibly distributed
  56 around about 0 and that integer bins are sensible for the
  57 histogram. Scores based on bits are often ideal for this (bits based
  58 scoring mechanisms is what this histogram object was originally
  59 designed for).
  60
  61
  62 =head1 CONTACT
  63
  64 The original code this was based on comes from the histogram module as
  65 part of the HMMer2 package. Look at http://hmmer.wustl.edu/
  66
  67 Its use in Bioperl is via the Compiled XS extension which is cared for
  68 by Ewan Birney (birney@sanger.ac.uk). Please contact Ewan first about
  69 the use of this module
  70
  71
  72 =head2 Reporting Bugs
  73
  74 Report bugs to the Bioperl bug tracking system to help us keep track
  75 the bugs and their resolution. Bug reports can be submitted via email
  76 or the web:
  77
  78     bioperl-bugs@bio.perl.org
  79     http://bio.perl.org/bioperl-bugs/
  80
  81
  82 =head1 APPENDIX
  83
  84 The rest of the documentation details each of the object
  85 methods. Internal methods are usually preceded with a _
  86
  87 =cut
  88
  89
  90 # Let the code begin...
  91
  92
  93 package Bio::SearchDist;
  94 use vars qw($AUTOLOAD @ISA);
  95 use strict;
  96
  97 # Object preamble - inheriets from Bio::Root::Object
  98
  99 use Bio::Root::Object;
 100
 101 BEGIN {
 102     eval {
 103         require bp_sw;
 104     };
 105     if ( $@ ) {
 106         print STDERR ("\nThe C-compiled engine for histogram object (bp_sw) has not been installed.\n Please read the installation instructions for bioperl for using the compiled extensions\n\n");
 107         exit(1);
 108     }
 109 }
 110
 111
 112 @ISA = qw(Bio::Root::Object Exporter);
 113 @EXPORT_OK = qw();
 114 # new() is inherited from Bio::Root::Object
 115
 116 # _initialize is where the heavy stuff will happen when new is called
 117
 118 sub _initialize {
 119   my($self,%args) = @_;
 120   my($min, $max, $lump) =
 121         $self->_rearrange([qw(MIN MAX LUMP)], %args);
 122
 123   my $make = $self->SUPER::_initialize(%args);
 124
 125   if( !  $min ) {
 126     $min = -100;
 127   }
 128
 129   if( !  $max ) {
 130     $max = +100;
 131   }
 132
 133   if( ! $lump ) {
 134     $lump = 50;
 135   }
 136
 137   $self->_engine(&Wise2::new_Histogram($min,$max,$lump));
 138
 139   return $make; # success - we hope!
 140 }
 141
 142 =head2 add_score
 143
 144  Title   : add_score
 145  Usage   : $dis->add_score(300);
 146  Function: Adds a single score to the distribution
 147  Returns : nothing
 148  Args    :
 149
 150
 151 =cut
 152
 153 sub add_score{
 154    my ($self,$score) = @_;
 155    my ($eng);
 156    $eng = $self->_engine();
 157    $eng->AddToHistogram($score);
 158
 159 }
 160
 161 =head2 fit_evd
 162
 163  Title   : fit_evd
 164  Usage   : $dis->fit_evd();
 165  Function: fits an evd to the current distribution
 166  Returns : 1 if it fits successfully, 0 if not
 167  Args    :
 168
 169
 170 =cut
 171
 172 sub fit_evd{
 173    my ($self,@args) = @_;
 174
 175    return $self->_engine()->ExtremeValueFitHistogram(10000);
 176
 177 }
 178
 179 =head2 evalue
 180
 181  Title   : evalue
 182  Usage   : $eval = $dis->evalue($score)
 183  Function: Returns the evalue of this score
 184  Returns : float
 185  Args    :
 186
 187
 188 =cut
 189
 190 sub evalue{
 191    my ($self,$score) = @_;
 192
 193    return $self->_engine()->evalue($score);
 194
 195 }
 196
 197
 198
 199 =head2 _engine
 200
 201  Title   : _engine
 202  Usage   : $obj->_engine($newval)
 203  Function: underlyine bp_sw:: histogram engine
 204  Returns : value of _engine
 205  Args    : newvalue (optional)
 206
 207
 208 =cut
 209
 210 sub _engine{
 211    my ($self,$value) = @_;
 212    if( defined $value) {
 213       $obj->{'_engine'} = $value;
 214     }
 215     return $obj->{'_engine'};
 216
 217 }
 218
 219
 220 ## End of Package
 221
 222 1;
 223 __END__
 224