Bio/Tools/Run/Phylo/Phylip/ProtDist.pm

   1 # BioPerl module for Bio::Tools::Run::Phylo::Phylip::ProtDist
   2 #
   3 # Created by
   4 #
   5 # Shawn Hoon
   6 #
   7 # You may distribute this module under the same terms as perl itself
   8
   9 # POD documentation - main docs before the code
  10
  11 =head1 NAME
  12
  13 Bio::Tools::Run::Phylo::Phylip::ProtDist - Wrapper for the phylip
  14 program protdist by Joseph Felsentein for creating a distance matrix
  15 comparing protein sequences from a multiple alignment file or a
  16 L<Bio::SimpleAlign> object and returns a hash ref to the table
  17
  18 =head1 SYNOPSIS
  19
  20   #Create a SimpleAlign object
  21   @params = ('ktuple' => 2, 'matrix' => 'BLOSUM');
  22   $factory = Bio::Tools::Run::Alignment::Clustalw->new(@params);
  23   $inputfilename = 't/data/cysprot.fa';
  24   $aln = $factory->align($inputfilename); # $aln is a SimpleAlign object.
  25
  26
  27   # Create the Distance Matrix using a default PAM matrix and id name
  28   # lengths limit of 30 note to use id name length greater than the
  29   # standard 10 in protdist, you will need to modify the protdist source
  30   # code
  31   @params = ('MODEL' => 'PAM');
  32   $protdist_factory = Bio::Tools::Run::Phylo::Phylip::ProtDist->new(@params);
  33   my $matrix  = $protdist_factory->create_distance_matrix($aln);
  34
  35   #finding the distance between two sequences
  36   my $distance = $matrix->{'protein_name_1'}{'protein_name_2'};
  37
  38   #Alternatively, one can create the matrix by passing in a file
  39   #name containing a multiple alignment in phylip format
  40   $protdist_factory = Bio::Tools::Run::Phylo::Phylip::ProtDist->new(@params);
  41   my $matrix  =
  42     $protdist_factory->create_distance_matrix('/home/shawnh/prot.phy');
  43
  44
  45 =head1 PARAMETERS FOR PROTDIST COMPUTATION
  46
  47 =head2 MODEL
  48
  49 Title           : MODEL
  50 Description     : (optional)
  51
  52                   This sets the model of amino acid substitution used
  53                   in the calculation of the distances.  3 different
  54                   models are supported:
  55                   PAM     Dayhoff PAM Matrix(default)
  56                   KIMURA  Kimura's Distance CAT
  57
  58                   Categories Distance Usage: @params =
  59                   ('model'=>'X');#where X is one of the values above
  60
  61                   Defaults to PAM For more information on the usage of
  62                   the different models, please refer to the
  63                   documentation
  64                   defaults to Equal
  65                   (0.25,0.25,0.25,0.25) found in the phylip package.
  66
  67 =head2 ALL SUBSEQUENT PARAMETERS WILL ONLY WORK IN CONJUNCTION WITH
  68 THE Categories Distance MODEL*
  69
  70 =cut
  71
  72 =head2 GENCODE
  73
  74   Title         : GENCODE
  75   Description   : (optional)
  76
  77                   This option allows the user to select among various
  78                   nuclear and mitochondrial genetic codes.
  79
  80                   Acceptable Values:
  81                   U           Universal
  82                   M           Mitochondrial
  83                   V           Vertebrate mitochondrial
  84                   F           Fly mitochondrial
  85                   Y           Yeast mitochondrial
  86   Usage         : @params = ('gencode'=>'X');
  87                   where X is one of the letters above
  88                   Defaults to U
  89
  90 =head2 CATEGORY
  91
  92 Title           : CATEGORY
  93 Description     : (optional)
  94
  95                   This option sets the categorization of amino acids
  96                   all have groups: (Glu Gln Asp Asn), (Lys Arg His),
  97                   (Phe Tyr Trp)  plus:
  98                   G   George/Hunt/Barker:
  99                           (Cys), (Met   Val  Leu  Ileu),
 100                           (Gly  Ala  Ser  Thr  Pro)
 101                   C   Chemical:
 102                           (Cys   Met), (Val  Leu  Ileu  Gly  Ala  Ser  Thr),
 103                           (Pro)
 104                   H   Hall:
 105                         (Cys), (Met   Val  Leu  Ileu), (Gly  Ala  Ser  Thr),
 106                         (Pro)
 107
 108   Usage         : @params = ('category'=>'X');
 109                   where X is one of the letters above
 110                   Defaults to G
 111
 112 =head2 PROBCHANGE
 113
 114   Title       : PROBCHANGE
 115   Description : (optional)
 116                  This option sets the ease of changing category of amino
 117                  acid.  (1.0 if no difficulty of changing,less if less
 118                  easy. Can't be negative)
 119
 120   Usage       : @params = ('probchange'=>X) where 0<=X<=1
 121                 Defaults to 0.4570
 122
 123 =head2 TRANS
 124
 125   Title       : TRANS
 126   Description : (optional)
 127                 This option sets transition/transversion ratio can be
 128                 any positive number
 129
 130   Usage        : @params = ('trans'=>X) where X >= 0
 131                  Defaults to 2
 132
 133 =head2 FREQ
 134
 135   Title       : FREQ
 136   Description : (optional)
 137                 This option sets the frequency of each base (A,C,G,T)
 138                 The sum of the frequency must sum to 1.
 139                 For example A,C,G,T = (0.25,0.5,0.125,0.125)
 140
 141   Usage       : @params = ('freq'=>('W','X','Y','Z')
 142                 where W + X + Y + Z = 1
 143                 Defaults to Equal (0.25,0.25,0.25,0.25)
 144
 145
 146 =head1 FEEDBACK
 147
 148 =head2 Mailing Lists
 149
 150 User feedback is an integral part of the evolution of this and other
 151 Bioperl modules. Send your comments and suggestions preferably to one
 152 of the Bioperl mailing lists.  Your participation is much appreciated.
 153
 154   bioperl-l@bioperl.org          - General discussion
 155   http://bio.perl.org/MailList.html             - About the mailing lists
 156
 157 =head2 Reporting Bugs
 158
 159 Report bugs to the Bioperl bug tracking system to help us keep track
 160  the bugs and their resolution.  Bug reports can be submitted via
 161  email or the web:
 162
 163   bioperl-bugs@bio.perl.org
 164   http://bio.perl.org/bioperl-bugs/
 165
 166 =head1 AUTHOR - Shawn Hoon
 167
 168 Email shawnh@fugu-sg.org
 169
 170 =head1 APPENDIX
 171
 172 The rest of the documentation details each of the object
 173 methods. Internal methods are usually preceded with a _
 174
 175 =cut
 176
 177 #'
 178
 179
 180 package Bio::Tools::Run::Phylo::Phylip::ProtDist;
 181
 182 use vars qw($AUTOLOAD @ISA $PROGRAM $PROGRAMDIR $PROGRAMNAME
 183             $TMPDIR $TMPOUTFILE @PROTDIST_PARAMS @OTHER_SWITCHES
 184             %OK_FIELD);
 185 use strict;
 186 use Bio::SimpleAlign;
 187 use Bio::AlignIO;
 188 use Bio::TreeIO;
 189 use Bio::Tools::Run::Phylo::Phylip::Base;
 190 use Cwd;
 191
 192 # inherit from Phylip::Base which has some methods for dealing with
 193 # Phylip specifics
 194 @ISA = qw(Bio::Tools::Run::Phylo::Phylip::Base);
 195
 196 # You will need to enable the protdist program. This
 197 # can be done in (at least) 3 ways:
 198 #
 199 # 1. define an environmental variable PHYLIPDIR:
 200 # export PHYLIPDIR=/home/shawnh/PHYLIP/bin
 201 #
 202 # 2. include a definition of an environmental variable CLUSTALDIR in
 203 # every script that will use Clustal.pm.
 204 # $ENV{PHYLIPDIR} = '/home/shawnh/PHYLIP/bin';
 205 #
 206 # 3. You can set the path to the program through doing:
 207 # my @params('program'=>'/usr/local/bin/protdist');
 208 # my $protdist_factory = Bio::Tools::Run::Phylo::Phylip::ProtDist->new(@params);
 209 #
 210
 211
 212 BEGIN {
 213     $PROGRAMNAME = 'protdist'  . ($^O =~ /mswin/i ?'.exe':'');
 214     if (defined $ENV{PHYLIPDIR}) {
 215         $PROGRAMDIR = $ENV{PHYLIPDIR} || '';
 216         $PROGRAM = Bio::Root::IO->catfile($PROGRAMDIR,
 217                                           'protdist'.($^O =~ /mswin/i ?'.exe':''));
 218     }
 219         @PROTDIST_PARAMS = qw(MODEL GENCODE CATEGORY PROBCHANGE TRANS FREQ);
 220         @OTHER_SWITCHES = qw(QUIET);
 221         foreach my $attr(@PROTDIST_PARAMS,@OTHER_SWITCHES) {
 222                 $OK_FIELD{$attr}++;
 223         }
 224 }
 225
 226 sub new {
 227     my ($class,@args) = @_;
 228     my $self = $class->SUPER::new(@args);
 229     # to facilitiate tempfile cleanup
 230     $self->io->_initialize_io();
 231
 232     my ($attr, $value);
 233     ($TMPDIR) = $self->io->tempdir(CLEANUP=>1);
 234     (undef,$TMPOUTFILE) = $self->io->tempfile(-dir => $TMPDIR);
 235     while (@args)  {
 236         $attr =   shift @args;
 237         $value =  shift @args;
 238         next if( $attr =~ /^-/ ); # don't want named parameters
 239         if ($attr =~/PROGRAM/i) {
 240             $self->executable($value);
 241             next;
 242         }
 243         if ($attr =~ /IDLENGTH/i){
 244             $self->idlength($value);
 245             next;
 246         }
 247         $self->$attr($value);
 248     }
 249     return $self;
 250 }
 251
 252 sub AUTOLOAD {
 253     my $self = shift;
 254     my $attr = $AUTOLOAD;
 255     $attr =~ s/.*:://;
 256     $attr = uc $attr;
 257     $self->throw("Unallowed parameter: $attr !") unless $OK_FIELD{$attr};
 258     $self->{$attr} = shift if @_;
 259     return $self->{$attr};
 260 }
 261
 262
 263 =head2 executable
 264
 265  Title   : executable
 266  Usage   : $obj->executable($newval)
 267  Function: Finds the full path to the 'protdist' executable
 268  Returns : string representing the full path to the exe
 269  Args    : [optional] name of executable to set path to
 270            [optional] boolean flag whether or not warn when exe is not found
 271
 272 =cut
 273
 274 sub executable{
 275    my ($self, $exe,$warn) = @_;
 276
 277    if( defined $exe ) {
 278      $self->{'_pathtoexe'} = $exe;
 279    }
 280
 281    unless( defined $self->{'_pathtoexe'} ) {
 282        if( $PROGRAM && -e $PROGRAM && -x $PROGRAM ) {
 283            $self->{'_pathtoexe'} = $PROGRAM;
 284        } else {
 285            my $exe;
 286            if( ( $exe = $self->io->exists_exe($PROGRAMNAME) ) &&
 287                -x $exe ) {
 288                $self->{'_pathtoexe'} = $exe;
 289            } else {
 290                $self->warn("Cannot find executable for $PROGRAMNAME") if $warn;
 291                $self->{'_pathtoexe'} = undef;
 292            }
 293        }
 294    }
 295    $self->{'_pathtoexe'};
 296 }
 297
 298 =head2 idlength
 299
 300  Title   : idlength
 301  Usage   : $obj->idlength ($newval)
 302  Function:
 303  Returns : value of idlength
 304  Args    : newvalue (optional)
 305
 306
 307 =cut
 308
 309 sub idlength{
 310    my $self = shift;
 311    if( @_ ) {
 312       my $value = shift;
 313       $self->{'idlength'} = $value;
 314     }
 315     return $self->{'idlength'};
 316
 317 }
 318
 319
 320 =head2  create_distance_matrix
 321
 322  Title   : create_distance_matrix
 323  Usage   :
 324         $inputfilename = 't/data/prot.phy';
 325         $matrix= $prodistfactory->create_distance_matrix($inputfilename);
 326 or
 327         $seq_array_ref = \@seq_array; @seq_array is array of Seq objs
 328         $aln = $protdistfactory->align($seq_array_ref);
 329         $matrix = $protdistfactory->create_distance_matrix($aln);
 330
 331  Function: Create a distance matrix from a SimpleAlign object or a multiple alignment file
 332  Example :
 333  Returns : Hash ref to a hash of a hash
 334  Args    : Name of a file containing a multiple alignment in Phylip format
 335            or an SimpleAlign object
 336
 337  Throws an exception if argument is not either a string (eg a
 338  filename) or a Bio::SimpleAlign object. If
 339  argument is string, throws exception if file corresponding to string
 340  name can not be found.
 341
 342 =cut
 343
 344 sub create_distance_matrix{
 345
 346     my ($self,$input) = @_;
 347     my ($infilename);
 348
 349 # Create input file pointer
 350         $infilename = $self->_setinput($input);
 351     if (!$infilename) {$self->throw("Problems setting up for protdist. Probably bad input data in $input !");}
 352
 353 # Create parameter string to pass to protdist program
 354     my $param_string = $self->_setparams();
 355 # run protdist
 356     my $aln = $self->_run($infilename,$param_string);
 357 }
 358
 359 #################################################
 360
 361 =head2  _run
 362
 363  Title   :  _run
 364  Usage   :  Internal function, not to be called directly
 365  Function:  makes actual system call to protdist program
 366  Example :
 367  Returns : Bio::Tree object
 368  Args    : Name of a file containing a set of multiple alignments in Phylip format
 369            and a parameter string to be passed to protdist
 370
 371
 372 =cut
 373
 374 sub _run {
 375     my ($self,$infile,$param_string) = @_;
 376     my $instring;
 377     $instring =  $infile."\n$param_string";
 378     $self->debug( "Program ".$self->executable." $param_string\n");
 379
 380         #open a pipe to run protdist to bypass interactive menus
 381     if ($self->quiet() || $self->verbose() < 0) {
 382         open(PROTDIST,"|".$self->executable.">/dev/null");
 383     }
 384     else {
 385         open(PROTDIST,"|".$self->executable);
 386     }
 387     print PROTDIST $instring;
 388     close(PROTDIST);
 389
 390         #get the results
 391     my $path = cwd;
 392     chomp($path);
 393     my $outfile = $self->io->catfile($path,$self->outfile);
 394
 395     $self->throw("protdist did not create matrix correctly ($outfile)")
 396         unless (-e $outfile);
 397
 398         #Create the distance matrix here
 399     my @values;
 400
 401     open(DIST, $outfile);
 402     while (<DIST>){
 403         next if (/^\s+\d+$/);
 404         my @line = split /\s+/,$_;
 405         push @values,[@line];
 406     }
 407     close(DIST);
 408         #list of sequences
 409     my @name = map{$_->[0]}@values;
 410     my %dist;
 411         #create the matrix using a hash of hash
 412     my $i = 0;
 413     foreach my $name (@name){
 414         my $j = 1;
 415         foreach my $n(@name){
 416                 $dist{$name}{$n} = $values[$i][$j];
 417                 $j++;
 418         }
 419         $i++;
 420     }
 421
 422     # Clean up the temporary files created along the way...
 423     unlink $outfile unless $self->save_tempfiles;
 424
 425     return \%dist;
 426 }
 427
 428
 429 =head2  _setinput()
 430
 431  Title   :  _setinput
 432  Usage   :  Internal function, not to be called directly
 433  Function:   Create input file for protdist program
 434  Example :
 435  Returns : name of file containing a multiple alignment in Phylip format
 436  Args    : SimpleAlign object reference or input file name
 437
 438
 439 =cut
 440
 441 sub _setinput {
 442     my ($self, $input) = @_;
 443     my ($alnfilename,$tfh);
 444
 445     # suffix is used to distinguish alignment files  from an align obkect
 446         #If $input is not a  reference it better be the name of a file with the sequence/
 447
 448     #  a phy formatted alignment file
 449         unless (ref $input) {
 450         # check that file exists or throw
 451         $alnfilename= $input;
 452         unless (-e $input) {return 0;}
 453                 return $alnfilename;
 454     }
 455
 456     #  $input may be a SimpleAlign Object
 457     if ($input->isa("Bio::SimpleAlign")) {
 458         #  Open temporary file for both reading & writing of BioSeq array
 459         ($tfh,$alnfilename) = $self->io->tempfile(-dir=>$TMPDIR);
 460         my $alnIO = Bio::AlignIO->new(-fh => $tfh,
 461                                       -format=>'phylip',
 462                                       -idlength=>$self->idlength());
 463         $alnIO->write_aln($input);
 464         $alnIO->close();
 465         close($tfh);
 466         return $alnfilename;
 467     }
 468     return 0;
 469 }
 470
 471 =head2  _setparams()
 472
 473  Title   :  _setparams
 474  Usage   :  Internal function, not to be called directly
 475  Function:   Create parameter inputs for protdist program
 476  Example :
 477  Returns : parameter string to be passed to protdist
 478  Args    : name of calling object
 479
 480 =cut
 481
 482 sub _setparams {
 483     my ($attr, $value, $self);
 484
 485     #do nothing for now
 486     $self = shift;
 487     my $param_string = "";
 488     my $cat = 0;
 489     foreach  my $attr ( @PROTDIST_PARAMS) {
 490         $value = $self->$attr();
 491         next unless (defined $value);
 492         if ($attr =~/MODEL/i){
 493             if ($value=~/CAT/i){
 494                 $cat = 1;
 495                 $param_string .= "P\nP\n";
 496                 next;
 497             }
 498             elsif($value=~/KIMURA/i){
 499                 $param_string .= "P\nY\n";
 500                 return $param_string;
 501             }
 502             else {
 503                 $param_string.="Y\n";
 504                 return $param_string;
 505             }
 506         }
 507         if ($cat == 1){
 508             if($attr =~ /GENCODE/i){
 509                 $self->throw("Unallowed value for genetic code") unless ($value =~ /[UMVFY]/);
 510                 $param_string .= "C\n$value\n";
 511             }
 512             if ($attr =~/CATEGORY/i){
 513                 $self->throw("Unallowed value for categorization of amino acids") unless ($value =~/[CHG]/);
 514                 $param_string .= "A\n$value\n";
 515             }
 516             if ($attr =~/PROBCHANGE/i){
 517                 if (($value =~ /\d+/)&&($value >= 0) && ($value < 1)){
 518                     $param_string .= "E\n$value\n";
 519                 }
 520                 else {
 521                     $self->throw("Unallowed value for probability change category");
 522                 }
 523             }
 524             if ($attr =~/TRANS/i){
 525                 if (($value=~/\d+/) && ($value >=0)){
 526                     $param_string .="T\n$value\n";
 527                 }
 528             }
 529             if ($attr =~ /FREQ/i){
 530                 my @freq = split(",",$value);
 531                 if ($freq[0] !~ /\d+/){ #a letter provided (sets frequencies equally to 0.25)
 532                     $param_string .="F\n".$freq[0]."\n";
 533                 }
 534                 elsif ($#freq ==  3) {#must have 4 digits for each base
 535                                           $param_string .="F\n";
 536                                           foreach my $f (@freq){
 537                                               $param_string.="$f\n";
 538                                           }
 539                                       }
 540                 else {
 541                     $self->throw("Unallowed value fo base frequencies");
 542                 }
 543             }
 544         }
 545     }
 546     $param_string .="Y\n";
 547
 548     return $param_string;
 549 }
 550
 551
 552
 553 =head1 Bio::Tools::Run::Wrapper methods
 554
 555 =cut
 556
 557 =head2 no_param_checks
 558
 559  Title   : no_param_checks
 560  Usage   : $obj->no_param_checks($newval)
 561  Function: Boolean flag as to whether or not we should
 562            trust the sanity checks for parameter values
 563  Returns : value of no_param_checks
 564  Args    : newvalue (optional)
 565
 566
 567 =cut
 568
 569 =head2 save_tempfiles
 570
 571  Title   : save_tempfiles
 572  Usage   : $obj->save_tempfiles($newval)
 573  Function:
 574  Returns : value of save_tempfiles
 575  Args    : newvalue (optional)
 576
 577
 578 =cut
 579
 580 =head2 outfile_name
 581
 582  Title   : outfile_name
 583  Usage   : my $outfile = $protdist->outfile_name();
 584  Function: Get/Set the name of the output file for this run
 585            (if you wanted to do something special)
 586  Returns : string
 587  Args    : [optional] string to set value to
 588
 589
 590 =cut
 591
 592
 593 =head2 tempdir
 594
 595  Title   : tempdir
 596  Usage   : my $tmpdir = $self->tempdir();
 597  Function: Retrieve a temporary directory name (which is created)
 598  Returns : string which is the name of the temporary directory
 599  Args    : none
 600
 601
 602 =cut
 603
 604 =head2 cleanup
 605
 606  Title   : cleanup
 607  Usage   : $codeml->cleanup();
 608  Function: Will cleanup the tempdir directory after a PAML run
 609  Returns : none
 610  Args    : none
 611
 612
 613 =cut
 614
 615 =head2 io
 616
 617  Title   : io
 618  Usage   : $obj->io($newval)
 619  Function:  Gets a L<Bio::Root::IO> object
 620  Returns : L<Bio::Root::IO>
 621  Args    : none
 622
 623
 624 =cut
 625
 626 1; # Needed to keep compiler happy