[Bio] / FigKernelScripts / make_proml_tree.pl Repository:
ViewVC logotype

Diff of /FigKernelScripts/make_proml_tree.pl

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.5, Fri Jan 12 18:44:25 2007 UTC revision 1.6, Thu Feb 1 20:49:10 2007 UTC
# Line 1  Line 1 
1  ########################################################################  #! /usr/bin/perl
2    #
3    #   make_proml_tree_2:  A command line interface to the proml program
4    #
5    my $usage = <<"End_of_Usage";
6    
7    Usage:
8    
9       make_proml_tree_2 < Alignment > Trees_with_likelihood_comment
10    
11    options:
12      -a Alpha         Alpha parameter for gamma distribution of rates
13      -b GammaBins     Number of rate bins in gamma approximation (2 - 9)
14      -c RatesCatFile  Rates and categories file (rates \\n categories)
15      -d Directory     Directory for temp files. If it exists, files are retained.
16      -g               Do global rearrangements
17      -h HMMFile       User-defined rate HMM file [Not functional]
18      -i FracInvar     Fraction of invariant sites
19      -j JumbleSeed    Random addition order seed (odd)
20      -l               User-defined branch lengths on user trees
21      -m Model         Amino acid change model (JTT (D) | PMB | PAM)
22      -n NJumbles      Number of random addition orders
23      -p Persistance   Number of sites with correlated rates
24      -s               Slower, more accurate calculation
25      -t Tmp           Directory for directory of temporary files (D = /tmp)
26      -u UserTreeFile  File of user trees
27      -v CoefOfVar     Coeficient of variation for gamma distribution for rates
28      -w Weights       Site weights file
29    
30    End_of_Usage
31    
32    
33    use proml;
34    use gjoseqlib;
35    use gjonewicklib;
36    
37  use FIG;  my %options = ( tree_format  => 'gjo' );
38  use tree_utilities;  
39  use Sys::Hostname;  while ( @ARGV && $ARGV[0] =~ /^-/ )
40    {
41  $usage = "usage: make_proml_tree [-a Alpha] Alignment Outfile OutTree";      my $flag = shift @ARGV;
42        my ( $file, $value, @trees );
43  my $host = hostname;  
44  my $dir = "$FIG_Config::temp/proml.$host.$$";      if    ( $flag =~ s/^-a// )
45  mkdir($dir,0777) || die "could not make $dir";      {
46            $value = $flag || shift @ARGV;
47  my $alpha = undef;          $value > 0 or usage( "Bad value of alpha parameter: $value\n" );
48  while ($ARGV[0] =~ /^-/)          $options{ alpha } = $value;
49  {      }
50      $_ = shift @ARGV;      elsif ( $flag =~ s/^-b// )
51      if ($_ =~ s/^-a//) { $alpha = ($_ || shift @ARGV) }      {
52      else               { die "Bad Flag: $_" }          $value = $flag || shift @ARGV;
53            $value > 0 && $value <= 9
54                or usage( "Bad value for gamma bins parameter: $value" );
55            $options{ gamma_bins } = $value;
56        }
57        elsif ( $flag =~ s/^-c// )
58        {
59            $file = $flag || shift @ARGV;
60            -f $file
61                and ( $value = read_categories_file( $file ) )
62                or  usage( "Bad category file: $file" );
63            $options{ categories } = $value;
64        }
65        elsif ( $flag =~ s/^-d// )
66        {
67            $value = $flag || shift @ARGV;
68            $value
69                or usage( "Bad value for directory for temp files: $value" );
70            $options{ tmp_dir } = $value;
71        }
72        elsif ( $flag =~ s/^-g// )
73        {
74            $options{ global } = 1;
75        }
76        elsif ( $flag =~ s/^-h// )
77        {
78            $file = $flag || shift @ARGV;
79            -f $file
80                or  usage( "Bad HMM file: $file" );
81        }
82        elsif ( $flag =~ s/^-i// )
83        {
84            $value = $flag || shift @ARGV;
85            $value >= 0 && $value < 1
86                or usage( "Bad value for fraction of invariant sites: $value" );
87            $options{ invar_frac } = $value;
88        }
89        elsif ( $flag =~ s/^-j// )
90        {
91            $value = $flag || shift @ARGV;
92            $value > 0 && ( $value % 2 == 1 )
93                or usage( "Bad value for jumble seed: $value" );
94            $options{ jumble_seed } = $value;
95        }
96        elsif ( $flag =~ s/^-l// )
97        {
98            $options{ user_lengths } = 1;
99        }
100        elsif ( $flag =~ s/^-m// )
101        {
102            $value = $flag || shift @ARGV;
103            $value
104                or usage( "Bad value for evolution model: $value" );
105            $options{ n_jumble } = $value;
106        }
107        elsif ( $flag =~ s/^-n// )
108        {
109            $value = $flag || shift @ARGV;
110            $value > 0
111                or usage( "Bad value for number of random addition orders: $value" );
112            $options{ n_jumble } = $value;
113        }
114        elsif ( $flag =~ s/^-p// )
115        {
116            $value = $flag || shift @ARGV;
117            $value > 0
118                or usage( "Bad value for site rate correlation persistance length: $value" );
119            $options{ persistance } = $value;
120        }
121        elsif ( $flag =~ s/^-s// )
122        {
123            $options{ slow } = 1;
124        }
125        elsif ( $flag =~ s/^-t// )
126        {
127            $value = $flag || shift @ARGV;
128            -d $value
129                or usage( "Bad value for tmp directory: $value" );
130            $options{ tmp } = $value;
131        }
132        elsif ( $flag =~ s/^-u// )
133        {
134            $file = $flag || shift @ARGV;
135            -f $file
136                and ( @trees = gjonewicklib::read_newick_trees( $file ) )
137                or  usage( "Bad tree file: $file" );
138            $options{ user_trees } = \@trees;
139        }
140        elsif ( $flag =~ s/^-v// )
141        {
142            $value = $flag || shift @ARGV;
143            $value >= 0
144                or usage( "Bad value for coeficient of variation of gamma distribution: $value" );
145            $options{ coef_of_var } = $value;
146        }
147        elsif ( $flag =~ s/^-w// )
148        {
149            my $file = $flag || shift @ARGV;
150            -f $file
151                and ( $value = read_weights_file( $file ) )
152                or  usage( "Bad weights file: $file" );
153            $options{ weights } = $value;
154        }
155        else
156        {
157            usage( "Bad flag: $flag" );
158        }
159    }
160    
161    my @align = gjoseqlib::read_fasta();
162    
163    my @results = proml::proml( \@align, \%options );
164    
165    foreach ( @results )
166    {
167        my ( $tree, $likelihood ) = @$_;
168        printf '[%12.4f]  ', $likelihood;
169        gjonewicklib::writeNewickTree( $tree );
170    }
171    
172    exit;
173    
174    
175    sub usage
176    {
177        foreach ( @_ ) { print STDERR "$_\n" }
178        print STDERR $usage;
179        exit;
180  }  }
181    
182    
183  (  sub read_categories_file
184   ($aliF    = shift @ARGV) &&  {
185   ($outF    = shift @ARGV) &&      open( CATS, "<$_[0]" ) or return undef;
186   ($treeF   = shift @ARGV)      $_ = <CATS>;
187  )      chomp;
188      || die $usage;      s/[\t ,]+/ /g;
189        my @rates = split;
190  &FIG::run("fasta_to_phylip < $aliF > $dir/infile");      my $categories = join( '', map { chomp; s/\D+//g; $_ } <CATS> );
191  open(PARMS,">$dir/parms") || die "could not open $dir/parms";      close( CATS );
192  $val = (int(rand() * 500000000) * 2) + 1;      [ \@rates, $categories ]
193  print PARMS "J\n$val\n1\n";  }
194  if    ($alpha) { print PARMS "R\nY\n",1 / sqrt($alpha),"\n5\n"; }  
195  else           { print PARMS "Y\n"; }  
196  close(PARMS);  sub read_weights_file
197  system "cd $dir; proml < parms > log";  {
198  system "mv $dir/outfile $outF; mv $dir/outtree $treeF; rm -r $dir";      open( WGTS, "<$_[0]" ) or return undef;
199  system "/bin/rm -r $dir";      my $weights = join( '', map { chomp; s/\D+//g; $_ } <WGTS> );
200        close( WGTS );
201        $weights
202    }
203    

Legend:
Removed from v.1.5  
changed lines
  Added in v.1.6

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3