~ubuntu-branches/ubuntu/raring/bioperl/raring

« back to all changes in this revision

Viewing changes to Bio/AlignIO/phylip.pm

  • Committer: Bazaar Package Importer
  • Author(s): Charles Plessy
  • Date: 2008-03-18 14:44:57 UTC
  • mfrom: (4 hardy)
  • mto: This revision was merged to the branch mainline in revision 6.
  • Revision ID: james.westby@ubuntu.com-20080318144457-1jjoztrvqwf0gruk
* debian/control:
  - Removed MIA Matt Hope (dopey) from the Uploaders field.
    Thank you for your work, Matt. I hope you are doing well.
  - Downgraded some recommended package to the 'Suggests' priority,
    according to the following discussion on Upstream's mail list.
    http://bioperl.org/pipermail/bioperl-l/2008-March/027379.html
    (Closes: #448890)
* debian/copyright converted to machine-readable format.

Show diffs side-by-side

added added

removed removed

Lines of Context:
1
 
# $Id: phylip.pm,v 1.26 2003/03/28 22:39:54 jason Exp $
 
1
# $Id: phylip.pm,v 1.36.4.3 2006/10/02 23:10:12 sendu Exp $
2
2
#
3
3
# BioPerl module for Bio::AlignIO::phylip
4
4
#
15
15
 
16
16
    use Bio::AlignIO;
17
17
    use Bio::SimpleAlign;
18
 
        #you can set the name length to something other than the default 10
19
 
        #if you use a version of phylip (hacked) that accepts ids > 10
 
18
    #you can set the name length to something other than the default 10
 
19
    #if you use a version of phylip (hacked) that accepts ids > 10
20
20
    my $phylipstream = new Bio::AlignIO(-format  => 'phylip',
21
 
                                        -fh      => \*STDOUT,
22
 
                                        -idlength=>30);
 
21
                                        -fh      => \*STDOUT,
 
22
                                        -idlength=>30);
23
23
    # convert data from one format to another
24
24
    my $gcgstream     =  new Bio::AlignIO(-format => 'msf',
25
 
                                          -file   => 't/data/cysprot1a.msf');
 
25
                                          -file   => 't/data/cysprot1a.msf');
26
26
 
27
27
    while( my $aln = $gcgstream->next_aln ) {
28
 
        $phylipstream->write_aln($aln);
 
28
        $phylipstream->write_aln($aln);
29
29
    }
30
30
 
31
 
    # do it again with phylip sequential format format 
 
31
    # do it again with phylip sequential format format
32
32
    $phylipstream->interleaved(0);
33
33
    # can also initialize the object like this
34
34
    $phylipstream = new Bio::AlignIO(-interleaved => 0,
35
 
                                     -format => 'phylip',
36
 
                                     -fh   => \*STDOUT,
37
 
                                     -idlength=>10);
 
35
                                     -format => 'phylip',
 
36
                                     -fh   => \*STDOUT,
 
37
                                     -idlength=>10);
38
38
    $gcgstream     =  new Bio::AlignIO(-format => 'msf',
39
 
                                       -file   => 't/data/cysprot1a.msf');    
 
39
                                       -file   => 't/data/cysprot1a.msf');
40
40
 
41
41
    while( my $aln = $gcgstream->next_aln ) {
42
 
        $phylipstream->write_aln($aln);
 
42
        $phylipstream->write_aln($aln);
43
43
    }
44
44
 
45
45
=head1 DESCRIPTION
56
56
=head2 Reporting Bugs
57
57
 
58
58
Report bugs to the Bioperl bug tracking system to help us keep track
59
 
 the bugs and their resolution.
60
 
 Bug reports can be submitted via email or the web:
 
59
the bugs and their resolution. Bug reports can be submitted via the
 
60
web:
61
61
 
62
 
  bioperl-bugs@bio.perl.org
63
 
  http://bugzilla.bioperl.org/
 
62
  http://bugzilla.open-bio.org/
