← Back to branch summary

~ubuntu-branches/ubuntu/oneiric/bioperl/oneiric

« back to all changes in this revision

Viewing changes to Bio/Align/ProteinStatistics.pm

  • Committer: Bazaar Package Importer
  • Author(s): Charles Plessy
  • Date: 2007-09-21 22:52:22 UTC
  • mfrom: (1.2.2 upstream)
  • Revision ID: james.westby@ubuntu.com-20070921225222-tt20m2yy6ycuy2d8
Tags: 1.5.2.102-1
* Developer release.
* Upgraded source package to debhelper 5 and standards-version 3.7.2.
* Added libmodule-build-perl and libtest-harness-perl to
  build-depends-indep.
* Disabled automatic CRAN download.
* Using quilt instead of .diff.gz to manage modifications.
* Updated Recommends list for the binary package.
* Moved the "production-quality" scripts to /usr/bin/.
* New maintainer: Debian-Med packaging team mailing list.
* New uploaders: Charles Plessy and Steffen Moeller.
* Updated Depends, Recommends and Suggests.
* Imported in Debian-Med's SVN repository on Alioth.
* Executing the regression tests during package building.
* Moved the Homepage: field out from the package's description.
* Updated watch file.

Show diffs side-by-side

added added

removed removed

Lines of Context:
Bio/Align/ProteinStatistics.pm
 
1
# $Id: ProteinStatistics.pm,v 1.7.4.1 2006/10/02 23:10:12 sendu Exp $
 
2
#
 
3
# BioPerl module for Bio::Align::ProteinStatistics
 
4
#
 
5
# Cared for by Jason Stajich <jason-at-bioperl.org>
 
6
#
 
7
# Copyright Jason Stajich
 
8
#
 
9
# You may distribute this module under the same terms as perl itself
 
10
 
 
11
# POD documentation - main docs before the code
 
12
 
 
13
=head1 NAME
 
14
 
 
15
Bio::Align::ProteinStatistics - Calculate Protein Alignment statistics (mostly distances)
 
16
 
 
17
=head1 SYNOPSIS
 
18
 
 
19
  use Bio::Align::ProteinStatistics;
 
20
  use Bio::AlignIO;
 
21
  my $in = new Bio::AlignIO(-format => 'fasta',
 
22
                            -file   => 'pep-104.fasaln');
 
23
  my $aln = $in->next_aln;
 
24
 
 
25
  my $pepstats = Bio::Align::ProteinStatistics->new();
 
26
  $kimura = $protstats->distance(-align => $aln,
 
27
                                 -method => 'Kimura');
 
28
  print $kimura->print_matrix;
 
29
 
 
30
 
 
31
=head1 DESCRIPTION
 
32
 
 
33
This object is for generating various statistics from a protein
 
34
alignment.  Mostly it is where pairwise protein distances can be
 
35
calculated.
 
36
 
 
37
=head1 REFERENCES 
 
38
 
 
39
D_Kimura - Kimura, M. 1983. The Neutral Theory of Molecular Evolution. CUP, 
 
40
           Cambridge.
 
41
 
 
42
=head1 FEEDBACK
 
43
 
 
44
=head2 Mailing Lists
 
45
 
 
46
User feedback is an integral part of the evolution of this and other
 
47
Bioperl modules. Send your comments and suggestions preferably to
 
48
the Bioperl mailing list.  Your participation is much appreciated.
 
49
 
 
50
  bioperl-l@bioperl.org                  - General discussion
 
51
  http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
 
52
 
 
53
=head2 Reporting Bugs
 
54
 
 
55
Report bugs to the Bioperl bug tracking system to help us keep track
 
56
of the bugs and their resolution. Bug reports can be submitted via the
 
57
web:
 
58
 
 
59
  http://bugzilla.open-bio.org/
 
60
 
 
61
=head1 AUTHOR - Jason Stajich
 
62
 
 
63
Email jason-at-bioperl.org
 
64
 
 
65
=head1 APPENDIX
 
66
 
 
67
The rest of the documentation details each of the object methods.
 
68
Internal methods are usually preceded with a _
 
69
 
 
70
=cut
 
71
 
 
72
 
 
73
# Let the code begin...
 
74
 
 
75
 
 
76
package Bio::Align::ProteinStatistics;
 
77
use vars qw(%DistanceMethods $Precision $DefaultGapPenalty);
 
78
use strict;
 
79
 
 
80
use Bio::Align::PairwiseStatistics;
 
81
use Bio::Matrix::PhylipDist;
 
82
 
 
83
%DistanceMethods = ('kimura|k' => 'Kimura',
 
84
                    );
 
85
$Precision = 5;
 
86
$DefaultGapPenalty = 0;
 
87
 
 
88
use base qw(Bio::Root::Root Bio::Align::StatisticsI);
 
