File Coverage

Bio/AlignIO/emboss.pm
Criterion Covered Total %
statement 84 88 95.4
branch 41 44 93.1
condition 23 26 88.4
subroutine 7 8 87.5
pod 2 2 100.0
total 157 168 93.4


line stmt bran cond sub pod time code
1             #
2             # BioPerl module for Bio::AlignIO::emboss
3             #
4             # Please direct questions and support issues to
5             #
6             # Cared for by Jason Stajich
7             #
8             # Copyright Jason Stajich
9             #
10             # You may distribute this module under the same terms as perl itself
11              
12             # POD documentation - main docs before the code
13              
14             =head1 NAME
15              
16             Bio::AlignIO::emboss - Parse EMBOSS alignment output (from applications water and needle)
17              
18             =head1 SYNOPSIS
19              
20             # do not use the object directly
21             use Bio::AlignIO;
22             # read in an alignment from the EMBOSS program water
23             my $in = Bio::AlignIO->new(-format => 'emboss',
24             -file => 'seq.water');
25             while( my $aln = $in->next_aln ) {
26             # do something with the alignment
27             }
28              
29             =head1 DESCRIPTION
30              
31             This object handles parsing and writing pairwise sequence alignments
32             from the EMBOSS suite.
33              
34             =head1 FEEDBACK
35              
36             =head2 Mailing Lists
37              
38             User feedback is an integral part of the evolution of this and other
39             Bioperl modules. Send your comments and suggestions preferably to
40             the Bioperl mailing list. Your participation is much appreciated.
41              
42             bioperl-l@bioperl.org - General discussion
43             http://bioperl.org/wiki/Mailing_lists - About the mailing lists
44              
45             =head2 Support
46              
47             Please direct usage questions or support issues to the mailing list:
48              
49             I
50              
51             rather than to the module maintainer directly. Many experienced and
52             reponsive experts will be able look at the problem and quickly
53             address it. Please include a thorough description of the problem
54             with code and data examples if at all possible.
55              
56             =head2 Reporting Bugs
57              
58             Report bugs to the Bioperl bug tracking system to help us keep track
59             of the bugs and their resolution. Bug reports can be submitted via the
60             web:
61              
62             https://github.com/bioperl/bioperl-live/issues
63              
64             =head1 AUTHOR - Jason Stajich
65              
66             Email jason@bioperl.org
67              
68             =head1 APPENDIX
69              
70             The rest of the documentation details each of the object methods.
71             Internal methods are usually preceded with a _
72              
73             =cut
74              
75              
76             # Let the code begin...
77              
78              
79             package Bio::AlignIO::emboss;
80 2     2   350 use vars qw($EMBOSSTitleLen $EMBOSSLineLen);
  2         3  
  2         98  
81 2     2   10 use strict;
  2         3  
  2         35  
82              
83 2     2   294 use Bio::LocatableSeq;
  2         102  
  2         56  
84              
85 2     2   56 use base qw(Bio::AlignIO);
  2         2  
  2         371  
