|
0
|
1 # $Id: qual.pm,v 1.22 2002/12/27 19:42:32 birney Exp $
|
|
|
2 #
|
|
|
3 # Copyright (c) 1997-9 bioperl, Chad Matsalla. All Rights Reserved.
|
|
|
4 # This module is free software; you can redistribute it and/or
|
|
|
5 # modify it under the same terms as Perl itself.
|
|
|
6 #
|
|
|
7 # Copyright Chad Matsalla
|
|
|
8 #
|
|
|
9 # You may distribute this module under the same terms as perl itself
|
|
|
10
|
|
|
11 # POD documentation - main docs before the code
|
|
|
12
|
|
|
13 =head1 NAME
|
|
|
14
|
|
|
15 Bio::SeqIO::Qual - .qual file input/output stream
|
|
|
16
|
|
|
17 =head1 SYNOPSIS
|
|
|
18
|
|
|
19 Do not use this module directly. Use it via the Bio::SeqIO class
|
|
|
20 (see L<Bio::SeqIO> for details).
|
|
|
21
|
|
|
22 =head1 DESCRIPTION
|
|
|
23
|
|
|
24 This object can transform .qual (similar to fasta) objects to and from
|
|
|
25 Bio::Seq::SeqWithQuality objects. See L<Bio::Seq::SeqWithQuality> for
|
|
|
26 details.
|
|
|
27
|
|
|
28 =head1 FEEDBACK
|
|
|
29
|
|
|
30 =head2 Mailing Lists
|
|
|
31
|
|
|
32 User feedback is an integral part of the evolution of this and other
|
|
|
33 Bioperl modules. Send your comments and suggestions preferably to one
|
|
|
34 of the Bioperl mailing lists. Your participation is much appreciated.
|
|
|
35
|
|
|
36 bioperl-l@bioperl.org - General discussion
|
|
|
37 http://www.bioperl.org/MailList.shtml - About the mailing lists
|
|
|
38
|
|
|
39 =head2 Reporting Bugs
|
|
|
40
|
|
|
41 Report bugs to the Bioperl bug tracking system to help us keep track
|
|
|
42 the bugs and their resolution. Bug reports can be submitted via email
|
|
|
43 or the web:
|
|
|
44
|
|
|
45 bioperl-bugs@bio.perl.org
|
|
|
46 http://bugzilla.bioperl.org/
|
|
|
47
|
|
|
48 =head1 AUTHOR Chad Matsalla
|
|
|
49
|
|
|
50 Chad Matsalla
|
|
|
51 bioinformatics@dieselwurks.com
|
|
|
52
|
|
|
53 =head1 CONTRIBUTORS
|
|
|
54
|
|
|
55 Jason Stajich, jason@bioperl.org
|
|
|
56
|
|
|
57 =head1 APPENDIX
|
|
|
58
|
|
|
59 The rest of the documentation details each of the object
|
|
|
60 methods. Internal methods are usually preceded with a _
|
|
|
61
|
|
|
62 =cut
|
|
|
63
|
|
|
64 # Let the code begin...
|
|
|
65
|
|
|
66 package Bio::SeqIO::qual;
|
|
|
67 use vars qw(@ISA);
|
|
|
68 use strict;
|
|
|
69 use Bio::SeqIO;
|
|
|
70 use Bio::Seq::SeqFactory;
|
|
|
71 require 'dumpvar.pl';
|
|
|
72
|
|
|
73 @ISA = qw(Bio::SeqIO);
|
|
|
74
|
|
|
75
|
|
|
76 sub _initialize {
|
|
|
77 my($self,@args) = @_;
|
|
|
78 $self->SUPER::_initialize(@args);
|
|
|
79 if( ! defined $self->sequence_factory ) {
|
|
|
80 $self->sequence_factory(new Bio::Seq::SeqFactory
|
|
|
81 (-verbose => $self->verbose(),
|
|
|
82 -type => 'Bio::Seq::PrimaryQual'));
|
|
|
83 }
|
|
|
84 }
|
|
|
85
|
|
|
86 =head2 next_seq()
|
|
|
87
|
|
|
88 Title : next_seq()
|
|
|
89 Usage : $scf = $stream->next_seq()
|
|
|
90 Function: returns the next scf sequence in the stream
|
|
|
91 Returns : Bio::Seq::PrimaryQual object
|
|
|
92 Notes : Get the next quality sequence from the stream.
|
|
|
93
|
|
|
94 =cut
|
|
|
95
|
|
|
96 sub next_seq {
|
|
|
97 my ($self,@args) = @_;
|
|
|
98 my ($qual,$seq);
|
|
|
99 my $alphabet;
|
|
|
100 local $/ = "\n>";
|
|
|
101
|
|
|
102 return unless my $entry = $self->_readline;
|
|
|
103
|
|
|
104 if ($entry eq '>') { # very first one
|
|
|
105 return unless $entry = $self->_readline;
|
|
|
106 }
|
|
|
107
|
|
|
108 # original: my ($top,$sequence) = $entry =~ /^(.+?)\n([^>]*)/s
|
|
|
109 my ($top,$sequence) = $entry =~ /^(.+?)\n([^>]*)/s
|
|
|
110 or $self->throw("Can't parse entry [$entry]");
|
|
|
111 my ($id,$fulldesc) = $top =~ /^\s*(\S+)\s*(.*)/
|
|
|
112 or $self->throw("Can't parse fasta header");
|
|
|
113 $id =~ s/^>//;
|
|
|
114 # create the seq object
|
|
|
115 $sequence =~ s/\n+/ /g;
|
|
|
116 return $self->sequence_factory->create
|
|
|
117 (-qual => $sequence,
|
|
|
118 -id => $id,
|
|
|
119 -primary_id => $id,
|
|
|
120 -display_id => $id,
|
|
|
121 -desc => $fulldesc
|
|
|
122 );
|
|
|
123 }
|
|
|
124
|
|
|
125 =head2 _next_qual
|
|
|
126
|
|
|
127 Title : _next_qual
|
|
|
128 Usage : $seq = $stream->_next_qual() (but do not do
|
|
|
129 that. Use $stream->next_seq() instead)
|
|
|
130 Function: returns the next quality in the stream
|
|
|
131 Returns : Bio::Seq::PrimaryQual object
|
|
|
132 Args : NONE
|
|
|
133 Notes : An internal method. Gets the next quality in
|
|
|
134 the stream.
|
|
|
135
|
|
|
136 =cut
|
|
|
137
|
|
|
138 sub _next_qual {
|
|
|
139 my $qual = next_primary_qual( $_[0], 1 );
|
|
|
140 return $qual;
|
|
|
141 }
|
|
|
142
|
|
|
143 =head2 next_primary_qual()
|
|
|
144
|
|
|
145 Title : next_primary_qual()
|
|
|
146 Usage : $seq = $stream->next_primary_qual()
|
|
|
147 Function: returns the next sequence in the stream
|
|
|
148 Returns : Bio::PrimaryQual object
|
|
|
149 Args : NONE
|
|
|
150
|
|
|
151 =cut
|
|
|
152
|
|
|
153 sub next_primary_qual {
|
|
|
154 # print("CSM next_primary_qual!\n");
|
|
|
155 my( $self, $as_next_qual ) = @_;
|
|
|
156 my ($qual,$seq);
|
|
|
157 local $/ = "\n>";
|
|
|
158
|
|
|
159 return unless my $entry = $self->_readline;
|
|
|
160
|
|
|
161 if ($entry eq '>') { # very first one
|
|
|
162 return unless $entry = $self->_readline;
|
|
|
163 }
|
|
|
164
|
|
|
165 my ($top,$sequence) = $entry =~ /^(.+?)\n([^>]*)/s
|
|
|
166 or $self->throw("Can't parse entry [$entry]");
|
|
|
167 my ($id,$fulldesc) = $top =~ /^\s*(\S+)\s*(.*)/
|
|
|
168 or $self->throw("Can't parse fasta header");
|
|
|
169 $id =~ s/^>//;
|
|
|
170 # create the seq object
|
|
|
171 $sequence =~ s/\n+/ /g;
|
|
|
172 if ($as_next_qual) {
|
|
|
173 $qual = Bio::Seq::PrimaryQual->new(-qual => $sequence,
|
|
|
174 -id => $id,
|
|
|
175 -primary_id => $id,
|
|
|
176 -display_id => $id,
|
|
|
177 -desc => $fulldesc
|
|
|
178 );
|
|
|
179 }
|
|
|
180 return $qual;
|
|
|
181 }
|
|
|
182
|
|
|
183 =head2 write_seq
|
|
|
184
|
|
|
185 Title : write_seq(-source => $source, -header => "some information")
|
|
|
186 Usage : $obj->write_seq( -source => $source,
|
|
|
187 -header => "some information");
|
|
|
188 Function: Write out an list of quality values to a fasta-style file.
|
|
|
189 Returns : Nothing.
|
|
|
190 Args : Requires: a reference to a SeqWithQuality object or a
|
|
|
191 PrimaryQual object as the -source. Optional: information
|
|
|
192 for the header.
|
|
|
193 Notes : If no -header is provided, $obj->id() will be used where
|
|
|
194 $obj is a reference to either a SeqWithQuality object or a
|
|
|
195 PrimaryQual object. If $source->id() fails, ">unknown" will be
|
|
|
196 the header. If the SeqWithQuality object has $source->length() of
|
|
|
197 "DIFFERENT" (read the pod, luke), write_seq will use the length
|
|
|
198 of the PrimaryQual object within the SeqWithQuality object.
|
|
|
199
|
|
|
200 =cut
|
|
|
201
|
|
|
202 sub write_seq {
|
|
|
203 my ($self,@args) = @_;
|
|
|
204 my ($source) = $self->_rearrange([qw(SOURCE)], @args);
|
|
|
205
|
|
|
206 if (!$source || ( !$source->isa('Bio::Seq::SeqWithQuality') &&
|
|
|
207 !$source->isa('Bio::Seq::PrimaryQual') )) {
|
|
|
208 $self->throw("You must pass a Bio::Seq::SeqWithQuality or a Bio::Seq::PrimaryQual object to write_seq as a parameter named \"source\"");
|
|
|
209 }
|
|
|
210 my $header = $source->id();
|
|
|
211 if (!$header) { $header = "unknown"; }
|
|
|
212 my @quals = $source->qual();
|
|
|
213 # ::dumpValue(\@quals);
|
|
|
214 $self->_print (">$header \n");
|
|
|
215 my (@slice,$max,$length);
|
|
|
216 $length = $source->length();
|
|
|
217 if ($length eq "DIFFERENT") {
|
|
|
218 $self->warn("You passed a SeqWithQuality object that contains a sequence and quality of differing lengths. Using the length of the PrimaryQual component of the SeqWithQuality object.");
|
|
|
219 $length = $source->qual_obj()->length();
|
|
|
220 }
|
|
|
221 # print("Printing $header to a file.\n");
|
|
|
222 for (my $count = 1; $count<=$length; $count+= 50) {
|
|
|
223 if ($count+50 > $length) { $max = $length; }
|
|
|
224 else { $max = $count+49; }
|
|
|
225 my @slice = @{$source->subqual($count,$max)};
|
|
|
226 $self->_print (join(' ',@slice), "\n");
|
|
|
227 }
|
|
|
228
|
|
|
229 $self->flush if $self->_flush_on_write && defined $self->_fh;
|
|
|
230 return 1;
|
|
|
231 }
|
|
|
232
|
|
|
233
|
|
|
234 1;
|
|
|
235 __END__
|