64
63
 
65
64
=head1 AUTHORS - Heikki Lehvaslaiho and Jason Stajich
66
65
 
67
 
Email: heikki@ebi.ac.uk
68
 
Email: jason@bioperl.org
 
66
Email: heikki at ebi.ac.uk
 
67
Email: jason at bioperl.org
69
68
 
70
69
=head1 APPENDIX
71
70
 
77
76
# Let the code begin...
78
77
 
79
78
package Bio::AlignIO::phylip;
80
 
use vars qw(@ISA $DEFAULTIDLENGTH $DEFAULTLINELEN);
 
79
use vars qw($DEFAULTIDLENGTH $DEFAULTLINELEN $DEFAULTTAGLEN);
81
80
use strict;
82
81
 
83
82
use Bio::SimpleAlign;
84
 
use Bio::AlignIO;
85
 
 
86
 
@ISA = qw(Bio::AlignIO);
87
 
 
88
 
BEGIN { 
 
83
use POSIX; # for the rounding call
 
84
 
 
85
use base qw(Bio::AlignIO);
 
86
 
 
87
BEGIN {
89
88
    $DEFAULTIDLENGTH = 10;
90
89
    $DEFAULTLINELEN = 60;
 
90
    $DEFAULTTAGLEN = 10;
91
91
}
92
92
 
93
93
=head2 new
100
100
 Function: Initialize a new L<Bio::AlignIO::phylip> reader or writer
101
101
 Returns : L<Bio::AlignIO> object
102
102
 Args    : [specific for writing of phylip format files]
103
 
           -idlength => integer - length of the id (will pad w/ 
104
 
                                                    spaces if needed) 
105
 
           -interleaved => boolean - whether or not write as interleaved 
 
103
           -idlength => integer - length of the id (will pad w/
 
104
                                                    spaces if needed)
 
105
           -interleaved => boolean - whether or not write as interleaved
106
106
                                     or sequential format
107
 
           -linelength  => integer of how long a sequence lines should be 
 
107
           -line_length  => integer of how long a sequence lines should be
108
108
           -idlinebreak => insert a line break after the sequence id
109
 
                           so that sequence starts on the next line 
 
109
                           so that sequence starts on the next line
 
110
           -flag_SI => whether or not write a "S" or "I" just after
 
111
                       the num.seq. and line len., in the first line
 
112
           -tag_length => integer of how long the tags have to be in
 
113
                         each line between the space separator. set it
 
114
                         to 0 to have 1 tag only.
 
115
           -wrap_sequential => boolean for whether or not sequential
 
116
                                   format should be broken up or a single line
 
117
                                   default is false (single line)
110
118
 
111
119
=cut
112
120
 
115
123
  $self->SUPER::_initialize(@args);
116
124
 
117
125
  my ($interleave,$linelen,$idlinebreak,
118
 
      $idlength) = $self->_rearrange([qw(INTERLEAVED 
119
 
                                         LINELENGTH
120
 
                                         IDLINEBREAK
121
 
                                         IDLENGTH)],@args);
 
126
      $idlength, $flag_SI, $tag_length,$ws) =
 
127
          $self->_rearrange([qw(INTERLEAVED
 
128
                                LINE_LENGTH
 
129
                                IDLINEBREAK
 
130
                                IDLENGTH
 
131
                                FLAG_SI
 
132
                                TAG_LENGTH
 
133
                                WRAP_SEQUENTIAL)],@args);
122
134
  $self->interleaved(1) if( $interleave || ! defined $interleave);
123
135
  $self->idlength($idlength || $DEFAULTIDLENGTH);
124
136
  $self->id_linebreak(1) if( $idlinebreak );
125
137
  $self->line_length($linelen) if defined $linelen && $linelen > 0;
 
138
  $self->flag_SI(1) if ( $flag_SI );
 
139
  $self->tag_length($tag_length) if ( $tag_length || $DEFAULTTAGLEN );
 
