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

Annotation of /FigKernelScripts/make_proml_tree.pl

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.6 - (view) (download) (as text)

1 : golsen 1.6 #! /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 : overbeek 1.1
37 : golsen 1.6 my %options = ( tree_format => 'gjo' );
38 :    
39 :     while ( @ARGV && $ARGV[0] =~ /^-/ )
40 :     {
41 :     my $flag = shift @ARGV;
42 :     my ( $file, $value, @trees );
43 :    
44 :     if ( $flag =~ s/^-a// )
45 :     {
46 :     $value = $flag || shift @ARGV;
47 :     $value > 0 or usage( "Bad value of alpha parameter: $value\n" );
48 :     $options{ alpha } = $value;
49 :     }
50 :     elsif ( $flag =~ s/^-b// )
51 :     {
52 :     $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 : overbeek 1.1 {
167 : golsen 1.6 my ( $tree, $likelihood ) = @$_;
168 :     printf '[%12.4f] ', $likelihood;
169 :     gjonewicklib::writeNewickTree( $tree );
170 : overbeek 1.1 }
171 :    
172 : golsen 1.6 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 :     {
185 :     open( CATS, "<$_[0]" ) or return undef;
186 :     $_ = <CATS>;
187 :     chomp;
188 :     s/[\t ,]+/ /g;
189 :     my @rates = split;
190 :     my $categories = join( '', map { chomp; s/\D+//g; $_ } <CATS> );
191 :     close( CATS );
192 :     [ \@rates, $categories ]
193 :     }
194 :    
195 :    
196 :     sub read_weights_file
197 :     {
198 :     open( WGTS, "<$_[0]" ) or return undef;
199 :     my $weights = join( '', map { chomp; s/\D+//g; $_ } <WGTS> );
200 :     close( WGTS );
201 :     $weights
202 :     }
203 : overbeek 1.1

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3