line |
stmt |
bran |
cond |
sub |
pod |
time |
code |
1
|
|
|
|
|
|
|
# |
2
|
|
|
|
|
|
|
# Copyright (c) 1997-9 bioperl, Chad Matsalla. All Rights Reserved. |
3
|
|
|
|
|
|
|
# This module is free software; you can redistribute it and/or |
4
|
|
|
|
|
|
|
# modify it under the same terms as Perl itself. |
5
|
|
|
|
|
|
|
# |
6
|
|
|
|
|
|
|
# Copyright Chad Matsalla |
7
|
|
|
|
|
|
|
# |
8
|
|
|
|
|
|
|
# You may distribute this module under the same terms as perl itself |
9
|
|
|
|
|
|
|
|
10
|
|
|
|
|
|
|
# POD documentation - main docs before the code |
11
|
|
|
|
|
|
|
|
12
|
|
|
|
|
|
|
=head1 NAME |
13
|
|
|
|
|
|
|
|
14
|
|
|
|
|
|
|
Bio::SeqIO::qual - .qual file input/output stream |
15
|
|
|
|
|
|
|
|
16
|
|
|
|
|
|
|
=head1 SYNOPSIS |
17
|
|
|
|
|
|
|
|
18
|
|
|
|
|
|
|
Do not use this module directly. Use it via the Bio::SeqIO class |
19
|
|
|
|
|
|
|
(see L for details). |
20
|
|
|
|
|
|
|
|
21
|
|
|
|
|
|
|
my $in_qual = Bio::SeqIO->new(-file => $qualfile, |
22
|
|
|
|
|
|
|
-format => 'qual', |
23
|
|
|
|
|
|
|
-width => $width, |
24
|
|
|
|
|
|
|
-verbose => $verbose); |
25
|
|
|
|
|
|
|
|
26
|
|
|
|
|
|
|
=head1 DESCRIPTION |
27
|
|
|
|
|
|
|
|
28
|
|
|
|
|
|
|
This object can transform .qual (similar to fasta) objects to and from |
29
|
|
|
|
|
|
|
Bio::Seq::Quality objects. See L for details. |
30
|
|
|
|
|
|
|
|
31
|
|
|
|
|
|
|
Like the fasta module, it can take an argument '-width' to change the |
32
|
|
|
|
|
|
|
number of values per line (defaults to 50). |
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 one |
40
|
|
|
|
|
|
|
of the Bioperl mailing lists. 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
|
|
|
|
|
|
|
the bugs and their resolution. Bug reports can be submitted via the web: |
60
|
|
|
|
|
|
|
|
61
|
|
|
|
|
|
|
https://github.com/bioperl/bioperl-live/issues |
62
|
|
|
|
|
|
|
|
63
|
|
|
|
|
|
|
=head1 AUTHOR Chad Matsalla |
64
|
|
|
|
|
|
|
|
65
|
|
|
|
|
|
|
Chad Matsalla |
66
|
|
|
|
|
|
|
bioinformatics@dieselwurks.com |
67
|
|
|
|
|
|
|
|
68
|
|
|
|
|
|
|
=head1 CONTRIBUTORS |
69
|
|
|
|
|
|
|
|
70
|
|
|
|
|
|
|
Jason Stajich, jason@bioperl.org |
71
|
|
|
|
|
|
|
|
72
|
|
|
|
|
|
|
=head1 APPENDIX |
73
|
|
|
|
|
|
|
|
74
|
|
|
|
|
|
|
The rest of the documentation details each of the object |
75
|
|
|
|
|
|
|
methods. Internal methods are usually preceded with a _ |
76
|
|
|
|
|
|
|
|
77
|
|
|
|
|
|
|
=cut |
78
|
|
|
|
|
|
|
|
79
|
|
|
|
|
|
|
# Let the code begin... |
80
|
|
|
|
|
|
|
|
81
|
|
|
|
|
|
|
package Bio::SeqIO::qual; |
82
|
2
|
|
|
2
|
|
14
|
use strict; |
|
2
|
|
|
|
|
4
|
|
|
2
|
|
|
|
|
77
|
|
83
|
2
|
|
|
2
|
|
493
|
use Bio::Seq::SeqFactory; |
|
2
|
|
|
|
|
4
|
|
|
2
|
|
|
|
|
49
|
|
84
|
2
|
|
|
2
|
|
630
|
use Dumpvalue; |
|
2
|
|
|
|
|
7425
|
|
|
2
|
|
|
|
|
76
|
|
85
|
|
|
|
|
|
|
|
86
|
|
|
|
|
|
|
my $dumper = Dumpvalue->new(); |
87
|
|
|
|
|
|
|
|
88
|
2
|
|
|
2
|
|
11
|
use base qw(Bio::SeqIO); |
|
2
|
|
|
|
|
3
|
|
|
2
|
|
|
|
|
1607
|
|
89
|
|
|
|
|
|
|
|
90
|
|
|
|
|
|
|
our $WIDTH = 25; |
91
|
|
|
|
|
|
|
|
92
|
|
|
|
|
|
|
sub _initialize { |
93
|
6
|
|
|
6
|
|
13
|
my($self,@args) = @_; |
94
|
6
|
|
|
|
|
21
|
$self->SUPER::_initialize(@args); |
95
|
6
|
|
|
|
|
20
|
my ($width) = $self->_rearrange([qw(WIDTH)], @args); |
96
|
6
|
50
|
|
|
|
15
|
$width && $self->width($width); |
97
|
6
|
50
|
|
|
|
22
|
if( ! defined $self->sequence_factory ) { |
98
|
6
|
|
|
|
|
16
|
$self->sequence_factory(Bio::Seq::SeqFactory->new |
99
|
|
|
|
|
|
|
(-verbose => $self->verbose(), |
100
|
|
|
|
|
|
|
-type => 'Bio::Seq::PrimaryQual')); |
101
|
|
|
|
|
|
|
} |
102
|
|
|
|
|
|
|
} |
103
|
|
|
|
|
|
|
|
104
|
|
|
|
|
|
|
=head2 next_seq() |
105
|
|
|
|
|
|
|
|
106
|
|
|
|
|
|
|
Title : next_seq() |
107
|
|
|
|
|
|
|
Usage : $scf = $stream->next_seq() |
108
|
|
|
|
|
|
|
Function: returns the next scf sequence in the stream |
109
|
|
|
|
|
|
|
Returns : Bio::Seq::PrimaryQual object |
110
|
|
|
|
|
|
|
Notes : Get the next quality sequence from the stream. |
111
|
|
|
|
|
|
|
|
112
|
|
|
|
|
|
|
=cut |
113
|
|
|
|
|
|
|
|
114
|
|
|
|
|
|
|
sub next_seq { |
115
|
19
|
|
|
19
|
1
|
551
|
my ($self,@args) = @_; |
116
|
19
|
|
|
|
|
30
|
my ($qual,$seq); |
117
|
19
|
|
|
|
|
0
|
my $alphabet; |
118
|
19
|
|
|
|
|
62
|
local $/ = "\n>"; |
119
|
|
|
|
|
|
|
|
120
|
19
|
100
|
|
|
|
55
|
return unless my $entry = $self->_readline; |
121
|
|
|
|
|
|
|
|
122
|
17
|
50
|
|
|
|
33
|
if ($entry eq '>') { # very first one |
123
|
0
|
0
|
|
|
|
0
|
return unless $entry = $self->_readline; |
124
|
|
|
|
|
|
|
} |
125
|
|
|
|
|
|
|
|
126
|
|
|
|
|
|
|
# original: my ($top,$sequence) = $entry =~ /^(.+?)\n([^>]*)/s |
127
|
17
|
50
|
|
|
|
316
|
my ($top,$sequence) = $entry =~ /^(.+?)\n([^>]*)/s |
128
|
|
|
|
|
|
|
or $self->throw("Can't parse entry [$entry]"); |
129
|
17
|
50
|
|
|
|
76
|
my ($id,$fulldesc) = $top =~ /^\s*(\S+)\s*(.*)/ |
130
|
|
|
|
|
|
|
or $self->throw("Can't parse fasta header"); |
131
|
17
|
|
|
|
|
30
|
$id =~ s/^>//; |
132
|
|
|
|
|
|
|
# create the seq object |
133
|
17
|
|
|
|
|
184
|
$sequence =~ s/\n+/ /g; |
134
|
17
|
|
|
|
|
45
|
return $self->sequence_factory->create |
135
|
|
|
|
|
|
|
(-qual => $sequence, |
136
|
|
|
|
|
|
|
-id => $id, |
137
|
|
|
|
|
|
|
-primary_id => $id, |
138
|
|
|
|
|
|
|
-display_id => $id, |
139
|
|
|
|
|
|
|
-desc => $fulldesc |
140
|
|
|
|
|
|
|
); |
141
|
|
|
|
|
|
|
} |
142
|
|
|
|
|
|
|
|
143
|
|
|
|
|
|
|
=head2 _next_qual |
144
|
|
|
|
|
|
|
|
145
|
|
|
|
|
|
|
Title : _next_qual |
146
|
|
|
|
|
|
|
Usage : $seq = $stream->_next_qual() (but do not do |
147
|
|
|
|
|
|
|
that. Use $stream->next_seq() instead) |
148
|
|
|
|
|
|
|
Function: returns the next quality in the stream |
149
|
|
|
|
|
|
|
Returns : Bio::Seq::PrimaryQual object |
150
|
|
|
|
|
|
|
Args : NONE |
151
|
|
|
|
|
|
|
Notes : An internal method. Gets the next quality in |
152
|
|
|
|
|
|
|
the stream. |
153
|
|
|
|
|
|
|
|
154
|
|
|
|
|
|
|
=cut |
155
|
|
|
|
|
|
|
|
156
|
|
|
|
|
|
|
sub _next_qual { |
157
|
0
|
|
|
0
|
|
0
|
my $qual = next_primary_qual( $_[0], 1 ); |
158
|
0
|
|
|
|
|
0
|
return $qual; |
159
|
|
|
|
|
|
|
} |
160
|
|
|
|
|
|
|
|
161
|
|
|
|
|
|
|
=head2 next_primary_qual() |
162
|
|
|
|
|
|
|
|
163
|
|
|
|
|
|
|
Title : next_primary_qual() |
164
|
|
|
|
|
|
|
Usage : $seq = $stream->next_primary_qual() |
165
|
|
|
|
|
|
|
Function: returns the next sequence in the stream |
166
|
|
|
|
|
|
|
Returns : Bio::PrimaryQual object |
167
|
|
|
|
|
|
|
Args : NONE |
168
|
|
|
|
|
|
|
|
169
|
|
|
|
|
|
|
=cut |
170
|
|
|
|
|
|
|
|
171
|
|
|
|
|
|
|
sub next_primary_qual { |
172
|
|
|
|
|
|
|
# print("CSM next_primary_qual!\n"); |
173
|
0
|
|
|
0
|
1
|
0
|
my( $self, $as_next_qual ) = @_; |
174
|
0
|
|
|
|
|
0
|
my ($qual,$seq); |
175
|
0
|
|
|
|
|
0
|
local $/ = "\n>"; |
176
|
|
|
|
|
|
|
|
177
|
0
|
0
|
|
|
|
0
|
return unless my $entry = $self->_readline; |
178
|
|
|
|
|
|
|
|
179
|
0
|
0
|
|
|
|
0
|
if ($entry eq '>') { # very first one |
180
|
0
|
0
|
|
|
|
0
|
return unless $entry = $self->_readline; |
181
|
|
|
|
|
|
|
} |
182
|
|
|
|
|
|
|
|
183
|
0
|
0
|
|
|
|
0
|
my ($top,$sequence) = $entry =~ /^(.+?)\n([^>]*)/s |
184
|
|
|
|
|
|
|
or $self->throw("Can't parse entry [$entry]"); |
185
|
0
|
0
|
|
|
|
0
|
my ($id,$fulldesc) = $top =~ /^\s*(\S+)\s*(.*)/ |
186
|
|
|
|
|
|
|
or $self->throw("Can't parse fasta header"); |
187
|
0
|
|
|
|
|
0
|
$id =~ s/^>//; |
188
|
|
|
|
|
|
|
# create the seq object |
189
|
0
|
|
|
|
|
0
|
$sequence =~ s/\n+/ /g; |
190
|
0
|
0
|
|
|
|
0
|
if ($as_next_qual) { |
191
|
0
|
|
|
|
|
0
|
$qual = Bio::Seq::PrimaryQual->new(-qual => $sequence, |
192
|
|
|
|
|
|
|
-id => $id, |
193
|
|
|
|
|
|
|
-primary_id => $id, |
194
|
|
|
|
|
|
|
-display_id => $id, |
195
|
|
|
|
|
|
|
-desc => $fulldesc |
196
|
|
|
|
|
|
|
); |
197
|
|
|
|
|
|
|
} |
198
|
0
|
|
|
|
|
0
|
return $qual; |
199
|
|
|
|
|
|
|
} |
200
|
|
|
|
|
|
|
|
201
|
|
|
|
|
|
|
|
202
|
|
|
|
|
|
|
=head2 width |
203
|
|
|
|
|
|
|
|
204
|
|
|
|
|
|
|
Title : width |
205
|
|
|
|
|
|
|
Usage : $obj->width($newval) |
206
|
|
|
|
|
|
|
Function: Get/Set the number of values per line for FASTA-like output |
207
|
|
|
|
|
|
|
Returns : value of width |
208
|
|
|
|
|
|
|
Args : newvalue (optional) |
209
|
|
|
|
|
|
|
|
210
|
|
|
|
|
|
|
|
211
|
|
|
|
|
|
|
=cut |
212
|
|
|
|
|
|
|
|
213
|
|
|
|
|
|
|
sub width{ |
214
|
12
|
|
|
12
|
1
|
20
|
my ($self,$value) = @_; |
215
|
12
|
50
|
|
|
|
20
|
if( defined $value) { |
216
|
0
|
|
|
|
|
0
|
$self->{'width'} = $value; |
217
|
|
|
|
|
|
|
} |
218
|
12
|
|
33
|
|
|
34
|
return $self->{'width'} || $WIDTH; |
219
|
|
|
|
|
|
|
} |
220
|
|
|
|
|
|
|
|
221
|
|
|
|
|
|
|
|
222
|
|
|
|
|
|
|
=head2 write_seq |
223
|
|
|
|
|
|
|
|
224
|
|
|
|
|
|
|
Title : write_seq |
225
|
|
|
|
|
|
|
Usage : $obj->write_seq( -source => $source, |
226
|
|
|
|
|
|
|
-header => "some information" |
227
|
|
|
|
|
|
|
-oneline => 0); |
228
|
|
|
|
|
|
|
Function: Write out a list of quality values to a fasta-style file. |
229
|
|
|
|
|
|
|
Returns : Nothing. |
230
|
|
|
|
|
|
|
Args : Requires a reference to a Bio::Seq::Quality object or a |
231
|
|
|
|
|
|
|
PrimaryQual object as the -source. Option 1: information |
232
|
|
|
|
|
|
|
for the header. Option 2: whether the quality score should |
233
|
|
|
|
|
|
|
be on a single line or not |
234
|
|
|
|
|
|
|
Notes : If no -header is provided, $obj->id() will be used where |
235
|
|
|
|
|
|
|
$obj is a reference to either a Quality object or a |
236
|
|
|
|
|
|
|
PrimaryQual object. If $source->id() fails, "unknown" will |
237
|
|
|
|
|
|
|
be the header. If the Quality object has $source->length() |
238
|
|
|
|
|
|
|
of "DIFFERENT" (read the pod, luke), write_seq will use the |
239
|
|
|
|
|
|
|
length of the PrimaryQual object within the Quality object. |
240
|
|
|
|
|
|
|
|
241
|
|
|
|
|
|
|
=cut |
242
|
|
|
|
|
|
|
|
243
|
|
|
|
|
|
|
sub write_seq { |
244
|
12
|
|
|
12
|
1
|
469
|
my ($self,@args) = @_; |
245
|
12
|
|
|
|
|
23
|
my $width = $self->width; |
246
|
12
|
|
|
|
|
41
|
my ($source, $head, $oneline) = $self->_rearrange([qw(SOURCE HEADER ONELINE)], @args); |
247
|
12
|
50
|
66
|
|
|
104
|
if (!$source || ( !$source->isa('Bio::Seq::Quality') && |
|
|
|
33
|
|
|
|
|
248
|
|
|
|
|
|
|
!$source->isa('Bio::Seq::PrimaryQual') )) { |
249
|
0
|
|
|
|
|
0
|
$self->throw("You must pass a Bio::Seq::Quality or a Bio::Seq::PrimaryQual". |
250
|
|
|
|
|
|
|
" object to write_seq() as a parameter named \"source\""); |
251
|
|
|
|
|
|
|
} |
252
|
12
|
100
|
100
|
|
|
77
|
my $header = ($source->can("header") && $source->header) ? |
|
|
100
|
66
|
|
|
|
|
253
|
|
|
|
|
|
|
$source->header : |
254
|
|
|
|
|
|
|
($source->can("id") && $source->id) ? |
255
|
|
|
|
|
|
|
$source->id : |
256
|
|
|
|
|
|
|
"unknown"; |
257
|
12
|
|
|
|
|
26
|
my @quals = $source->qual(); |
258
|
|
|
|
|
|
|
# ::dumpValue(\@quals); |
259
|
12
|
50
|
|
|
|
51
|
my $desc = $source->desc if $source->can('desc'); |
260
|
12
|
|
50
|
|
|
47
|
$desc ||= ''; |
261
|
12
|
|
|
|
|
55
|
$self->_print (">$header $desc\n"); |
262
|
12
|
|
|
|
|
16
|
my (@slice,$max,$length); |
263
|
12
|
|
|
|
|
30
|
$length = $source->length(); |
264
|
|
|
|
|
|
|
|
265
|
12
|
50
|
33
|
|
|
34
|
if ( not(defined($oneline)) || $oneline == 0) { |
266
|
|
|
|
|
|
|
# $width quality values per line |
267
|
12
|
|
|
|
|
26
|
for (my $count = 1; $count<=$length; $count+= $width) { |
268
|
235
|
100
|
|
|
|
316
|
if ($count+$width > $length) { |
269
|
12
|
|
|
|
|
12
|
$max = $length; |
270
|
|
|
|
|
|
|
} else { |
271
|
223
|
|
|
|
|
218
|
$max = $count+$width-1; |
272
|
|
|
|
|
|
|
} |
273
|
235
|
|
|
|
|
217
|
my @slice = @{$source->subqual($count,$max)}; |
|
235
|
|
|
|
|
342
|
|
274
|
235
|
|
|
|
|
940
|
$self->_print (join(' ',@slice), "\n"); |
275
|
|
|
|
|
|
|
} |
276
|
|
|
|
|
|
|
} else { |
277
|
|
|
|
|
|
|
# quality values on a single line |
278
|
0
|
|
|
|
|
0
|
my @slice = @{$source->qual}; |
|
0
|
|
|
|
|
0
|
|
279
|
0
|
|
|
|
|
0
|
$self->_print (join(' ',@slice), "\n"); |
280
|
|
|
|
|
|
|
} |
281
|
|
|
|
|
|
|
|
282
|
12
|
50
|
33
|
|
|
28
|
$self->flush if $self->_flush_on_write && defined $self->_fh; |
283
|
12
|
|
|
|
|
50
|
return 1; |
284
|
|
|
|
|
|
|
} |
285
|
|
|
|
|
|
|
|
286
|
|
|
|
|
|
|
|
287
|
|
|
|
|
|
|
1; |
288
|
|
|
|
|
|
|
|
289
|
|
|
|
|
|
|
__END__ |