140
  $self->wrap_sequential($ws ? 1 : 0);
126
141
  1;
127
142
}
128
143
 
134
149
           Throws an exception if trying to read in PHYLIP
135
150
           sequential format.
136
151
 Returns : L<Bio::SimpleAlign> object
137
 
 Args    : 
 
152
 Args    :
138
153
 
139
154
=cut
140
155
 
143
158
    my $entry;
144
159
    my ($seqcount, $residuecount, %hash, $name,$str,
145
160
        @names,$seqname,$start,$end,$count,$seq);
146
 
    
 
161
 
147
162
    my $aln =  Bio::SimpleAlign->new(-source => 'phylip');
148
 
    $entry = $self->_readline and 
 
163
    $entry = $self->_readline and
149
164
        ($seqcount, $residuecount) = $entry =~ /\s*(\d+)\s+(\d+)/;
150
165
    return 0 unless $seqcount and $residuecount;
151
 
    
 
166
 
152
167
    # first alignment section
153
168
    my $idlen = $self->idlength;
154
169
    $count = 0;
155
170
    my $iter = 1;
156
 
    my $non_interleaved = ! $self->interleaved ;
157
 
    
 
171
    my $interleaved = $self->interleaved;
158
172
    while( $entry = $self->_readline) {
159
 
        last if( $entry =~ /^\s?$/ && ! $non_interleaved );
 
173
        last if( $entry =~ /^\s?$/ && $interleaved );
160
174
 
161
 
        if( $entry =~ /^\s+(\d+)\s+(\d+)\s*$/) { 
 
175
        if( $entry =~ /^\s+(\d+)\s+(\d+)\s*$/) {
162
176
            $self->_pushback($entry);
163
177
            last;
164
178
        }
165
179
        if( $entry =~ /^\s+(.+)$/ ) {
 
180
            $interleaved = 0;
166
181
            $str = $1;
167
 
            $non_interleaved = 1;
168
182
            $str =~ s/\s//g;
169
 
            unless( ! $non_interleaved ) {
170
 
                $count = scalar @names;
171
 
                $hash{$count} .= $str;
172
 
            } else { 
173
 
                $hash{$iter++} .= $str;
174
 
                $iter = 1 if $iter > $count;
175
 
            }
176
 
        } elsif( $entry =~ /^(.{$idlen})\s+(.*)\s$/ ||           
 
183
            $count = scalar @names;
 
184
            $hash{$count} .= $str;
 
185
 
 
186
        } elsif( $entry =~ /^(.{$idlen})\s+(.*)\s$/ ||
177
187
                 $entry =~ /^(.{$idlen})(\S{$idlen}\s+.+)\s$/ # Handle weirdnes s when id is too long
178
188
                 ) {
179
189
            $name = $1;
180
190
            $str = $2;
181
191
            $name =~ s/[\s\/]/_/g;
182
192
            $name =~ s/_+$//; # remove any trailing _'s
 
193
 
183
194
            push @names, $name;
184
195
            $str =~ s/\s//g;
185
196
            $count = scalar @names;
186
197
            $hash{$count} = $str;
 
198
        } elsif( $interleaved ) {
 
199
            if( $entry =~ /^(\S+)\s+(.+)/ ||
 
200
                $entry =~ /^(.{$idlen})(.*)\s$/ ) {
 
201
                $name = $1;
 
202
                $str = $2;
 
203
                $name =~ s/[\s\/]/_/g;
 
204
                $name =~ s/_+$//; # remove any trailing _'s
 
205
                push @names, $name;
 
206
                $str =~ s/\s//g;
 
207
                $count = scalar @names;
 
208
                $hash{$count} = $str;
 
209
            } else {
 
210
                $self->debug("unmatched line: $entry");
 
211
            }
187
212
        }
188
 
        $self->throw("Not a valid interleaved PHYLIP file!") if $count > $seqcount; 
 
213
        $self->throw("Not a valid interleaved PHYLIP file!") if $count > $seqcount;
189
214
    }
190
 
    
191
 
    unless( $non_interleaved ) {    
 
215
 
 
216
    if( $interleaved ) {
192
217
        # interleaved sections
193
218
        $count = 0;
194
219
        while( $entry = $self->_readline) {
195
 
            
196
220
            # finish current entry
197
221
            if($entry =~/\s*\d+\s+\d+/){
198
222
                $self->_pushback($entry);
205
229
                $count++;
206
230
                $hash{$count} .= $str;
207
231
            };
208
 
            $self->throw("Not a valid interleaved PHYLIP file!") if $count > $seqcount; 
 
232
            $self->throw("Not a valid interleaved PHYLIP file! [$count,$seqcount] ($entry)") if $count > $seqcount;
209
233
        }
210
234
    }
211
235
    return 0 if scalar @names < 1;
212
 
    
 
236
 
213
237
    # sequence creation
214
238
    $count = 0;
215
239
    foreach $name ( @names ) {
226
250
            $end = length($str);
227
251
        }
228
252
        # consistency test
229
 
        $self->throw("Length of sequence [$seqname] is not [$residuecount] it is ".CORE::length($hash{$count})."! ") 
230
 
            unless CORE::length($hash{$count}) == $residuecount; 
231
 
        
 
253
        $self->throw("Length of sequence [$seqname] is not [$residuecount] it is ".CORE::length($hash{$count})."! ")
 
254
            unless CORE::length($hash{$count}) == $residuecount;
 
255
 
232
256
       $seq = new Bio::LocatableSeq('-seq'=>$hash{$count},
233
257
                                    '-id'=>$seqname,
234
258
                                    '-start'=>$start,
235
259
                                    '-end'=>$end,
236
 
                                    );
237
 
        
238
 
       $aln->add_seq($seq);
 
260
                                   );
 
261
        $aln->add_seq($seq);
239
262
 
240
263
   }