86              
87             BEGIN {
88 2     2   2 $EMBOSSTitleLen = 13;
89 2         1420 $EMBOSSLineLen = 50;
90             }
91              
92             sub _initialize {
93 6     6   8 my($self,@args) = @_;
94 6         17 $self->SUPER::_initialize(@args);
95 6         12 $self->{'_type'} = undef;
96             }
97              
98             =head2 next_aln
99              
100             Title : next_aln
101             Usage : $aln = $stream->next_aln()
102             Function: returns the next alignment in the stream.
103             Returns : L object - returns 0 on end of file
104             or on error
105             Args : NONE
106              
107             =cut
108              
109             sub next_aln {
110 7     7 1 12 my ($self) = @_;
111 7         5 my $seenbegin = 0;
112             my %data = ( 'seq1' => {
113             'start'=> undef,
114             'end'=> undef,
115             'name' => '',
116             'data' => '' },
117             'seq2' => {
118             'start'=> undef,
119             'end'=> undef,
120             'name' => '',
121             'data' => '' },
122             'align' => '',
123 7         36 'type' => $self->{'_type'}, # to restore type from
124             # previous aln if possible
125             );
126 7         8 my %names;
127 7         16 while( defined($_ = $self->_readline) ) {
128 226 100 100     833 next if( /^\#?\s+$/ || /^\#+\s*$/ );
129 139 100 100     1149 if( /^\#(\=|\-)+\s*$/) {
    100 66        
    100 100        
    100          
    100          
130 13 100       28 last if( $seenbegin);
131             } elsif( /(Local|Global):\s*(\S+)\s+vs\s+(\S+)/ ||
132             /^\#\s+Program:\s+(\S+)/ )
133             {
134 6         14 my ($name1,$name2) = ($2,$3);
135 6 100       9 if( ! defined $name1 ) { # Handle EMBOSS 2.2.X
136 4         7 $data{'type'} = $1;
137 4         5 $name1 = $name2 = '';
138             } else {
139 2 100       5 $data{'type'} = $1 eq 'Local' ? 'water' : 'needle';
140             }
141 6         10 $data{'seq1'}->{'name'} = $name1;
142 6         6 $data{'seq2'}->{'name'} = $name2;
143              
144 6         13 $self->{'_type'} = $data{'type'};
145              
146             } elsif( /Score:\s+(\S+)/ ) {
147 7         18 $data{'score'} = $1;
148             } elsif( /^\#\s+(1|2):\s+(\S+)/ && ! $data{"seq$1"}->{'name'} ) {
149 9         10 my $nm = $2;
150 9         13 $nm = substr($nm,0,$EMBOSSTitleLen); # emboss has a max seq length
151 9 100       13 if( $names{$nm} ) {
152 1         2 $nm .= "-". $names{$nm};
153             }
154 9         12 $names{$nm}++;
155 9         21 $data{"seq$1"}->{'name'} = $nm;
156             } elsif( $data{'seq1'}->{'name'} &&
157             /^\Q$data{'seq1'}->{'name'}/ ) {
158 51         39 my $count = 0;
159 51         38 $seenbegin = 1;
160 51         31 my @current;
161 51         67 while( defined ($_) ) {
162 153         113 my $align_other = '';
163 153         78 my $delayed;
164 153 100 100     320 if($count == 0 || $count == 2 ) {
165 102         180 my @l = split;
166 102         65 my ($seq,$align,$start,$end);
167 102 100 100     302 if( $count == 2 && $data{'seq2'}->{'name'} eq '' ) {
    100          
168             # weird boundary condition
169 4         5 ($start,$align,$end) = @l;
170             } elsif( @l == 3 ) {
171 4         4 $align = '';
172 4         5 ($seq,$start,$end) = @l
173             } else {
174 94         107 ($seq,$start,$align,$end) = @l;
175             }
176              
177 102 100       199 my $seqname = sprintf("seq%d", ($count == 0) ? '1' : '2');
178 102         127 $data{$seqname}->{'data'} .= $align;
179 102   66     151 $data{$seqname}->{'start'} ||= $start;
180 102         96 $data{$seqname}->{'end'} = $end;
181 102   100     245 $current[$count] = [ $start,$align || ''];
182             } else {
183 51         104 s/^\s+//;
184 51         215 s/\s+$//;
185 51         70 $data{'align'} .= $_;
186             }
187              
188 153 100       227 BOTTOM:
189             last if( $count++ == 2);
190 102         127 $_ = $self->_readline();
191             }
192              
193 51 100       107 if( $data{'type'} eq 'needle' ) {
194             # which ever one is shorter we want to bring it up to
195             # length. Man this stinks.
196 24         21 my ($s1,$s2) = ($data{'seq1'}, $data{'seq2'});
197              
198 24         23 my $d = length($current[0]->[1]) - length($current[2]->[1]);
199 24 100       65 if( $d < 0 ) { # s1 is smaller, need to add some
    100          
200             # compare the starting points for this alignment line
201 3 50       6 if( $current[0]->[0] <= 1 ) {
202 0         0 $s1->{'data'} = ('-' x abs($d)) . $s1->{'data'};
203 0         0 $data{'align'} = (' 'x abs($d)).$data{'align'};
204             } else {
205 3         5 $s1->{'data'} .= '-' x abs($d);
206 3         11 $data{'align'} .= ' 'x abs($d);
207             }
208             } elsif( $d > 0) { # s2 is smaller, need to add some
209 6 100       291 if( $current[2]->[0] <= 1 ) {
210 3         10 $s2->{'data'} = ('-' x abs($d)) . $s2->{'data'};
211 3         11 $data{'align'} = (' 'x abs($d)).$data{'align'};
212             } else {
213 3         6 $s2->{'data'} .= '-' x abs($d);
214 3         10 $data{'align'} .= ' 'x abs($d);
215             }
216             }
217             }
218              
219             }
220             }
221 7 50       11 return unless $seenbegin;
222             my $aln = Bio::SimpleAlign->new(-verbose => $self->verbose(),
223             -score => $data{'score'},
224 7         39 -source => "EMBOSS-".$data{'type'});
225              
226 7         12 foreach my $seqname ( qw(seq1 seq2) ) {
227 14 50       22 return unless ( defined $data{$seqname} );
228 14   66     27 $data{$seqname}->{'name'} ||= $seqname;
229             my $seq = Bio::LocatableSeq->new
230             ('-seq' => $data{$seqname}->{'data'},
231             '-display_id' => $data{$seqname}->{'name'},
232             '-start' => $data{$seqname}->{'start'},
233 14         41 '-end' => $data{$seqname}->{'end'},
234             '-alphabet' => $self->alphabet,
235             );
236 14         30 $aln->add_seq($seq);
237             }
238 7         35 return $aln;
239             }
240              
241             =head2 write_aln
242              
243             Title : write_aln
244             Usage : $stream->write_aln(@aln)
245             Function: writes the $aln object into the stream in emboss format
246             Returns : 1 for success and 0 for error
247             Args : L object
248              
249              
250             =cut
251              
252             sub write_aln {
253 0     0 1   my ($self,@aln) = @_;
254              
255 0           $self->throw("Sorry: writing emboss output is not currently available! \n");
256             }
257              
258             1;