scripts/run_protdist.PLS

   1 #!/usr/bin/perl -w
   2 use strict;
   3 #!/usr/bin/perl -w
   4 use strict;
   5
   6 =head1 NAME
   7
   8 run_neighbor - run Phylip's 'protdist' program through Bioperl
   9
  10 =head1 SYNOPSIS
  11
  12 run_protdist [-i inputfile] [-o outfilename]
  13
  14 =head1 DESCRIPTION
  15
  16 Provide an alignment file to run protdist on.  File should be named
  17 either .aln or .phy.  This is required so that we can determine if we
  18 need to convert a clustalw alignment into phylip.  You are welcome to
  19 extend the script to work on other MSA formats which bioperl supports.
  20 This is intended to be used in very simple manual pipelines.
  21
  22 The input file should be named in the form of file.phy or file.aln
  23 the program expects a file in the form of (\S+)\.(\S+).
  24
  25 This will run the application 'protdist' using the 'KIMURA' formula to
  26 build a a protein distance matrix.  Those with phylip3.6 will want to
  27 make some changes if they want to use JTT.  I'm happy to help add this
  28 in as a cmd-line argument if it is requested.
  29
  30 =head1 AUTHOR
  31
  32 Jason Stajich, jason-AT-open-bio-DOT-org
  33
  34 =cut
  35
  36 use Bio::AlignIO;
  37 use Bio::Tools::Run::Phylo::Phylip::ProtDist;
  38 use Getopt::Long;
  39
  40 my @params = ( 'MODEL' => 'KIMURA',
  41                'quiet' => 1,
  42                );
  43
  44 my ($out,$file);
  45
  46 GetOptions(
  47            'o|out:s'   => \$out,
  48            'i|in:s'    =>  \$file,
  49            'h|help'    => sub { exec('perldoc',$0); exit(0) }
  50            );
  51
  52 ($file) ||= shift @ARGV;
  53
  54 my ($stem,$ext) = ($file =~ /(\S+)\.(\S+)$/);
  55 $stem ||= $file;
  56
  57 my $outfh;
  58 if( $out ) {
  59     open($outfh, ">$out") || die($!);
  60 } else {
  61     open($outfh, ">$stem.matrix") || die($!);
  62 }
  63
  64 if( $ext eq 'aln' ) {
  65     my $inaln = new Bio::AlignIO(-file => $file,
  66                                 -format => 'clustalw');
  67     $file = "$stem.phy";
  68
  69     my $outreformat = new Bio::AlignIO(-file => ">$file",
  70                                        -interleaved => 1,
  71                                        -format => 'phylip');
  72     while( my $aln = $inaln->next_aln ) {
  73         $outreformat->write_aln($aln);
  74     }
  75     $outreformat->close();
  76     $outreformat = undef;
  77     $inaln       = undef;
  78 }
  79
  80 my $factory = new Bio::Tools::Run::Phylo::Phylip::ProtDist(@params);
  81 my (@matrix) = $factory->create_distance_matrix($file);
  82
  83 foreach my $mat ( @matrix ) {
  84     print $outfh $mat->print_matrix;
  85 }
  86