241
264
   return $aln;
257
280
    my $count = 0;
258
281
    my $wrapped = 0;
259
282
    my $maxname;
 
283
    my $width = $self->line_length();
260
284
    my ($length,$date,$name,$seq,$miss,$pad,
261
 
        %hash,@arr,$tempcount,$index,$idlength);
262
 
    
 
285
        %hash,@arr,$tempcount,$index,$idlength,$flag_SI,$line_length, $tag_length);
 
286
 
263
287
    foreach my $aln (@aln) {
264
 
        if( ! $aln || ! $aln->isa('Bio::Align::AlignI')  ) { 
 
288
        if( ! $aln || ! $aln->isa('Bio::Align::AlignI')  ) {
265
289
            $self->warn("Must provide a Bio::Align::AlignI object when calling write_aln");
266
290
            next;
267
291
        }
268
 
        $self->throw("All sequences in the alignment must be the same length") 
 
292
        $self->throw("All sequences in the alignment must be the same length")
269
293
            unless $aln->is_flush(1) ;
270
294
 
 
295
        $flag_SI = $self->flag_SI();
271
296
        $aln->set_displayname_flat(); # plain
272
297
        $length  = $aln->length();
273
 
        $self->_print (sprintf(" %s %s\n", $aln->no_sequences, $aln->length));
 
298
        if ($flag_SI) {
 
299
            if ($self->interleaved() ) {
 
300
                $self->_print (sprintf(" %s %s I\n", $aln->no_sequences, $aln->length));
 
301
            } else {
 
302
                $self->_print (sprintf(" %s %s S\n", $aln->no_sequences, $aln->length));
 
303
            }
 
304
        } else {
 
305
            $self->_print (sprintf(" %s %s\n", $aln->no_sequences, $aln->length));
 
306
        }
274
307
 
275
 
        $idlength = $self->idlength();  
 
308
        $idlength = $self->idlength();
 
309
        $line_length = $self->line_length();
 
310
        $tag_length = $self->tag_length();
276
311
        foreach $seq ( $aln->each_seq() ) {
277
312
            $name = $aln->displayname($seq->get_nse);
278
313
            $name = substr($name, 0, $idlength) if length($name) > $idlength;
279
 
            $name = sprintf("%-".$idlength."s",$name);      
 
314
            $name = sprintf("%-".$idlength."s",$name);
280
315
            if( $self->interleaved() ) {
281
316
                $name .= '   ' ;
282
 
            } elsif( $self->id_linebreak) { 
283
 
                $name .= "\n"; 
 
317
            } elsif( $self->id_linebreak) {
 
318
                $name .= "\n";
284
319
            }
285
320
 
286
 
      #phylip needs dashes not dots 
287
 
      my $seq = $seq->seq();
288
 
      $seq=~s/\./-/g;
 
321
            #phylip needs dashes not dots
 
322
            my $seq = $seq->seq();
 
323
            $seq =~ s/\./-/g;
289
324
            $hash{$name} = $seq;
290
325
            push(@arr,$name);
291
326
        }
292
327
 
293
328
        if( $self->interleaved() ) {
294
 
            while( $count < $length ) { 
295
 
                
 
329
            my $numtags;
 
330
            if ($tag_length <= $line_length) {
 
331
                $numtags = floor($line_length/$tag_length);
 
332
                $line_length = $tag_length*$numtags;
 
333
            } else {
 
334
                $numtags = 1;
 
335
            }
 
336
            while( $count < $length ) {
 
337
 
296
338
                # there is another block to go!
297
339
                foreach $name ( @arr ) {
298
340
                    my $dispname = $name;
299
341
                    $dispname = '' if $wrapped;
300
342
                    $self->_print (sprintf("%".($idlength+3)."s",$dispname));
301
343
                    $tempcount = $count;
302
 
                    $index = 0;
303
 
                    while( ($tempcount + $idlength < $length) && ($index < 5)  ) {
 
344
                    $index = 0;
 
345
                    $self->debug("residue count: $count\n") if ($count%100000 == 0);
 
346
                    while( ($tempcount + $tag_length < $length) &&
 
347
                           ($index < $numtags)  ) {
304
348
                        $self->_print (sprintf("%s ",substr($hash{$name},
305
349
                                                            $tempcount,
306
 
                                                            $idlength)));
307
 
                        $tempcount += $idlength;
 
350
                                                            $tag_length)));
 
351
                        $tempcount += $tag_length;
308
352
                        $index++;
309
353
                    }
310
354
                    # last
311
 
                    if( $index < 5) {
 
355
                    if( $index < $numtags) {
312
356
                        # space to print!
313
357
                        $self->_print (sprintf("%s ",substr($hash{$name},
314
358
                                                            $tempcount)));
315
 
                        $tempcount += $idlength;
 
359
                        $tempcount += $tag_length;
316
360
                    }
317
361
                    $self->_print ("\n");
318
362
                }
319
363
                $self->_print ("\n");
320
364
                $count = $tempcount;
321
365
                $wrapped = 1;
322
 
            }                   
 
366
            }
323
367
        } else {
324
368
            foreach $name ( @arr ) {
325
369
                my $dispname = $name;
326
 
                $dispname = '' if $wrapped;
327
 
                $self->_print (sprintf("%s%s\n",$dispname,$hash{$name}));
328
 
            }   
 
370
                my $line = sprintf("%s%s\n",$dispname,$hash{$name});
 
371
                if( $self->wrap_sequential ) {
 
372
                    $line =~ s/(.{1,$width})/$1\n/g;
 
373
                }
 
374
                $self->_print ($line);
 
375
            }
329
376
        }
