[Bio] / FigKernelPackages / protdist_neighbor.pm Repository:
ViewVC logotype

Annotation of /FigKernelPackages/protdist_neighbor.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : golsen 1.1 package protdist_neighbor;
2 :    
3 :     #===============================================================================
4 :     # A perl interface to the protdist and neighbor programs in the PHYLIP
5 :     # program package
6 :     #
7 :     # $tree = protdist_neighbor( \@alignment, \%options )
8 :     # $tree = protdist_neighbor( \@alignment, %options )
9 :     # $tree = protdist_neighbor( \%options ) # alignment must be included as option
10 :     # $tree = protdist_neighbor( %options ) # alignment must be included as option
11 :     #
12 :     # [ [ dist11, dist12, ... ], [ dist21, dist22, ... ], ... ] = protdist( \@alignment, \%options )
13 :     # [ [ dist11, dist12, ... ], [ dist21, dist22, ... ], ... ] = protdist( \@alignment, %options )
14 :     # [ [ dist11, dist12, ... ], [ dist21, dist22, ... ], ... ] = protdist( \%options )
15 :     # [ [ dist11, dist12, ... ], [ dist21, dist22, ... ], ... ] = protdist( %options )
16 :     #
17 :     # @alignment = array of id_seq pairs, or id_definition_seq triples
18 :     #
19 :     #===============================================================================
20 :     #
21 :     # options:
22 :     #
23 :     # For protdist:
24 :     # alignment => \@alignment the way to supply the alignment as an option, rather than first param
25 :     # alpha => float alpha parameter of gamma distribution (0.5 - inf)
26 :     # categories => [ [ rates ], site_categories ]
27 :     # coef_of_var => float 1/sqrt(alpha) for gamma distribution (D = 0)
28 :     # invar_frac => 0 - 1 fraction of site that are invariant
29 :     # model => model evolution model JTT (D) | PMB | PAM
30 :     # persistance => float persistance length of rate category
31 :     # rate_hmm => [ [ rates ], [ probabilies ] ]
32 :     # weights => site_weights
33 :     #
34 :     # For neighbor (not really):
35 :     # jumble_seed => odd int jumble random seed
36 :     #
37 :     # Other:
38 :     # keep_duplicates => bool do not remove duplicate sequences (D = false) [NOT IMPLIMENTED]
39 :     # protdist => protdist allows fully defined path
40 :     # neighbor => neighbor allows fully defined path
41 :     # tmp => directory directory for tmp_dir (D = /tmp or .)
42 :     # tmp_dir => directory directory for temporary files (D = $tmp/protdist_neighbor.$$)
43 :     # tree_format => overbeek | gjo | fig format of output tree
44 :     #
45 :     # tmp_dir is created and deleted unless its name is supplied, and it already
46 :     # exists.
47 :     #
48 :     #===============================================================================
49 :    
50 :    
51 :     use strict;
52 :     use gjonewicklib qw( gjonewick_to_overbeek
53 :     newick_is_unrooted
54 :     newick_relabel_nodes
55 :     newick_tree_length
56 :     overbeek_to_gjonewick
57 :     parse_newick_tree_str
58 :     strNewickTree
59 :     uproot_newick
60 :     );
61 :    
62 :    
63 :     sub protdist_neighbor
64 :     {
65 :     my $align;
66 :     if ( ref( $_[0] ) eq 'ARRAY' )
67 :     {
68 :     $align = shift @_;
69 :     ( $align && ( ref( $align ) eq 'ARRAY' ) )
70 :     || ( ( print STDERR "protdist_neighbor::protdist_neighbor() called without alignment\n" )
71 :     && ( return () )
72 :     );
73 :     }
74 :    
75 :     my %options;
76 :     if ( $_[0] )
77 :     {
78 :     %options = ( ref( $_[0]) eq 'HASH' ) ? %{ $_[0] } : @_;
79 :     }
80 :    
81 :     #---------------------------------------------------------------------------
82 :     # Work on a copy of the alignment. Id is always first, seq is always last
83 :     #---------------------------------------------------------------------------
84 :    
85 :     $align ||= $options{ alignment } || $options{ align };
86 :    
87 :     my ( $seq, $id );
88 :     my %id;
89 :     my %local_id;
90 :     my $local_id = 'seq0000000';
91 :     my @align = map { $id = $_->[0];
92 :     $local_id++;
93 :     $id{ $local_id } = $id;
94 :     $local_id{ $id } = $local_id;
95 :     $seq = $_->[-1];
96 :     $seq =~ s/[BJOUZ]/X/gi; # Bad letters go to X
97 :     $seq =~ s/[^A-Z]/-/gi; # Anything else becomes -
98 :     [ $local_id, $seq ]
99 :     } @$align;
100 :    
101 :     #---------------------------------------------------------------------------
102 :     # Process protdist_neighbor options:
103 :     #---------------------------------------------------------------------------
104 :    
105 :     my $categories = $options{ categories }; # [ [ cat_rates ], site_cats ]
106 :     if ( $categories )
107 :     {
108 :     if ( ref( $categories ) ne 'ARRAY'
109 :     || ! ( ( @$categories == 2 ) || ( ( @$categories == 3 ) && ( shift @$categories ) ) )
110 :     || ref( $categories->[0] ) ne 'ARRAY'
111 :     )
112 :     {
113 :     print STDERR "proml::proml categories option value must be [ [ cat_rate1, ... ], site_categories ]\n";
114 :     return ();
115 :     }
116 :    
117 :     # Rate values cannot have very many decimal places or proml can't read it:
118 :    
119 :     @{$categories->[0]} = map { sprintf "%.6f", $_ } @{$categories->[0]};
120 :     }
121 :    
122 :     my $coef_of_var = $options{ coef_of_var }
123 :     || ( $options{ alpha } && ( $options{ alpha } > 0) && ( 1 / sqrt( $options{ alpha } ) ) )
124 :     || 0;
125 :     if ( $coef_of_var < 0 )
126 :     {
127 :     print STDERR "protdist_neighbor::protdist_neighbor coef_of_var option value must be >= 0\n";
128 :     return ();
129 :     }
130 :    
131 :     my $invar_frac = $options{ invar_frac } || 0;
132 :     if ( $invar_frac && ( $invar_frac < 0 || $invar_frac >= 1 ) )
133 :     {
134 :     print STDERR "protdist_neighbor::protdist_neighbor invar_frac option value must be >= 0 and < 1\n";
135 :     return ();
136 :     }
137 :    
138 :     my $jumble_seed = int( $options{ jumble_seed } ) || 0;
139 :     if ( $jumble_seed && ( ( $jumble_seed < 0) || ( $jumble_seed % 2 != 1 ) ) )
140 :     {
141 :     print STDERR "protdist_neighbor::protdist_neighbor jumble_seed option value must be an odd number > 0\n";
142 :     return ();
143 :     }
144 :    
145 :     my $model = ( $options{ model } =~ m/PAM/i ) ? 'PAM'
146 :     : ( $options{ model } =~ m/Dayhoff/i ) ? 'PAM'
147 :     : ( $options{ model } =~ m/PMB/i ) ? 'PMB'
148 :     : ( $options{ model } =~ m/Henikoff/i ) ? 'PMB'
149 :     : ( $options{ model } =~ m/Tillier/i ) ? 'PMB'
150 :     : ( $options{ model } =~ m/JTT/i ) ? 'JTT'
151 :     : ( $options{ model } =~ m/Jones/i ) ? 'JTT'
152 :     : ( $options{ model } =~ m/Taylor/i ) ? 'JTT'
153 :     : ( $options{ model } =~ m/Thornton/i ) ? 'JTT'
154 :     : ( $options{ model } =~ m/Kimura/i ) ? 'Kimura'
155 :     : 'JTT';
156 :    
157 :     my $persistance = $options{ persistance } || 0;
158 :     if ( $persistance && ( $persistance <= 1 ) )
159 :     {
160 :     print STDERR "protdist_neighbor::protdist_neighbor persistance option value must be > 1\n";
161 :     return ();
162 :     }
163 :    
164 :     my $weights = $options{ weights };
165 :    
166 :    
167 :     #---------------------------------------------------------------------------
168 :     # Options that are not protdist_neighbor options per se:
169 :     #---------------------------------------------------------------------------
170 :    
171 :     my $protdist = $options{ protdist } || 'protdist';
172 :    
173 :     my $neighbor = $options{ neighbor } || 'neighbor';
174 :    
175 :     my $tmp = $options{ tmp };
176 :    
177 :     my $tmp_dir = $options{ tmp_dir };
178 :    
179 :     my $tree_format = $options{ tree_format } =~ m/overbeek/i ? 'overbeek'
180 :     : $options{ tree_format } =~ m/gjo/i ? 'gjonewick'
181 :     : $options{ tree_format } =~ m/fig/i ? 'fig'
182 :     : 'overbeek'; # Default
183 :    
184 :     my $save_tmp = $tmp_dir && -d $tmp_dir;
185 :     if ( $tmp_dir )
186 :     {
187 :     if ( -d $tmp_dir ) { $save_tmp = 1 }
188 :     else { mkdir $tmp_dir }
189 :     }
190 :     else
191 :     {
192 :     $tmp = $tmp && -d $tmp ? $tmp
193 :     : -d '/tmp' ? '/tmp'
194 :     : '.';
195 :     my $int = int( 1000000000 * rand);
196 :     $tmp_dir = "$tmp/protdist_neighbor.$$.$int";
197 :     mkdir $tmp_dir;
198 :     }
199 :    
200 :     #---------------------------------------------------------------------------
201 :     # Write the files and run the program:
202 :     #---------------------------------------------------------------------------
203 :    
204 :     my $cwd = $ENV{ cwd } || `pwd`;
205 :     chomp $cwd;
206 :     chdir $tmp_dir;
207 :    
208 :     unlink 'outfile' if -f 'outfile'; # Just checking
209 :     unlink 'outtree' if -f 'outtree'; # ditto
210 :    
211 :     # protdist 3.66 has a serious bug when weights and categories are both
212 :     # used. So we pack our own data for this combination. Seems to be
213 :     # ineffectual.
214 :    
215 :     if ( $categories && $weights )
216 :     {
217 :     my $mask = $weights;
218 :     $mask =~ tr/0/\177/c; # Everything except 0 becomes X'FF'
219 :     $mask =~ tr/0/\000/; # 0 becomes X'FF'
220 :     @align = map { my ( $id, $seq ) = @$_;
221 :     [ $id, pack_by_mask( $seq, $mask ) ]
222 :     }
223 :     @align;
224 :     # Make a copy so that we do not clobber the calling program's copy
225 :     $categories = [ @$categories ];
226 :     $categories->[1] = pack_by_mask( $categories->[1], $mask );
227 :     $weights = undef;
228 :     }
229 :    
230 :     &write_infile( @align ) or print STDERR "protdist_neighbor::protdist_neighbor: Could write infile\n"
231 :     and chdir $cwd
232 :     and return ();
233 :    
234 :     open( PROTD, ">protdist_cmd" ) or print STDERR "protdist_neighbor::protdist_neighbor: Could not open command file for $protdist\n"
235 :     and chdir $cwd
236 :     and return ();
237 :    
238 :    
239 :     # Start writing optoins for protdist:
240 :    
241 :     if ( $categories )
242 :     {
243 :     &write_categories( $categories->[1] ) or print STDERR "protdist_neighbor::protdist_neighbor: Could not write categories\n"
244 :     and chdir $cwd
245 :     and return ();
246 :     print PROTD "C\n",
247 :     scalar @{$categories->[0]}, "\n",
248 :     join( ' ', map { sprintf( "%.6f", $_ ) } @{ $categories->[0] } ), "\n";
249 :     }
250 :    
251 :     if ( $invar_frac || $coef_of_var )
252 :     {
253 :     print PROTD "G\n";
254 :     print PROTD "G\n" if $invar_frac;
255 :     print PROTD "A\n", "$persistance\n" if $persistance;
256 :     }
257 :    
258 :     print PROTD "P\n" if $model =~ m/PMB/i;
259 :     print PROTD "P\nP\n" if $model =~ m/PAM/i;
260 :     print PROTD "P\nP\nP\n" if $model =~ m/Kimura/i;
261 :    
262 :     if ( $weights )
263 :     {
264 :     &write_weights( $weights ) or print STDERR "protdist_neighbor::protdist_neighbor: Could not write weights\n"
265 :     and chdir $cwd
266 :     and return ();
267 :     print PROTD "W\n";
268 :     }
269 :    
270 :     # All the options are written, try to lauch the run:
271 :    
272 :     print PROTD "Y\n";
273 :    
274 :     # Becuase of the options interface, these values have to be supplied after
275 :     # the Y:
276 :    
277 :     if ( $invar_frac || $coef_of_var )
278 :     {
279 :     print PROTD "$coef_of_var\n";
280 :     print PROTD "$invar_frac\n" if $invar_frac;
281 :     }
282 :    
283 :     close PROTD;
284 :    
285 :     system "$protdist < protdist_cmd > /dev/null; /bin/mv -f outfile infile";
286 :     # system "$protdist < protdist_cmd > /dev/null"; chdir $cwd; return;
287 :     # system "$protdist < protdist_cmd > /dev/null; /bin/cp infile protdist.infile; /bin/mv -f outfile infile";
288 :    
289 :     # Move on to neighbor:
290 :    
291 :     open( NEIGH, ">neigh_cmd" ) or print STDERR "protdist_neighbor::protdist_neighbor: Could not open neighbor command file\n"
292 :     and chdir $cwd
293 :     and return ();
294 :    
295 :    
296 :     # Start sending optoins ot program:
297 :    
298 :     print NEIGH "J\n", "$jumble_seed\n" if $jumble_seed;
299 :    
300 :     # All the options are written, try to launch the run:
301 :    
302 :     print NEIGH "Y\n";
303 :    
304 :     close NEIGH;
305 :    
306 :     system "$neighbor < neigh_cmd > /dev/null";
307 :    
308 :     my ( $tree ) = gjonewicklib::read_newick_tree( 'outtree' );
309 :     $tree or print STDERR "protdist_neighbor::protdist_neighbor: Could read neighbor outtree file\n"
310 :     and chdir $cwd
311 :     and return ();
312 :    
313 :     # We are done, go back to the original directory:
314 :    
315 :     chdir $cwd;
316 :    
317 :     # Returned trees have our labels:
318 :    
319 :     gjonewicklib::newick_relabel_nodes( $tree, \%id );
320 :    
321 :     if ( $tree_format =~ m/overbeek/i )
322 :     {
323 :     $tree = gjonewicklib::gjonewick_to_overbeek( $tree );
324 :     }
325 :    
326 :     system "/bin/rm -r $tmp_dir" if ! $save_tmp;
327 :    
328 :     return $tree;
329 :     }
330 :    
331 :     #===============================================================================
332 :     # $distances = protdist( \@align, \%options )
333 :     #===============================================================================
334 :     sub protdist
335 :     {
336 :     my $align;
337 :     if ( ref( $_[0] ) eq 'ARRAY' )
338 :     {
339 :     $align = shift @_;
340 :     ( $align && ( ref( $align ) eq 'ARRAY' ) )
341 :     || ( ( print STDERR "protdist_neighbor::protdist_neighbor() called without alignment\n" )
342 :     && ( return undef )
343 :     );
344 :     }
345 :    
346 :     my %options;
347 :     if ( $_[0] )
348 :     {
349 :     %options = ( ref( $_[0]) eq 'HASH' ) ? %{ $_[0] } : @_;
350 :     }
351 :    
352 :     #---------------------------------------------------------------------------
353 :     # Work on a copy of the alignment. Id is always first, seq is always last
354 :     #---------------------------------------------------------------------------
355 :    
356 :     $align ||= $options{ alignment } || $options{ align };
357 :    
358 :     my ( $seq, $id );
359 :     my %id;
360 :     my %local_id;
361 :     my $local_id = 'seq0000000';
362 :     my @align = map { $id = $_->[0];
363 :     $local_id++;
364 :     $id{ $local_id } = $id;
365 :     $local_id{ $id } = $local_id;
366 :     $seq = $_->[-1];
367 :     $seq =~ s/[BJOUZ]/X/gi; # Bad letters go to X
368 :     $seq =~ s/[^A-Z]/-/gi; # Anything else becomes -
369 :     [ $local_id, $seq ]
370 :     } @$align;
371 :    
372 :     #---------------------------------------------------------------------------
373 :     # Process protdist options:
374 :     #---------------------------------------------------------------------------
375 :    
376 :     my $categories = $options{ categories }; # [ [ cat_rates ], site_cats ]
377 :     if ( $categories )
378 :     {
379 :     if ( ref( $categories ) ne 'ARRAY'
380 :     || ! ( ( @$categories == 2 ) || ( ( @$categories == 3 ) && ( shift @$categories ) ) )
381 :     || ref( $categories->[0] ) ne 'ARRAY'
382 :     )
383 :     {
384 :     print STDERR "proml::proml categories option value must be [ [ cat_rate1, ... ], site_categories ]\n";
385 :     return undef;
386 :     }
387 :    
388 :     # Rate values cannot have very many decimal places or proml can't read it:
389 :    
390 :     @{$categories->[0]} = map { sprintf "%.6f", $_ } @{$categories->[0]};
391 :     }
392 :    
393 :     my $coef_of_var = $options{ coef_of_var }
394 :     || ( $options{ alpha } && ( $options{ alpha } > 0) && ( 1 / sqrt( $options{ alpha } ) ) )
395 :     || 0;
396 :     if ( $coef_of_var < 0 )
397 :     {
398 :     print STDERR "protdist_neighbor::protdist_neighbor coef_of_var option value must be >= 0\n";
399 :     return undef;
400 :     }
401 :    
402 :     my $invar_frac = $options{ invar_frac } || 0;
403 :     if ( $invar_frac && ( $invar_frac < 0 || $invar_frac >= 1 ) )
404 :     {
405 :     print STDERR "protdist_neighbor::protdist_neighbor invar_frac option value must be >= 0 and < 1\n";
406 :     return undef;
407 :     }
408 :    
409 :     my $model = ( $options{ model } =~ m/PAM/i ) ? 'PAM'
410 :     : ( $options{ model } =~ m/Dayhoff/i ) ? 'PAM'
411 :     : ( $options{ model } =~ m/PMB/i ) ? 'PMB'
412 :     : ( $options{ model } =~ m/Henikoff/i ) ? 'PMB'
413 :     : ( $options{ model } =~ m/Tillier/i ) ? 'PMB'
414 :     : ( $options{ model } =~ m/JTT/i ) ? 'JTT'
415 :     : ( $options{ model } =~ m/Jones/i ) ? 'JTT'
416 :     : ( $options{ model } =~ m/Taylor/i ) ? 'JTT'
417 :     : ( $options{ model } =~ m/Thornton/i ) ? 'JTT'
418 :     : ( $options{ model } =~ m/Kimura/i ) ? 'Kimura'
419 :     : 'JTT';
420 :    
421 :     my $persistance = $options{ persistance } || 0;
422 :     if ( $persistance && ( $persistance <= 1 ) )
423 :     {
424 :     print STDERR "protdist_neighbor::protdist_neighbor persistance option value must be > 1\n";
425 :     return undef;
426 :     }
427 :    
428 :     my $weights = $options{ weights };
429 :    
430 :    
431 :     #---------------------------------------------------------------------------
432 :     # Options that are not protdist_neighbor options per se:
433 :     #---------------------------------------------------------------------------
434 :    
435 :     my $protdist = $options{ protdist } || 'protdist';
436 :    
437 :     my $neighbor = $options{ neighbor } || 'neighbor';
438 :    
439 :     my $tmp = $options{ tmp };
440 :    
441 :     my $tmp_dir = $options{ tmp_dir };
442 :    
443 :     my $tree_format = $options{ tree_format } =~ m/overbeek/i ? 'overbeek'
444 :     : $options{ tree_format } =~ m/gjo/i ? 'gjonewick'
445 :     : $options{ tree_format } =~ m/fig/i ? 'fig'
446 :     : 'overbeek'; # Default
447 :    
448 :     my $save_tmp = $tmp_dir && -d $tmp_dir;
449 :     if ( $tmp_dir )
450 :     {
451 :     if ( -d $tmp_dir ) { $save_tmp = 1 }
452 :     else { mkdir $tmp_dir }
453 :     }
454 :     else
455 :     {
456 :     $tmp = $tmp && -d $tmp ? $tmp
457 :     : -d '/tmp' ? '/tmp'
458 :     : '.';
459 :     my $int = int( 1000000000 * rand);
460 :     $tmp_dir = "$tmp/protdist_neighbor.$$.$int";
461 :     mkdir $tmp_dir;
462 :     }
463 :    
464 :     #---------------------------------------------------------------------------
465 :     # Write the files and run the program:
466 :     #---------------------------------------------------------------------------
467 :    
468 :     my $cwd = $ENV{ cwd } || `pwd`;
469 :     chomp $cwd;
470 :     chdir $tmp_dir;
471 :    
472 :     unlink 'outfile' if -f 'outfile'; # Just checking
473 :     unlink 'outtree' if -f 'outtree'; # ditto
474 :    
475 :     # protdist 3.66 has a serious bug when weights and categories are both
476 :     # used. So we pack our own data for this combination. Seems to be
477 :     # ineffectual.
478 :    
479 :     if ( $categories && $weights )
480 :     {
481 :     my $mask = $weights;
482 :     $mask =~ tr/0/\177/c; # Everything except 0 becomes X'FF'
483 :     $mask =~ tr/0/\000/; # 0 becomes X'FF'
484 :     @align = map { my ( $id, $seq ) = @$_;
485 :     [ $id, pack_by_mask( $seq, $mask ) ]
486 :     }
487 :     @align;
488 :     # Make a copy so that we do not clobber the calling program's copy
489 :     $categories = [ @$categories ];
490 :     $categories->[1] = pack_by_mask( $categories->[1], $mask );
491 :     $weights = undef;
492 :     }
493 :    
494 :     &write_infile( @align ) or print STDERR "protdist_neighbor::protdist_neighbor: Could write infile\n"
495 :     and chdir $cwd
496 :     and return undef;
497 :    
498 :     open( PROTD, ">protdist_cmd" ) or print STDERR "protdist_neighbor::protdist_neighbor: Could not open command file for $protdist\n"
499 :     and chdir $cwd
500 :     and return undef;
501 :    
502 :    
503 :     # Start writing optoins for protdist:
504 :    
505 :     if ( $categories )
506 :     {
507 :     &write_categories( $categories->[1] ) or print STDERR "protdist_neighbor::protdist_neighbor: Could not write categories\n"
508 :     and chdir $cwd
509 :     and return undef;
510 :     print PROTD "C\n",
511 :     scalar @{$categories->[0]}, "\n",
512 :     join( ' ', map { sprintf( "%.6f", $_ ) } @{ $categories->[0] } ), "\n";
513 :     }
514 :    
515 :     if ( $invar_frac || $coef_of_var )
516 :     {
517 :     print PROTD "G\n";
518 :     print PROTD "G\n" if $invar_frac;
519 :     print PROTD "A\n", "$persistance\n" if $persistance;
520 :     }
521 :    
522 :     print PROTD "P\n" if $model =~ m/PMB/i;
523 :     print PROTD "P\nP\n" if $model =~ m/PAM/i;
524 :     print PROTD "P\nP\nP\n" if $model =~ m/Kimura/i;
525 :    
526 :     if ( $weights )
527 :     {
528 :     &write_weights( $weights ) or print STDERR "protdist_neighbor::protdist_neighbor: Could not write weights\n"
529 :     and chdir $cwd
530 :     and return undef;
531 :     print PROTD "W\n";
532 :     }
533 :    
534 :     # All the options are written, try to lauch the run:
535 :    
536 :     print PROTD "Y\n";
537 :    
538 :     # Becuase of the options interface, these values have to be supplied after
539 :     # the Y:
540 :    
541 :     if ( $invar_frac || $coef_of_var )
542 :     {
543 :     print PROTD "$coef_of_var\n";
544 :     print PROTD "$invar_frac\n" if $invar_frac;
545 :     }
546 :    
547 :     close PROTD;
548 :    
549 :     system "$protdist < protdist_cmd > /dev/null";
550 :    
551 :     my $distances = read_distances();
552 :     $distances or print STDERR "protdist_neighbor::protdist: Could not read 'outfile'\n"
553 :     and chdir $cwd
554 :     and return undef;
555 :    
556 :     # We are done, go back to the original directory:
557 :    
558 :     chdir $cwd;
559 :    
560 :     system "/bin/rm -r $tmp_dir" if ! $save_tmp;
561 :    
562 :     return $distances;
563 :     }
564 :    
565 :    
566 :     #-------------------------------------------------------------------------------
567 :     # Auxiliary functions:
568 :     #-------------------------------------------------------------------------------
569 :    
570 :     sub pack_by_mask
571 :     {
572 :     my ( $seq, $mask ) = @_;
573 :     $seq &= $mask; # Mask the string
574 :     $seq =~ tr/\000//d; # Compress out X'00'
575 :     $seq;
576 :     }
577 :    
578 :    
579 :     sub write_infile
580 :     {
581 :     open( INFILE, '>infile' ) or return 0;
582 :     print INFILE scalar @_, ' ', length( $_[0]->[1] ), "\n";
583 :     foreach ( @_ ) { printf INFILE "%-10s %s\n", @$_ }
584 :     close( INFILE );
585 :     }
586 :    
587 :    
588 :     sub write_categories
589 :     {
590 :     my $categories = shift;
591 :     open( CATEGORIES, '>categories' ) or return 0;
592 :     print CATEGORIES "$categories\n";
593 :     close( CATEGORIES );
594 :     }
595 :    
596 :    
597 :     sub write_weights
598 :     {
599 :     my $weights = shift;
600 :     open( WEIGHTS, '>weights' ) or return 0;
601 :     print WEIGHTS "$weights\n";
602 :     close( WEIGHTS );
603 :     }
604 :    
605 :    
606 :     sub read_distances
607 :     {
608 :     my @distances;
609 :     open( DISTS, '<outfile' ) or return undef;
610 :     $_ = <DISTS> or return undef;
611 :     my ( $ndist ) = $_ =~ m/(\d+)/;
612 :     for ( my $i = 0; $i < $ndist; $i++ )
613 :     {
614 :     $_ = <DISTS> or return undef;
615 :     my ( undef, @row ) = split;
616 :     while ( @row < $ndist )
617 :     {
618 :     $_ = <DISTS> or return undef;
619 :     push @row, split;
620 :     }
621 :    
622 :     push @distances, \@row;
623 :     }
624 :     close( DISTS );
625 :     \@distances;
626 :     }
627 :    
628 :     1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3