89
 
 
90
=head2 new
 
91
 
 
92
 Title   : new
 
93
 Usage   : my $obj = new Bio::Align::ProteinStatistics();
 
94
 Function: Builds a new Bio::Align::ProteinStatistics object 
 
95
 Returns : an instance of Bio::Align::ProteinStatistics
 
96
 Args    :
 
97
 
 
98
 
 
99
=cut
 
100
 
 
101
sub new {
 
102
  my($class,@args) = @_;
 
103
 
 
104
  my $self = $class->SUPER::new(@args);
 
105
  $self->pairwise_stats( new Bio::Align::PairwiseStatistics());
 
106
 
 
107
  return $self;
 
108
}
 
109
 
 
110
=head2 distance
 
111
 
 
112
 Title   : distance
 
113
 Usage   : my $distance_mat = $stats->distance(-align  => $aln, 
 
114
                                               -method => $method);
 
115
 Function: Calculates a distance matrix for all pairwise distances of
 
116
           sequences in an alignment.
 
117
 Returns : L<Bio::Matrix::PhylipDist> object
 
118
 Args    : -align  => Bio::Align::AlignI object
 
119
           -method => String specifying specific distance method 
 
120
                      (implementing class may assume a default)
 
121
 
 
122
=cut
 
123
 
 
124
sub distance{
 
125
   my ($self,@args) = @_;
 
126
   my ($aln,$method) = $self->_rearrange([qw(ALIGN METHOD)],@args);
 
127
   if( ! defined $aln || ! ref ($aln) || ! $aln->isa('Bio::Align::AlignI') ) { 
 
128
       $self->throw("Must supply a valid Bio::Align::AlignI for the -align parameter in distance");
 
129
   }
 
130
   $method ||= 'Kimura';
 
131
   foreach my $m ( keys %DistanceMethods ) {
 
132
       if(defined $m &&  $method =~ /$m/i ) {
 
133
           my $mtd = "D_$DistanceMethods{$m}";
 
134
           return $self->$mtd($aln);
 
135
       }
 
136
   }
 
137
   $self->warn("Unrecognized distance method $method must be one of [".
 
138
               join(',',$self->available_distance_methods())."]");
 
139
   return;
 
140
}
 
141
 
 
142
=head2 available_distance_methods
 
143
 
 
144
 Title   : available_distance_methods
 
145
 Usage   : my @methods = $stats->available_distance_methods();
 
146
 Function: Enumerates the possible distance methods
 
147
 Returns : Array of strings
 
148
 Args    : none
 
149
 
 
150
 
 
151
=cut
 
152
 
 
153
sub available_distance_methods{
 
154
   my ($self,@args) = @_;
 
155
   return values %DistanceMethods;
 
156
}
 
157
 
 
158
=head2 D - distance methods
 
159
 
 
160
 
 
161
=cut
 
162
 
 
163
 
 
164
=head2 D_Kimura
 
165
 
 
166
 Title   : D_Kimura
 
167
 Usage   : my $matrix = $pepstats->D_Kimura($aln);
 
168
 Function: Calculate Kimura protein distance (Kimura 1983) which 
 
169
           approximates PAM distance
 
170
           D = -ln ( 1 - p - 0.2 * p^2 )
 
171
 Returns : L<Bio::Matrix::PhylipDist>
 
172
 Args    : L<Bio::Align::AlignI>
 
173
 
 
174
 
 
175
=cut
 
176
 
 
177
# Kimura, M. 1983. The Neutral Theory of Molecular Evolution. CUP, Cambridge.
 
