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

« back to all changes in this revision

Viewing changes to Bio/AlignIO/fasta.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: fasta.pm,v 1.14 2003/07/19 22:35:44 jason Exp $
 
1
# $Id: fasta.pm,v 1.27.4.1 2006/10/02 23:10:12 sendu Exp $
2
2
#
3
3
# BioPerl module for Bio::AlignIO::fasta
4
 
 
5
 
#       based on the Bio::SeqIO::fasta module
6
 
#       by Ewan Birney <birney@sanger.ac.uk>
7
 
#       and Lincoln Stein  <lstein@cshl.org>
8
 
#
9
 
#       and the SimpleAlign.pm module of Ewan Birney
10
4
#
11
5
# Copyright Peter Schattner
12
6
#
13
7
# You may distribute this module under the same terms as perl itself
14
 
# _history
15
 
# September 5, 2000
16
8
# POD documentation - main docs before the code
17
9
 
18
10
=head1 NAME
19
11
 
20
 
Bio::AlignIO::fasta - FastA MSA Sequence input/output stream
 
12
Bio::AlignIO::fasta - fasta MSA Sequence input/output stream
21
13
 
22
14
=head1 SYNOPSIS
23
15
 
24
 
Do not use this module directly.  Use it via the L<Bio::AlignIO> class.
 
16
Do not use this module directly.  Use it via the L<Bio::AlignIO> 
 
17
class.
25
18
 
26
19
=head1 DESCRIPTION
27
20
 
28
21
This object can transform L<Bio::SimpleAlign> objects to and from
29
 
fasta flat file databases.  This is for the fasta sequence format NOT
30
 
FastA analysis program.  To process the pairwise alignments from a
 
22
fasta flat file databases.  This is for the fasta alignment format, not
 
23
for the FastA sequence analysis program.  To process the alignments from
31
24
FastA (FastX, FastN, FastP, tFastA, etc) use the Bio::SearchIO module.
32
25
 
33
26
=head1 FEEDBACK
35
28
=head2 Reporting Bugs
36
29
 
37
30
Report bugs to the Bioperl bug tracking system to help us keep track
38
 
the bugs and their resolution.  Bug reports can be submitted via email
39
 
or the web:
40
 
 
41
 
  bioperl-bugs@bio.perl.org
42
 
  http://bugzilla.bioperl.org/
43
 
 
44
 
=head1 AUTHORS - Peter Schattner
45
 
 
46
 
Email: schattner@alum.mit.edu
47
 
 
 
31
the bugs and their resolution.  Bug reports can be submitted via the
 
32
web:
 
33
 
 
34
  http://bugzilla.open-bio.org/
 
35
 
 
36
=head1 AUTHORS
 
37
 
 
38
Peter Schattner
48
39
 
49
40
=head1 APPENDIX
50
41
 
56
47
# Let the code begin...
57
48
 
58
49
package Bio::AlignIO::fasta;
59
 
use vars qw(@ISA);
 
50
use vars qw($WIDTH);
60
51
use strict;
61
52
 
62
 
use Bio::AlignIO;
63
 
use Bio::SimpleAlign;
64
 
 
65
 
@ISA = qw(Bio::AlignIO);
66
 
 
 
53
 
 
54
use base qw(Bio::AlignIO);
 
55
$WIDTH = 60;
67
56
 
68
57
=head2 next_aln
69
58
 
70
59
 Title   : next_aln
71
 
 Usage   : $aln = $stream->next_aln()
 
60
 Usage   : $aln = $stream->next_aln
72
61
 Function: returns the next alignment in the stream.
73
 
 Returns : L<Bio::Align::AlignI> object - returns 0 on end of file
74
 
            or on error
75
 
 Args    : NONE
 
62
 Returns : Bio::Align::AlignI object - returns 0 on end of file
 
63
                or on error
 
64
 Args    : -width => optional argument to specify the width sequence
 
65
           will be written (60 chars by default)
 
66
 
 
67
See L<Bio::Align::AlignI>
76
68
 
77
69
=cut
78
70
 