330
377
    }
331
378
    $self->flush if $self->_flush_on_write && defined $self->_fh;
346
393
sub interleaved{
347
394
   my ($self,$value) = @_;
348
395
   my $previous = $self->{'_interleaved'};
349
 
   if( defined $value ) { 
 
396
   if( defined $value ) {
350
397
       $self->{'_interleaved'} = $value;
351
398
   }
352
399
   return $previous;
353
400
}
354
401
 
 
402
=head2 flag_SI
 
403
 
 
404
 Title   : flag_SI
 
405
 Usage   : my $flag = $obj->flag_SI
 
406
 Function: Get/Set if the Sequential/Interleaved flag has to be shown
 
407
           after the number of sequences and sequence length
 
408
 Example :
 
409
 Returns : boolean
 
410
 Args    : boolean
 
411
 
 
412
 
 
413
=cut
 
414
 
 
415
sub flag_SI{
 
416
   my ($self,$value) = @_;
 
417
   my $previous = $self->{'_flag_SI'};
 
418
   if( defined $value ) {
 
419
       $self->{'_flag_SI'} = $value;
 
420
   }
 
421
   return $previous;
 
422
}
 
423
 
355
424
=head2 idlength
356
425
 
357
426
 Title   : idlength
358
 
 Usage   : my $idlength = $obj->interleaved
359
 
 Function: Get/Set value of id length 
360
 
 Returns : string 
361
 
 Args    : string 
 
427
 Usage   : my $idlength = $obj->idlength
 
428
 Function: Get/Set value of id length
 
429
 Returns : string
 
430
 Args    : string
362
431
 
363
432
 
364
433
=cut
375
444
 
376
445
 Title   : line_length
377
446
 Usage   : $obj->line_length($newval)
378
 
 Function: 
 
447
 Function:
379
448
 Returns : value of line_length
380
449
 Args    : newvalue (optional)
381
450
 
385
454
sub line_length{
386
455
   my ($self,$value) = @_;
387
456
   if( defined $value) {
388
 
      $self->{'line_length'} = $value;
389
 
    }
390
 
    return $self->{'line_length'} || $DEFAULTLINELEN;
391
 
 
392
 
}
 
457
      $self->{'_line_length'} = $value;
 
458
    }
 