178
 
 
179
sub D_Kimura{
 
180
   my ($self,$aln) = @_;
 
181
   return 0 unless $self->_check_arg($aln);
 
182
   # ambiguities ignored at this point
 
183
   my (@seqs,@names,@values,%dist);
 
184
   my $seqct = 0;
 
185
   foreach my $seq ( $aln->each_seq) {
 
186
       push @names, $seq->display_id;
 
187
       push @seqs, uc $seq->seq();
 
188
       $seqct++;
 
189
   }
 
190
   my $len = $aln->length;
 
191
   my $precisionstr = "%.$Precision"."f";
 
192
 
 
193
   for( my $i = 0; $i < $seqct-1; $i++ ) {
 
194
       # (diagonals) distance is 0 for same sequence
 
195
       $dist{$names[$i]}->{$names[$i]} = [$i,$i];
 
196
       $values[$i][$i] = sprintf($precisionstr,0);
 
197
       for( my $j = $i+1; $j < $seqct; $j++ ) {
 
198
           my ($scored,$match) = (0,0);
 
199
           for( my $k=0; $k < $len; $k++ ) {
 
200
               my $m1 = substr($seqs[$i],$k,1);
 
201
               my $m2 = substr($seqs[$j],$k,1);
 
202
               if( $m1 ne '-' && $m2 ne '-' ) {
 
203
                   # score is number of scored bases (alignable bases)
 
204
                   # it could have also come from 
 
205
                   # my $L = $self->pairwise_stats->number_of_comparable_bases($pairwise);
 
206
                   # match is number of matches weighting ambiguity bases
 
207
                   # as well
 
208
                   $match += _check_ambiguity_protein($m1,$m2);
 
209
                   $scored++;
 
210
               }
 
211
           }
 
212
           # From Felsenstein's PHYLIP documentation:
 
213
           # This is very quick to do but has some obvious
 
214
           # limitations. It does not take into account which amino
 
215
           # acids differ or to what amino acids they change, so some
 
216
           # information is lost. The units of the distance measure
 
217
           # are fraction of amino acids differing, as also in the
 
218
           # case of the PAM distance. If the fraction of amino acids
 
219
           # differing gets larger than 0.8541 the distance becomes
 
220
           # infinite.
 
221
 
 
222
           my $D = 1 - ( $match / $scored );
 
223
           if( $D < 0.8541 ) {
 
224
               $D = - log ( 1 - $D - (0.2 * ($D ** 2)));
 
225
               $values[$j][$i] = $values[$i][$j] = sprintf($precisionstr,$D);
 
226
           } else { 
 
227
               $values[$j][$i] = $values[$i][$j] = '    NaN';
 
228
           }
 
229
           # fwd and rev lookup
 
230
           $dist{$names[$i]}->{$names[$j]} = [$i,$j];
 
231
           $dist{$names[$j]}->{$names[$i]} = [$i,$j];      
 
232
 
 
233
           # (diagonals) distance is 0 for same sequence
 
234
           $dist{$names[$j]}->{$names[$j]} = [$j,$j];      
 
235
           $values[$j][$j] = sprintf($precisionstr,0); 
 
236
 
 
237
       }
 
238
   }
 
239
   return Bio::Matrix::PhylipDist->new(-program => 'bioperl_PEPstats',
 
240
                                       -matrix  => \%dist,
 
241
                                       -names   => \@names,
 
242
                                       -values  => \@values); 
 
243
   
 
244
}
 
245
 
 
246
# some methods from EMBOSS distmat
 
247
sub _check_ambiguity_protein
 
248
{
 
249
    my ($t1,$t2) = @_;
 
250
    my $n = 0;
 
251
 
 
252
    if( $t1 ne 'X' && $t1 eq $t2 ) { 
 
253
        $n = 1.0;
 
254
    } elsif(  ((($t1 eq 'B' && $t2 eq 'DN') ||
 
255
               ($t2 eq 'B' && $t2 eq 'DN'))) ||
 
256
              
 
257
              ( ($t1 eq 'Z' && $t2 eq 'EQ') ||
 
258
                ($t2 eq 'Z' && $t1 eq 'EQ'))) {
 
259
        $n = 0.5;
 
260
    } elsif ( $t1 eq 'X' && $t2 eq 'X' ) {
 
261
        $n = 0.0025;
 
262
    } elsif(  $t1 eq 'X' || $t2 eq 'X' ) {
 
263
        $n = 0.05;
 
264
    }
 
265
    return $n;
 
266
}
 
267
 
 
268
=head2 Data Methods
 
269
 
 
270
=cut
 
271
 
 
272
=head2 pairwise_stats
 
273
 
 
274
 Title   : pairwise_stats
 
275
 Usage   : $obj->pairwise_stats($newval)
 
276
 Function: 
 
277
 Returns : value of pairwise_stats
 
278
 Args    : newvalue (optional)
 
279
 
 
280
 
 
281
=cut
 
282
 
 
283
sub pairwise_stats{
 
284
   my ($self,$value) = @_;
 
285
   if( defined $value) {
 
286
      $self->{'_pairwise_stats'} = $value;
 
287
    }
 
288
    return $self->{'_pairwise_stats'};
 
289
 
 
290
}
 
291
 
 
292
sub _check_arg {
 
293
    my($self,$aln ) = @_;
 
294
    if( ! defined $aln || ! $aln->isa('Bio::Align::AlignI') ) {
 
295
        $self->warn("Must provide a Bio::Align::AlignI compliant object to Bio::Align::DNAStatistics");
 
296
        return 0;
 
297
    } elsif( $aln->get_seq_by_pos(1)->alphabet ne 'protein' ) { 
 
298
        $self->warn("Must provide a protein alignment to Bio::Align::ProteinStatistics, you provided a " . $aln->get_seq_by_pos(1)->alphabet);
 
299
        return 0;
 
300
    }
 
301
    return 1;
 
302
}
 
303
 
 
304
1;