79
71
sub next_aln {
80
 
    my $self = shift;
81
 
    my $entry;
82
 
    my ($start,$end,$name,$seqname,$seq,$seqchar,$tempname,$tempdesc,
83
 
        %align,$desc);
84
 
    my $aln =  Bio::SimpleAlign->new();
85
 
    my $maxlen;
86
 
    while(defined ($entry = $self->_readline) ) {
87
 
        if( $entry =~ s/^>(\S+)\s*// ) {
88
 
            $tempname  = $1;
89
 
            chomp($entry);
90
 
            $tempdesc  = $entry;
91
 
            if( defined $name ) {
92
 
                # put away last name and sequence
93
 
                if( $name =~ /(\S+)\/(\d+)-(\d+)/ ) {
94
 
                    $seqname = $1;
95
 
                    $start = $2;
96
 
                    $end = $3;
97
 
                } else {
98
 
                    $seqname=$name;
99
 
                    $start = 1;
100
 
                    $end = length($seqchar); #ps 9/6/00
101
 
                }
102
 
#               print STDERR  "Going to add with $seqchar $seqname\n";
103
 
                $seq = new Bio::LocatableSeq('-seq'        =>$seqchar,
104
 
                                             '-display_id' =>$seqname,
105
 
                                             '-description'=>$desc,
106
 
                                             '-start'      =>$start,
107
 
                                             '-end'        =>$end,
 
72
        my $self = shift;
 
73
        my ($width) = $self->_rearrange([qw(WIDTH)],@_);
 
74
        $self->width($width || $WIDTH);
 
75
 
 
76
        my ($start, $end, $name, $seqname, $seq, $seqchar, $entry, 
 
77
                 $tempname, $tempdesc, %align, $desc, $maxlen);
 
78
        my $aln = Bio::SimpleAlign->new();
 
79
 
 
80
        while (defined ($entry = $self->_readline) ) {
 
81
                chomp $entry;
 
82
                if ( $entry =~ s/^>\s*(\S+)\s*// ) {
 
83
                        $tempname  = $1;
 
84
                        chomp($entry);
 
85
                        $tempdesc = $entry;
 
86
                        if ( defined $name ) {
 
87
                                # put away last name and sequence
 
88
                                if ( $name =~ /(\S+)\/(\d+)-(\d+)/ ) {
 
89
                                        $seqname = $1;
 
90
                                        $start = $2;
 
91
                                        $end = $3;
 
92
                                } else {
 
93
                                        $seqname = $name;
 
94
                                        $start = 1;
 
95
                                        $end = $self->_get_len($seqchar);
 
96
                                }
 
97
                                $seq = new Bio::LocatableSeq(
 
98
                                                  -seq         => $seqchar,
 
99
                                             -display_id  => $seqname,
 
100
                                             -description => $desc,
 
101
                                             -start       => $start,
 
102
                                             -end         => $end,
108
103
                                             );
 
104
                                $aln->add_seq($seq);
 
105
                                $self->debug("Reading $seqname\n");
 
106
                        }
 
107
                        $desc = $tempdesc;      
 
108
                        $name = $tempname;
 
109
                        $desc = $entry;
 
110
                        $seqchar  = "";
 
111
                        next;
 
112
                }
 
113
                # removed redundant symbol validation
 
114
                # this is already done in Bio::PrimarySeq
 
115
                $seqchar .= $entry;
 
116
        }
 
117
 
 
118
        #  Next two lines are to silence warnings that
 
119
        #  otherwise occur at EOF when using <$fh>
 
120
        $name = "" if (!defined $name);
 
121
        $seqchar="" if (!defined $seqchar);
 
122
 
 
123
        #  Put away last name and sequence
 
124
        if ( $name =~ /(\S+)\/(\d+)-(\d+)/ ) {
 
125
                $seqname = $1;
 
126
                $start = $2;
 
127
                $end = $3;
 
128
        } else {
 
129
                $seqname = $name;
 
130
                $start = 1;
 
131
                $end = $self->_get_len($seqchar);
 
132
        }
 
133
 
 
134
        #  If $end <= 0, we have either reached the end of
 
135
        #  file in <> or we have encountered some other error
 
136
        if ( $end <= 0 ) { 
 
137
                undef $aln; 
 
138
                return $aln;
 
139
        }
 
140
 
 
141
        # This logic now also reads empty lines at the 
 
142
        # end of the file. Skip this is seqchar and seqname is null
 
143
        unless ( length($seqchar) == 0 && length($seqname) == 0 ) {
 
144
                $seq = new Bio::LocatableSeq(-seq         => $seqchar,
 
145
                                                                                          -display_id  => $seqname,
 
146
                                                                                          -description => $desc,
 
147
                                                                                          -start       => $start,
 
148
                                                                                          -end         => $end,
 
149
                                                                                         );
109
150
                $aln->add_seq($seq);
110
 
            }
111
 
            $desc = $tempdesc;  
112
 
            $name = $tempname;
113
 
            $desc = $entry;
114
 
            $seqchar  = "";
115
 
            next;
116
 
        }
117
 
        $entry =~ s/[^A-Za-z\.\-]//g;
118
 
        $seqchar .= $entry;     
119
 
    }
120
 
#
121
 
#  Next two lines are to silence warnings that
122
 
#  otherwise occur at EOF when using <$fh>
123
 
 
124
 
   if (!defined $name) {$name="";}
125
 
   if (!defined $seqchar) {$seqchar="";}
126
 
 
127
 
#  Put away last name and sequence
128
 
    if( $name =~ /(\S+)\/(\d+)-(\d+)/ ) {
129
 
        $seqname = $1;
130
 
        $start = $2;
131
 
        $end = $3;
132
 
    } else {
133
 
        $seqname=$name;
134
 
        $start = 1;
135
 
        $end = length($seqchar);   #ps 9/6/00
136
 
#       $end = length($align{$name});
137
 
    }
138
 
 
139
 
 
140
 
#  If $end <= 0, we have either reached the end of
141
 
#  file in <> or we have encountered some other error
142
 
#
143
 
   if ($end <= 0) { undef $aln; return $aln;}
144
 
 
145
 
# This logic now also reads empty lines at the 
146
 
# end of the file. Skip this is seqchar and seqname is null
147
 
 
148
 
    if( length($seqchar) == 0 && length($seqname) == 0 ) {
149
 
        # skip
150
 
    } else {
151
 
#       print STDERR "end to add with $seqchar $seqname\n";
152
 
        $seq = new Bio::LocatableSeq('-seq'        => $seqchar,
153
 
                                     '-display_id' => $seqname,
154
 
                                     '-description'=> $desc,
155
 
                                     '-start'      => $start,
156
 
                                     '-end'        => $end,
157
 
                                     );
158
 
        
159
 
        $aln->add_seq($seq);
160
 
    }
161
 
    my $alnlen = $aln->length;
162
 
    foreach my $seq ( $aln->each_seq ) {
163
 
        if( $seq->length < $alnlen ) {
164
 
            my ($diff) = ($alnlen - $seq->length);
165
 
            $seq->seq( $seq->seq() . "-" x $diff);
166
 
        }
167
 
    }
168
 
 
169
 
    return $aln;
170
 
 
 
151
                $self->debug("Reading $seqname\n");
 
152
        }
 
153
        my $alnlen = $aln->length;
 
154
        foreach my $seq ( $aln->each_seq ) {
 
155
                if ( $seq->length < $alnlen ) {
 
156
                        my ($diff) = ($alnlen - $seq->length);
 
157
                        $seq->seq( $seq->seq() . "-" x $diff);
 
158
                }
 
159
        }
 
160
        return $aln;
171
161
}
172
 
        
173
162
 
174
163
=head2 write_aln
175
164
 
179
168
 Returns : 1 for success and 0 for error
180
169
 Args    : L<Bio::Align::AlignI> object
181
170
 
 
171
See L<Bio::Align::AlignI>
182
172
 
183
173
=cut
184
174
 
185
175
sub write_aln {
186
176
    my ($self,@aln) = @_;
 
177
    my $width = $self->width;
187
178
    my ($seq,$desc,$rseq,$name,$count,$length,$seqsub);
188
179
 
189
180
    foreach my $aln (@aln) {
191
182
            $self->warn("Must provide a Bio::Align::AlignI object when calling write_aln");
192
183
            next;
193
184
        }
 
185
        if( $self->force_displayname_flat ) {
 
186
            $aln->set_displayname_flat(1);
 
187
        }
194
188
        foreach $rseq ( $aln->each_seq() ) {
195
189
            $name = $aln->displayname($rseq->get_nse());
196
190
            $seq  = $rseq->seq();
197
191
            $desc = $rseq->description || '';
198
192
            $self->_print (">$name $desc\n") or return ;        
199
 
            $count =0;
 
193
            $count = 0;
200
194
            $length = length($seq);
201
 
            while( ($count * 60 ) < $length ) {
202
 
                $seqsub = substr($seq,$count*60,60);
203
 
                $self->_print ("$seqsub\n") or return ;
204
 
                $count++;
 
195
            if(defined $seq && $length > 0) {
 
196
                $seq =~ s/(.{1,$width})/$1\n/g;
 
197
            } else {
 
198
                $seq = "\n";
205
199
            }
 
200
            $self->_print($seq);
206
201
        }
207
202
    }
208
203
    $self->flush if $self->_flush_on_write && defined $self->_fh;
209
204
    return 1;
210
205
}
211
206
 
 
207
=head2 _get_len
 
208
 
 
209
 Title   : _get_len
 
210
 Usage   : 
 
211
 Function: determine number of alphabetic chars
 
212
 Returns : integer
 
213
 Args    : sequence string
 
214
 
 
215
=cut
 
216
 
 
217
sub _get_len {
 
218
        my ($self,$seq) = @_;
 
219
        $seq =~ s/[^A-Z]//gi;
 
220
        return CORE::length($seq);
 
221
}
 
222
 
 
223
=head2 width
 
224
 
 
225
 Title   : width
 
226
 Usage   : $obj->width($newwidth)
 
227
           $width = $obj->width;
 
228
 Function: Get/set width of alignment
 
229
 Returns : integer value of width 
 
230
 Args    : on set, new value (a scalar or undef, optional)
 
231
 
 
232
 
 
233
=cut
 
234
 
 
235
sub width{
 
236
    my $self = shift;
 
237
 
 
238
    return $self->{'_width'} = shift if @_;
 
239
    return $self->{'_width'} || $WIDTH;
 
240
}
 
241
 
212
242
1;