459
    return $self->{'_line_length'} || $DEFAULTLINELEN;
 
460
 
 
461
}
 
462
 
 
463
=head2 tag_length
 
464
 
 
465
 Title   : tag_length
 
466
 Usage   : $obj->tag_length($newval)
 
467
 Function:
 
468
 Example : my $tag_length = $obj->tag_length
 
469
 Returns : value of the length for each space-separated tag in a line
 
470
 Args    : newvalue (optional) - set to zero to have one tag per line
 
471
 
 
472
 
 
473
=cut
 
474
 
 
475
sub tag_length{
 
476
   my ($self,$value) = @_;
 
477
   if( defined $value) {
 
478
      $self->{'_tag_length'} = $value;
 
479
    }
 
480
    return $self->{'_tag_length'} || $DEFAULTTAGLEN;
 
481
}
 
482
 
393
483
 
394
484
=head2 id_linebreak
395
485
 
396
486
 Title   : id_linebreak
397
487
 Usage   : $obj->id_linebreak($newval)
398
 
 Function: 
 
488
 Function:
399
489
 Returns : value of id_linebreak
400
490
 Args    : newvalue (optional)
401
491
 
410
500
    return $self->{'_id_linebreak'} || 0;
411
501
}
412
502
 
 
503
 
 
504
=head2 wrap_sequential
 
505
 
 
506
 Title   : wrap_sequential
 
507
 Usage   : $obj->wrap_sequential($newval)
 
508
 Function:
 
509
 Returns : value of wrap_sequential
 
510
 Args    : newvalue (optional)
 
511
 
 
512
 
 
513
=cut
 
514
 
 
515
sub wrap_sequential{
 
516
   my ($self,$value) = @_;
 
517
   if( defined $value) {
 
518
      $self->{'_wrap_sequential'} = $value;
 
519
    }
 
520
    return $self->{'_wrap_sequential'} || 0;
 
521
}
 
522
 
413
523
1;