|
0
|
1 # $Id: SeqWithQuality.pm,v 1.17 2002/12/19 22:02:38 matsallac Exp $
|
|
|
2 #
|
|
|
3 # BioPerl module for Bio::Seq::QualI
|
|
|
4 #
|
|
|
5 # Cared for by Chad Matsalla <bioinformatics@dieselwurks.com
|
|
|
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::Seq::SeqWithQuality - Bioperl object packaging a sequence with its quality
|
|
|
16
|
|
|
17 =head1 SYNOPSIS
|
|
|
18
|
|
|
19 use Bio::PrimarySeq;
|
|
|
20 use Bio::Seq::PrimaryQual;
|
|
|
21
|
|
|
22 # make from memory
|
|
|
23 my $qual = Bio::Seq::SeqWithQuality->new
|
|
|
24 ( -qual => '10 20 30 40 50 50 20 10',
|
|
|
25 -seq => 'ATCGATCG',
|
|
|
26 -id => 'human_id',
|
|
|
27 -accession_number => 'AL000012',
|
|
|
28 );
|
|
|
29
|
|
|
30 # make from objects
|
|
|
31 # first, make a PrimarySeq object
|
|
|
32 my $seqobj = Bio::PrimarySeq->new
|
|
|
33 ( -seq => 'atcgatcg',
|
|
|
34 -id => 'GeneFragment-12',
|
|
|
35 -accession_number => 'X78121',
|
|
|
36 -alphabet => 'dna'
|
|
|
37 );
|
|
|
38
|
|
|
39 # now make a PrimaryQual object
|
|
|
40 my $qualobj = Bio::Seq::PrimaryQual->new
|
|
|
41 ( -qual => '10 20 30 40 50 50 20 10',
|
|
|
42 -id => 'GeneFragment-12',
|
|
|
43 -accession_number => 'X78121',
|
|
|
44 -alphabet => 'dna'
|
|
|
45 );
|
|
|
46
|
|
|
47 # now make the SeqWithQuality object
|
|
|
48 my $swqobj = Bio::Seq::SeqWithQuality->new
|
|
|
49 ( -seq => $seqobj,
|
|
|
50 -qual => $qualobj
|
|
|
51 );
|
|
|
52 # done!
|
|
|
53
|
|
|
54 $swqobj->id(); # the id of the SeqWithQuality object
|
|
|
55 # may not match the the id of the sequence or
|
|
|
56 # of the quality (check the pod, luke)
|
|
|
57 $swqobj->seq(); # the sequence of the SeqWithQuality object
|
|
|
58 $swqobj->qual(); # the quality of the SeqWithQuality object
|
|
|
59
|
|
|
60 # to get out parts of the sequence.
|
|
|
61
|
|
|
62 print "Sequence ", $seqobj->id(), " with accession ",
|
|
|
63 $seqobj->accession, " and desc ", $seqobj->desc, "\n";
|
|
|
64
|
|
|
65 $string2 = $seqobj->subseq(1,40);
|
|
|
66
|
|
|
67 =head1 DESCRIPTION
|
|
|
68
|
|
|
69 This object stores base quality values together with the sequence string.
|
|
|
70
|
|
|
71 =head1 FEEDBACK
|
|
|
72
|
|
|
73 =head2 Mailing Lists
|
|
|
74
|
|
|
75 User feedback is an integral part of the evolution of this and other
|
|
|
76 Bioperl modules. Send your comments and suggestions preferably to one
|
|
|
77 of the Bioperl mailing lists. Your participation is much appreciated.
|
|
|
78
|
|
|
79 bioperl-l@bioperl.org - General discussion
|
|
|
80 http://bio.perl.org/MailList.html - About the mailing lists
|
|
|
81
|
|
|
82 =head2 Reporting Bugs
|
|
|
83
|
|
|
84 Report bugs to the Bioperl bug tracking system to help us keep track
|
|
|
85 the bugs and their resolution. Bug reports can be submitted via email
|
|
|
86 or the web:
|
|
|
87
|
|
|
88 bioperl-bugs@bio.perl.org
|
|
|
89 http://bugzilla.bioperl.org/
|
|
|
90
|
|
|
91 =head1 AUTHOR - Chad Matsalla
|
|
|
92
|
|
|
93 Email bioinformatics@dieselwurks.com
|
|
|
94
|
|
|
95 =head1 CONTRIBUTORS
|
|
|
96
|
|
|
97 Jason Stajich, jason@bioperl.org
|
|
|
98
|
|
|
99 =head1 APPENDIX
|
|
|
100
|
|
|
101 The rest of the documentation details each of the object methods.
|
|
|
102 Internal methods are usually preceded with a _
|
|
|
103
|
|
|
104 =cut
|
|
|
105
|
|
|
106
|
|
|
107 package Bio::Seq::SeqWithQuality;
|
|
|
108
|
|
|
109 use vars qw(@ISA);
|
|
|
110
|
|
|
111 use strict;
|
|
|
112 use Bio::Root::Root;
|
|
|
113 use Bio::Seq::QualI;
|
|
|
114 use Bio::PrimarySeqI;
|
|
|
115 use Bio::PrimarySeq;
|
|
|
116 use Bio::Seq::PrimaryQual;
|
|
|
117
|
|
|
118 @ISA = qw(Bio::Root::Root Bio::PrimarySeqI Bio::Seq::QualI);
|
|
|
119
|
|
|
120 =head2 new()
|
|
|
121
|
|
|
122 Title : new()
|
|
|
123 Usage : $qual = Bio::Seq::SeqWithQuality ->new
|
|
|
124 ( -qual => '10 20 30 40 50 50 20 10',
|
|
|
125 -seq => 'ATCGATCG',
|
|
|
126 -id => 'human_id',
|
|
|
127 -accession_number => 'AL000012',
|
|
|
128 -trace_indices => '0 5 10 15 20 25 30 35'
|
|
|
129 );
|
|
|
130 Function: Returns a new Bio::Seq::SeqWithQual object from basic
|
|
|
131 constructors.
|
|
|
132 Returns : a new Bio::Seq::PrimaryQual object
|
|
|
133 Notes : Arguments:
|
|
|
134 -qual can be a quality string (see Bio::Seq::PrimaryQual for more
|
|
|
135 information on this) or a reference to a Bio::Seq::PrimaryQual
|
|
|
136 object.
|
|
|
137 -seq can be a sequence string (see Bio::PrimarySeq for more
|
|
|
138 information on this) or a reference to a Bio::PrimaryQual object.
|
|
|
139 -seq, -id, -accession_number, -primary_id, -desc, -id behave like
|
|
|
140 this:
|
|
|
141 1. if they are provided on construction of the
|
|
|
142 Bio::Seq::SeqWithQuality they will be set as the descriptors for
|
|
|
143 the object unless changed by one of the following mechanisms:
|
|
|
144 a) $obj->set_common_descriptors() is used and both the -seq and
|
|
|
145 the -qual object have the same descriptors. These common
|
|
|
146 descriptors will then become the descriptors for the
|
|
|
147 Bio::Seq::SeqWithQual object.
|
|
|
148 b) the descriptors are manually set using the seq(), id(),
|
|
|
149 desc(), or accession_number(), primary_id(),
|
|
|
150 2. if no descriptors are provided, the new() constructor will see
|
|
|
151 if the descriptor used in the PrimarySeq and in the
|
|
|
152 PrimaryQual objects match. If they do, they will become
|
|
|
153 the descriptors for the SeqWithQuality object.
|
|
|
154
|
|
|
155 To eliminate ambiguity, I strongly suggest you set the
|
|
|
156 descriptors manually on construction of the object. Really.
|
|
|
157 -trace_indices : a space_delimited list of trace indices
|
|
|
158 (where would the peaks be drawn if this list of qualities
|
|
|
159 was to be plotted?)
|
|
|
160
|
|
|
161 =cut
|
|
|
162
|
|
|
163 sub new {
|
|
|
164 my ($class, @args) = @_;
|
|
|
165 my $self = $class->SUPER::new(@args);
|
|
|
166 # default: turn OFF the warnings
|
|
|
167 $self->{supress_warnings} = 1;
|
|
|
168 my($qual,$seq,$id,$acc,$pid,$desc,$given_id,$alphabet,$trace_indices) =
|
|
|
169 $self->_rearrange([qw(
|
|
|
170 QUAL
|
|
|
171 SEQ
|
|
|
172 DISPLAY_ID
|
|
|
173 ACCESSION_NUMBER
|
|
|
174 PRIMARY_ID
|
|
|
175 DESC
|
|
|
176 ID
|
|
|
177 ALPHABET
|
|
|
178 TRACE_INDICES
|
|
|
179 )],
|
|
|
180 @args);
|
|
|
181 # first, deal with the sequence and quality information
|
|
|
182 if ( defined $id && defined $given_id ) {
|
|
|
183 if( $id ne $given_id ) {
|
|
|
184 $self->throw("Provided both id and display_id constructor functions. [$id] [$given_id]");
|
|
|
185 }
|
|
|
186 }
|
|
|
187 if( defined $given_id ) {
|
|
|
188 $self->display_id($given_id);
|
|
|
189 $id = $given_id;
|
|
|
190 }
|
|
|
191 if (!$seq) {
|
|
|
192 my $id;
|
|
|
193 unless ($self->{supress_warnings} == 1) {
|
|
|
194 $self->warn("You did not provide sequence information during the construction of a Bio::Seq::SeqWithQuality object. Sequence components for this object will be empty.");
|
|
|
195 }
|
|
|
196 if (!$alphabet) {
|
|
|
197 $self->throw("If you want me to create a PrimarySeq object for your empty sequence <boggle> you must specify a -alphabet to satisfy the constructor requirements for a Bio::PrimarySeq object with no sequence. Read the POD for it, luke.");
|
|
|
198 }
|
|
|
199 $self->{seq_ref} = Bio::PrimarySeq->new
|
|
|
200 (
|
|
|
201 -seq => "",
|
|
|
202 -accession_number => $acc,
|
|
|
203 -primary_id => $pid,
|
|
|
204 -desc => $desc,
|
|
|
205 -display_id => $id,
|
|
|
206 -alphabet => $alphabet
|
|
|
207 );
|
|
|
208 }
|
|
|
209 elsif (ref($seq) eq "Bio::PrimarySeq" ) {
|
|
|
210 $self->{seq_ref} = $seq;
|
|
|
211 }
|
|
|
212
|
|
|
213 else {
|
|
|
214 my $seqobj = Bio::PrimarySeq->new
|
|
|
215 (
|
|
|
216 -seq => $seq,
|
|
|
217 -accession_number => $acc,
|
|
|
218 -primary_id => $pid,
|
|
|
219 -desc => $desc,
|
|
|
220 -display_id => $id,
|
|
|
221 );
|
|
|
222 $self->{seq_ref} = $seqobj;
|
|
|
223 }
|
|
|
224
|
|
|
225 if (!$qual) {
|
|
|
226 $self->{qual_ref} = Bio::Seq::PrimaryQual->new
|
|
|
227 (
|
|
|
228 -qual => "",
|
|
|
229 -accession_number => $acc,
|
|
|
230 -primary_id => $pid,
|
|
|
231 -desc => $desc,
|
|
|
232 -display_id => $id,
|
|
|
233 );
|
|
|
234 }
|
|
|
235 elsif (ref($qual) eq "Bio::Seq::PrimaryQual") {
|
|
|
236 $self->{qual_ref} = $qual;
|
|
|
237 }
|
|
|
238 else {
|
|
|
239 my $qualobj = Bio::Seq::PrimaryQual->new
|
|
|
240 (
|
|
|
241 -qual => $qual,
|
|
|
242 -accession_number => $acc,
|
|
|
243 -primary_id => $pid,
|
|
|
244 -desc => $desc,
|
|
|
245 -display_id => $id,
|
|
|
246 -trace_indices => $trace_indices
|
|
|
247 );
|
|
|
248 $self->{qual_ref} = $qualobj;
|
|
|
249 }
|
|
|
250
|
|
|
251 # now try to set the descriptors for this object
|
|
|
252 $self->_set_descriptors($qual,$seq,$id,$acc,$pid,$desc,$given_id,$alphabet);
|
|
|
253 $self->length();
|
|
|
254
|
|
|
255 return $self;
|
|
|
256 }
|
|
|
257
|
|
|
258 =head2 _common_id()
|
|
|
259
|
|
|
260 Title : _common_id()
|
|
|
261 Usage : $common_id = $self->_common_id();
|
|
|
262 Function: Compare the display_id of {qual_ref} and {seq_ref}.
|
|
|
263 Returns : Nothing if they don't match. If they do return
|
|
|
264 {seq_ref}->display_id()
|
|
|
265 Args : None.
|
|
|
266
|
|
|
267 =cut
|
|
|
268
|
|
|
269 #'
|
|
|
270 sub _common_id {
|
|
|
271 my $self = shift;
|
|
|
272 return if (!$self->{seq_ref} || !$self->{qual_ref});
|
|
|
273 my $sid = $self->{seq_ref}->display_id();
|
|
|
274 return if (!$sid);
|
|
|
275 return if (!$self->{qual_ref}->display_id());
|
|
|
276 return $sid if ($sid eq $self->{qual_ref}->display_id());
|
|
|
277 # should this become a warning?
|
|
|
278 # print("ids $sid and $self->{qual_ref}->display_id() do not match. Bummer.\n");
|
|
|
279 }
|
|
|
280
|
|
|
281 =head2 _common_display_id()
|
|
|
282
|
|
|
283 Title : _common_id()
|
|
|
284 Usage : $common_id = $self->_common_display_id();
|
|
|
285 Function: Compare the display_id of {qual_ref} and {seq_ref}.
|
|
|
286 Returns : Nothing if they don't match. If they do return
|
|
|
287 {seq_ref}->display_id()
|
|
|
288 Args : None.
|
|
|
289
|
|
|
290 =cut
|
|
|
291
|
|
|
292 #'
|
|
|
293 sub _common_display_id {
|
|
|
294 my $self = shift;
|
|
|
295 $self->common_id();
|
|
|
296 }
|
|
|
297
|
|
|
298 =head2 _common_accession_number()
|
|
|
299
|
|
|
300 Title : _common_accession_number()
|
|
|
301 Usage : $common_id = $self->_common_accession_number();
|
|
|
302 Function: Compare the accession_number() of {qual_ref} and {seq_ref}.
|
|
|
303 Returns : Nothing if they don't match. If they do return
|
|
|
304 {seq_ref}->accession_number()
|
|
|
305 Args : None.
|
|
|
306
|
|
|
307 =cut
|
|
|
308
|
|
|
309 #'
|
|
|
310 sub _common_accession_number {
|
|
|
311 my $self = shift;
|
|
|
312 return if ($self->{seq_ref} || $self->{qual_ref});
|
|
|
313 my $acc = $self->{seq_ref}->accession_number();
|
|
|
314 # if (!$acc) { print("the seqref has no acc.\n"); }
|
|
|
315 return if (!$acc);
|
|
|
316 # if ($acc eq $self->{qual_ref}->accession_number()) { print("$acc matches ".$self->{qual_ref}->accession_number()."\n"); }
|
|
|
317 return $acc if ($acc eq $self->{qual_ref}->accession_number());
|
|
|
318 # should this become a warning?
|
|
|
319 # print("accession numbers $acc and $self->{qual_ref}->accession_number() do not match. Bummer.\n");
|
|
|
320 }
|
|
|
321
|
|
|
322 =head2 _common_primary_id()
|
|
|
323
|
|
|
324 Title : _common_primary_id()
|
|
|
325 Usage : $common_primard_id = $self->_common_primary_id();
|
|
|
326 Function: Compare the primary_id of {qual_ref} and {seq_ref}.
|
|
|
327 Returns : Nothing if they don't match. If they do return
|
|
|
328 {seq_ref}->primary_id()
|
|
|
329 Args : None.
|
|
|
330
|
|
|
331 =cut
|
|
|
332
|
|
|
333 #'
|
|
|
334 sub _common_primary_id {
|
|
|
335 my $self = shift;
|
|
|
336 return if ($self->{seq_ref} || $self->{qual_ref});
|
|
|
337 my $pid = $self->{seq_ref}->primary_id();
|
|
|
338 return if (!$pid);
|
|
|
339 return $pid if ($pid eq $self->{qual_ref}->primary_id());
|
|
|
340 # should this become a warning?
|
|
|
341 # print("primary_ids $pid and $self->{qual_ref}->primary_id() do not match. Bummer.\n");
|
|
|
342
|
|
|
343 }
|
|
|
344
|
|
|
345 =head2 _common_desc()
|
|
|
346
|
|
|
347 Title : _common_desc()
|
|
|
348 Usage : $common_desc = $self->_common_desc();
|
|
|
349 Function: Compare the desc of {qual_ref} and {seq_ref}.
|
|
|
350 Returns : Nothing if they don't match. If they do return
|
|
|
351 {seq_ref}->desc()
|
|
|
352 Args : None.
|
|
|
353
|
|
|
354 =cut
|
|
|
355
|
|
|
356 #'
|
|
|
357 sub _common_desc {
|
|
|
358 my $self = shift;
|
|
|
359 return if ($self->{seq_ref} || $self->{qual_ref});
|
|
|
360 my $des = $self->{seq_ref}->desc();
|
|
|
361 return if (!$des);
|
|
|
362 return $des if ($des eq $self->{qual_ref}->desc());
|
|
|
363 # should this become a warning?
|
|
|
364 # print("descriptions $des and $self->{qual_ref}->desc() do not match. Bummer.\n");
|
|
|
365
|
|
|
366 }
|
|
|
367
|
|
|
368 =head2 set_common_descriptors()
|
|
|
369
|
|
|
370 Title : set_common_descriptors()
|
|
|
371 Usage : $self->set_common_descriptors();
|
|
|
372 Function: Compare the descriptors (id,accession_number,display_id,
|
|
|
373 primary_id, desc) for the PrimarySeq and PrimaryQual objects
|
|
|
374 within the SeqWithQuality object. If they match, make that
|
|
|
375 descriptor the descriptor for the SeqWithQuality object.
|
|
|
376 Returns : Nothing.
|
|
|
377 Args : None.
|
|
|
378
|
|
|
379 =cut
|
|
|
380
|
|
|
381 sub set_common_descriptors {
|
|
|
382 my $self = shift;
|
|
|
383 return if ($self->{seq_ref} || $self->{qual_ref});
|
|
|
384 &_common_id();
|
|
|
385 &_common_display_id();
|
|
|
386 &_common_accession_number();
|
|
|
387 &_common_primary_id();
|
|
|
388 &_common_desc();
|
|
|
389 }
|
|
|
390
|
|
|
391 =head2 alphabet()
|
|
|
392
|
|
|
393 Title : alphabet();
|
|
|
394 Usage : $molecule_type = $obj->alphabet();
|
|
|
395 Function: Get the molecule type from the PrimarySeq object.
|
|
|
396 Returns : What what PrimarySeq says the type of the sequence is.
|
|
|
397 Args : None.
|
|
|
398
|
|
|
399 =cut
|
|
|
400
|
|
|
401 sub alphabet {
|
|
|
402 my $self = shift;
|
|
|
403 return $self->{seq_ref}->alphabet();
|
|
|
404 }
|
|
|
405
|
|
|
406 =head2 display_id()
|
|
|
407
|
|
|
408 Title : display_id()
|
|
|
409 Usage : $id_string = $obj->display_id();
|
|
|
410 Function: Returns the display id, aka the common name of the Quality
|
|
|
411 object.
|
|
|
412 The semantics of this is that it is the most likely string to be
|
|
|
413 used as an identifier of the quality sequence, and likely to have
|
|
|
414 "human" readability. The id is equivalent to the ID field of the
|
|
|
415 GenBank/EMBL databanks and the id field of the Swissprot/sptrembl
|
|
|
416 database. In fasta format, the >(\S+) is presumed to be the id,
|
|
|
417 though some people overload the id to embed other information.
|
|
|
418 Bioperl does not use any embedded information in the ID field,
|
|
|
419 and people are encouraged to use other mechanisms (accession
|
|
|
420 field for example, or extending the sequence object) to solve
|
|
|
421 this. Notice that $seq->id() maps to this function, mainly for
|
|
|
422 legacy/convience issues.
|
|
|
423 This method sets the display_id for the SeqWithQuality object.
|
|
|
424 Returns : A string
|
|
|
425 Args : If a scalar is provided, it is set as the new display_id for
|
|
|
426 the SeqWithQuality object.
|
|
|
427 Status : Virtual
|
|
|
428
|
|
|
429 =cut
|
|
|
430
|
|
|
431 sub display_id {
|
|
|
432 my ($obj,$value) = @_;
|
|
|
433 if( defined $value) {
|
|
|
434 $obj->{'display_id'} = $value;
|
|
|
435 }
|
|
|
436 return $obj->{'display_id'};
|
|
|
437
|
|
|
438 }
|
|
|
439
|
|
|
440 =head2 accession_number()
|
|
|
441
|
|
|
442 Title : accession_number()
|
|
|
443 Usage : $unique_biological_key = $obj->accession_number();
|
|
|
444 Function: Returns the unique biological id for a sequence, commonly
|
|
|
445 called the accession_number. For sequences from established
|
|
|
446 databases, the implementors should try to use the correct
|
|
|
447 accession number. Notice that primary_id() provides the unique id
|
|
|
448 for the implemetation, allowing multiple objects to have the same
|
|
|
449 accession number in a particular implementation. For sequences
|
|
|
450 with no accession number, this method should return "unknown".
|
|
|
451 This method sets the accession_number for the SeqWithQuality
|
|
|
452 object.
|
|
|
453 Returns : A string (the value of accession_number)
|
|
|
454 Args : If a scalar is provided, it is set as the new accession_number
|
|
|
455 for the SeqWithQuality object.
|
|
|
456 Status : Virtual
|
|
|
457
|
|
|
458
|
|
|
459 =cut
|
|
|
460
|
|
|
461 sub accession_number {
|
|
|
462 my( $obj, $acc ) = @_;
|
|
|
463
|
|
|
464 if (defined $acc) {
|
|
|
465 $obj->{'accession_number'} = $acc;
|
|
|
466 } else {
|
|
|
467 $acc = $obj->{'accession_number'};
|
|
|
468 $acc = 'unknown' unless defined $acc;
|
|
|
469 }
|
|
|
470 return $acc;
|
|
|
471 }
|
|
|
472
|
|
|
473 =head2 primary_id()
|
|
|
474
|
|
|
475 Title : primary_id()
|
|
|
476 Usage : $unique_implementation_key = $obj->primary_id();
|
|
|
477 Function: Returns the unique id for this object in this implementation.
|
|
|
478 This allows implementations to manage their own object ids in a
|
|
|
479 way the implementaiton can control clients can expect one id to
|
|
|
480 map to one object. For sequences with no accession number, this
|
|
|
481 method should return a stringified memory location.
|
|
|
482 This method sets the primary_id for the SeqWithQuality
|
|
|
483 object.
|
|
|
484 Returns : A string. (the value of primary_id)
|
|
|
485 Args : If a scalar is provided, it is set as the new primary_id for
|
|
|
486 the SeqWithQuality object.
|
|
|
487
|
|
|
488 =cut
|
|
|
489
|
|
|
490 sub primary_id {
|
|
|
491 my ($obj,$value) = @_;
|
|
|
492 if ($value) {
|
|
|
493 $obj->{'primary_id'} = $value;
|
|
|
494 }
|
|
|
495 return $obj->{'primary_id'};
|
|
|
496
|
|
|
497 }
|
|
|
498
|
|
|
499 =head2 desc()
|
|
|
500
|
|
|
501 Title : desc()
|
|
|
502 Usage : $qual->desc($newval); _or_
|
|
|
503 $description = $qual->desc();
|
|
|
504 Function: Get/set description text for this SeqWithQuality object.
|
|
|
505 Returns : A string. (the value of desc)
|
|
|
506 Args : If a scalar is provided, it is set as the new desc for the
|
|
|
507 SeqWithQuality object.
|
|
|
508
|
|
|
509 =cut
|
|
|
510
|
|
|
511 sub desc {
|
|
|
512 # a mechanism to set the disc for the SeqWithQuality object.
|
|
|
513 # probably will be used most often by set_common_features()
|
|
|
514 my ($obj,$value) = @_;
|
|
|
515 if( defined $value) {
|
|
|
516 $obj->{'desc'} = $value;
|
|
|
517 }
|
|
|
518 return $obj->{'desc'};
|
|
|
519 }
|
|
|
520
|
|
|
521 =head2 id()
|
|
|
522
|
|
|
523 Title : id()
|
|
|
524 Usage : $id = $qual->id();
|
|
|
525 Function: Return the ID of the quality. This should normally be (and
|
|
|
526 actually is in the implementation provided here) just a synonym
|
|
|
527 for display_id().
|
|
|
528 Returns : A string. (the value of id)
|
|
|
529 Args : If a scalar is provided, it is set as the new id for the
|
|
|
530 SeqWithQuality object.
|
|
|
531
|
|
|
532 =cut
|
|
|
533
|
|
|
534 sub id {
|
|
|
535 my ($self,$value) = @_;
|
|
|
536 if (!$self) { $self->throw("no value for self in $value"); }
|
|
|
537 if( defined $value ) {
|
|
|
538 return $self->display_id($value);
|
|
|
539 }
|
|
|
540 return $self->display_id();
|
|
|
541 }
|
|
|
542
|
|
|
543 =head2 seq
|
|
|
544
|
|
|
545 Title : seq()
|
|
|
546 Usage : $string = $obj->seq(); _or_
|
|
|
547 $obj->seq("atctatcatca");
|
|
|
548 Function: Returns the sequence that is contained in the imbedded in the
|
|
|
549 PrimarySeq object within the SeqWithQuality object
|
|
|
550 Returns : A scalar (the seq() value for the imbedded PrimarySeq object.)
|
|
|
551 Args : If a scalar is provided, the SeqWithQuality object will
|
|
|
552 attempt to set that as the sequence for the imbedded PrimarySeq
|
|
|
553 object. Otherwise, the value of seq() for the PrimarySeq object
|
|
|
554 is returned.
|
|
|
555 Notes : This is probably not a good idea because you then should call
|
|
|
556 length() to make sure that the sequence and quality are of the
|
|
|
557 same length. Even then, how can you make sure that this sequence
|
|
|
558 belongs with that quality? I provided this to give you rope to
|
|
|
559 hang yourself with. Tie it to a strong device and use a good
|
|
|
560 knot.
|
|
|
561
|
|
|
562 =cut
|
|
|
563
|
|
|
564 sub seq {
|
|
|
565 my ($self,$value) = @_;
|
|
|
566 if( defined $value) {
|
|
|
567 $self->{seq_ref}->seq($value);
|
|
|
568 $self->length();
|
|
|
569 }
|
|
|
570 return $self->{seq_ref}->seq();
|
|
|
571 }
|
|
|
572
|
|
|
573 =head2 qual()
|
|
|
574
|
|
|
575 Title : qual()
|
|
|
576 Usage : @quality_values = @{$obj->qual()}; _or_
|
|
|
577 $obj->qual("10 10 20 40 50");
|
|
|
578 Function: Returns the quality as imbedded in the PrimaryQual object
|
|
|
579 within the SeqWithQuality object.
|
|
|
580 Returns : A reference to an array containing the quality values in the
|
|
|
581 PrimaryQual object.
|
|
|
582 Args : If a scalar is provided, the SeqWithQuality object will
|
|
|
583 attempt to set that as the quality for the imbedded PrimaryQual
|
|
|
584 object. Otherwise, the value of qual() for the PrimaryQual
|
|
|
585 object is returned.
|
|
|
586 Notes : This is probably not a good idea because you then should call
|
|
|
587 length() to make sure that the sequence and quality are of the
|
|
|
588 same length. Even then, how can you make sure that this sequence
|
|
|
589 belongs with that quality? I provided this to give you a strong
|
|
|
590 board with which to flagellate yourself.
|
|
|
591
|
|
|
592 =cut
|
|
|
593
|
|
|
594 sub qual {
|
|
|
595 my ($self,$value) = @_;
|
|
|
596
|
|
|
597 if( defined $value) {
|
|
|
598 $self->{qual_ref}->qual($value);
|
|
|
599 # update the lengths
|
|
|
600 $self->length();
|
|
|
601 }
|
|
|
602 return $self->{qual_ref}->qual();
|
|
|
603 }
|
|
|
604
|
|
|
605
|
|
|
606
|
|
|
607 =head2 trace_indices()
|
|
|
608
|
|
|
609 Title : trace_indices()
|
|
|
610 Usage : @trace_indice_values = @{$obj->trace_indices()}; _or_
|
|
|
611 $obj->trace_indices("10 10 20 40 50");
|
|
|
612 Function: Returns the trace_indices as imbedded in the Primaryqual object
|
|
|
613 within the SeqWithQualiity object.
|
|
|
614 Returns : A reference to an array containing the trace_indice values in the
|
|
|
615 PrimaryQual object.
|
|
|
616 Args : If a scalar is provided, the SeqWithuQuality object will
|
|
|
617 attempt to set that as the trace_indices for the imbedded PrimaryQual
|
|
|
618 object. Otherwise, the value of trace_indices() for the PrimaryQual
|
|
|
619 object is returned.
|
|
|
620 Notes : This is probably not a good idea because you then should call
|
|
|
621 length() to make sure that the sequence and trace_indices are of the
|
|
|
622 same length. Even then, how can you make sure that this sequence
|
|
|
623 belongs with that trace_indicex? I provided this to give you a strong
|
|
|
624 board with which to flagellate yourself.
|
|
|
625
|
|
|
626 =cut
|
|
|
627
|
|
|
628 sub trace_indices {
|
|
|
629 my ($self,$value) = @_;
|
|
|
630
|
|
|
631 if( defined $value) {
|
|
|
632 $self->{qual_ref}->trace_indices($value);
|
|
|
633 # update the lengths
|
|
|
634 $self->length();
|
|
|
635 }
|
|
|
636 return $self->{qual_ref}->trace_indices();
|
|
|
637 }
|
|
|
638
|
|
|
639
|
|
|
640
|
|
|
641
|
|
|
642 =head2 length()
|
|
|
643
|
|
|
644 Title : length()
|
|
|
645 Usage : $length = $seqWqual->length();
|
|
|
646 Function: Get the length of the SeqWithQuality sequence/quality.
|
|
|
647 Returns : Returns the length of the sequence and quality if they are
|
|
|
648 both the same. Returns "DIFFERENT" if they differ.
|
|
|
649 Args : None.
|
|
|
650
|
|
|
651 =cut
|
|
|
652
|
|
|
653 sub length {
|
|
|
654 my $self = shift;
|
|
|
655 if (!$self->{seq_ref}) {
|
|
|
656 unless ($self->{supress_warnings} == 1) {
|
|
|
657 $self->warn("Can't find {seq_ref} here in length().");
|
|
|
658 }
|
|
|
659 return;
|
|
|
660 }
|
|
|
661 if (!$self->{qual_ref}) {
|
|
|
662 unless ($self->{supress_warnings} == 1) {
|
|
|
663 $self->warn("Can't find {qual_ref} here in length().");
|
|
|
664 }
|
|
|
665 return;
|
|
|
666 }
|
|
|
667 my $seql = $self->{seq_ref}->length();
|
|
|
668
|
|
|
669 if ($seql != $self->{qual_ref}->length()) {
|
|
|
670 unless ($self->{supress_warnings} == 1) {
|
|
|
671 $self->warn("Sequence length (".$seql.") is different from quality length (".$self->{qual_ref}->length().") in the SeqWithQuality object. This can only lead to problems later.");
|
|
|
672 }
|
|
|
673 $self->{'length'} = "DIFFERENT";
|
|
|
674 }
|
|
|
675 else {
|
|
|
676 $self->{'length'} = $seql;
|
|
|
677 }
|
|
|
678 return $self->{'length'};
|
|
|
679 }
|
|
|
680
|
|
|
681
|
|
|
682 =head2 qual_obj
|
|
|
683
|
|
|
684 Title : qual_obj($different_obj)
|
|
|
685 Usage : $qualobj = $seqWqual->qual_obj(); _or_
|
|
|
686 $qualobj = $seqWqual->qual_obj($ref_to_primaryqual_obj);
|
|
|
687 Function: Get the PrimaryQual object that is imbedded in the
|
|
|
688 SeqWithQuality object or if a reference to a PrimaryQual object
|
|
|
689 is provided, set this as the PrimaryQual object imbedded in the
|
|
|
690 SeqWithQuality object.
|
|
|
691 Returns : A reference to a Bio::Seq::SeqWithQuality object.
|
|
|
692
|
|
|
693 =cut
|
|
|
694
|
|
|
695 sub qual_obj {
|
|
|
696 my ($self,$value) = @_;
|
|
|
697 if (defined($value)) {
|
|
|
698 if (ref($value) eq "Bio::Seq::PrimaryQual") {
|
|
|
699 $self->{qual_ref} = $value;
|
|
|
700
|
|
|
701 $self->debug("You successfully changed the PrimaryQual object within a SeqWithQuality object. ID's for the SeqWithQuality object may now not be what you expect. Use something like set_common_descriptors() to fix them if you care,");
|
|
|
702 }
|
|
|
703 else {
|
|
|
704 $self->debug("You tried to change the PrimaryQual object within a SeqWithQuality object but you passed a reference to an object that was not a Bio::Seq::PrimaryQual object. Thus your change failed. Sorry.\n");
|
|
|
705 }
|
|
|
706 }
|
|
|
707 return $self->{qual_ref};
|
|
|
708 }
|
|
|
709
|
|
|
710
|
|
|
711 =head2 seq_obj
|
|
|
712
|
|
|
713 Title : seq_obj()
|
|
|
714 Usage : $seqobj = $seqWqual->qual_obj(); _or_
|
|
|
715 $seqobj = $seqWqual->seq_obj($ref_to_primary_seq_obj);
|
|
|
716 Function: Get the PrimarySeq object that is imbedded in the
|
|
|
717 SeqWithQuality object or if a reference to a PrimarySeq object is
|
|
|
718 provided, set this as the PrimarySeq object imbedded in the
|
|
|
719 SeqWithQuality object.
|
|
|
720 Returns : A reference to a Bio::PrimarySeq object.
|
|
|
721
|
|
|
722 =cut
|
|
|
723
|
|
|
724 sub seq_obj {
|
|
|
725 my ($self,$value) = @_;
|
|
|
726 if( defined $value) {
|
|
|
727 if (ref($value) eq "Bio::PrimarySeq") {
|
|
|
728 $self->debug("You successfully changed the PrimarySeq object within a SeqWithQuality object. ID's for the SeqWithQuality object may now not be what you expect. Use something like set_common_descriptors() to fix them if you care,");
|
|
|
729 } else {
|
|
|
730 $self->debug("You tried to change the PrimarySeq object within a SeqWithQuality object but you passed a reference to an object that was not a Bio::PrimarySeq object. Thus your change failed. Sorry.\n");
|
|
|
731 }
|
|
|
732 }
|
|
|
733 return $self->{seq_ref};
|
|
|
734 }
|
|
|
735
|
|
|
736 =head2 _set_descriptors
|
|
|
737
|
|
|
738 Title : _set_descriptors()
|
|
|
739 Usage : $seqWqual->_qual_obj($qual,$seq,$id,$acc,$pid,$desc,$given_id,
|
|
|
740 $alphabet);
|
|
|
741 Function: Set the descriptors for the SeqWithQuality object. Try to
|
|
|
742 match the descriptors in the PrimarySeq object and in the
|
|
|
743 PrimaryQual object if descriptors were not provided with
|
|
|
744 construction.
|
|
|
745 Returns : Nothing.
|
|
|
746 Args : $qual,$seq,$id,$acc,$pid,$desc,$given_id,$alphabet as found
|
|
|
747 in the new() method.
|
|
|
748 Notes : Really only intended to be called by the new() method. If
|
|
|
749 you want to invoke a similar function try
|
|
|
750 set_common_descriptors().
|
|
|
751
|
|
|
752 =cut
|
|
|
753
|
|
|
754
|
|
|
755 sub _set_descriptors {
|
|
|
756 my ($self,$qual,$seq,$id,$acc,$pid,$desc,$given_id,$alphabet) = @_;
|
|
|
757 my ($c_id,$c_acc,$c_pid,$c_desc);
|
|
|
758 if (!$self->display_id()) {
|
|
|
759 if ($c_id = $self->_common_id() ) { $self->display_id($c_id); }
|
|
|
760 else {
|
|
|
761 if ($self->{seq_ref}) {
|
|
|
762 # print("Using seq_ref to set id to ".$self->{seq_ref}->display_id()."\n");
|
|
|
763 # ::dumpValue($self->{seq_ref});
|
|
|
764 $self->display_id($self->{seq_ref}->id());
|
|
|
765 }
|
|
|
766 elsif ($self->{qual_ref}) {
|
|
|
767 $self->display_id($self->{qual_ref}->id());
|
|
|
768 }
|
|
|
769 }
|
|
|
770 }
|
|
|
771 if ($acc) { $self->accession_number($acc); }
|
|
|
772 elsif ($c_acc = $self->_common_accession_number() ) { $self->accession_number($c_acc); }
|
|
|
773 if ($pid) { $self->primary_id($pid); }
|
|
|
774 elsif ($c_pid = $self->_common_primary_id() ) { $self->primary_id($c_pid); }
|
|
|
775 if ($desc) { $self->desc($desc); }
|
|
|
776 elsif ($c_desc = $self->_common_desc() ) { $self->desc($c_desc); }
|
|
|
777 }
|
|
|
778
|
|
|
779 =head2 subseq($start,$end)
|
|
|
780
|
|
|
781 Title : subseq($start,$end)
|
|
|
782 Usage : $subsequence = $obj->subseq($start,$end);
|
|
|
783 Function: Returns the subseq from start to end, where the first base
|
|
|
784 is 1 and the number is inclusive, ie 1-2 are the first two
|
|
|
785 bases of the sequence.
|
|
|
786 Returns : A string.
|
|
|
787 Args : Two positions.
|
|
|
788
|
|
|
789 =cut
|
|
|
790
|
|
|
791 sub subseq {
|
|
|
792 my ($self,@args) = @_;
|
|
|
793 # does a single value work?
|
|
|
794 return $self->{seq_ref}->subseq(@args);
|
|
|
795 }
|
|
|
796
|
|
|
797 =head2 baseat($position)
|
|
|
798
|
|
|
799 Title : baseat($position)
|
|
|
800 Usage : $base_at_position_6 = $obj->baseat("6");
|
|
|
801 Function: Returns a single base at the given position, where the first
|
|
|
802 base is 1 and the number is inclusive, ie 1-2 are the first two
|
|
|
803 bases of the sequence.
|
|
|
804 Returns : A scalar.
|
|
|
805 Args : A position.
|
|
|
806
|
|
|
807 =cut
|
|
|
808
|
|
|
809 sub baseat {
|
|
|
810 my ($self,$val) = @_;
|
|
|
811 return $self->{seq_ref}->subseq($val,$val);
|
|
|
812 }
|
|
|
813
|
|
|
814 =head2 subqual($start,$end)
|
|
|
815
|
|
|
816 Title : subqual($start,$end)
|
|
|
817 Usage : @qualities = @{$obj->subqual(10,20);
|
|
|
818 Function: returns the quality values from $start to $end, where the
|
|
|
819 first value is 1 and the number is inclusive, ie 1-2 are the
|
|
|
820 first two bases of the sequence. Start cannot be larger than
|
|
|
821 end but can be equal.
|
|
|
822 Returns : A reference to an array.
|
|
|
823 Args : a start position and an end position
|
|
|
824
|
|
|
825 =cut
|
|
|
826
|
|
|
827 sub subqual {
|
|
|
828 my ($self,@args) = @_;
|
|
|
829 return $self->{qual_ref}->subqual(@args);
|
|
|
830 }
|
|
|
831
|
|
|
832 =head2 qualat($position)
|
|
|
833
|
|
|
834 Title : qualat($position)
|
|
|
835 Usage : $quality = $obj->qualat(10);
|
|
|
836 Function: Return the quality value at the given location, where the
|
|
|
837 first value is 1 and the number is inclusive, ie 1-2 are the
|
|
|
838 first two bases of the sequence. Start cannot be larger than
|
|
|
839 end but can be equal.
|
|
|
840 Returns : A scalar.
|
|
|
841 Args : A position.
|
|
|
842
|
|
|
843 =cut
|
|
|
844
|
|
|
845 sub qualat {
|
|
|
846 my ($self,$val) = @_;
|
|
|
847 return $self->{qual_ref}->qualat($val);
|
|
|
848 }
|
|
|
849
|
|
|
850 =head2 sub_trace_index($start,$end)
|
|
|
851
|
|
|
852 Title : sub_trace_index($start,$end)
|
|
|
853 Usage : @trace_indices = @{$obj->sub_trace_index(10,20);
|
|
|
854 Function: returns the trace index values from $start to $end, where the
|
|
|
855 first value is 1 and the number is inclusive, ie 1-2 are the
|
|
|
856 first two bases of the sequence. Start cannot be larger than
|
|
|
857 end but can be e_trace_index.
|
|
|
858 Returns : A reference to an array.
|
|
|
859 Args : a start position and an end position
|
|
|
860
|
|
|
861 =cut
|
|
|
862
|
|
|
863 sub sub_trace_index {
|
|
|
864 my ($self,@args) = @_;
|
|
|
865 return $self->{qual_ref}->sub_trace_index(@args);
|
|
|
866 }
|
|
|
867
|
|
|
868 =head2 trace_index_at($position)
|
|
|
869
|
|
|
870 Title : trace_index_at($position)
|
|
|
871 Usage : $trace_index = $obj->trace_index_at(10);
|
|
|
872 Function: Return the trace_index value at the given location, where the
|
|
|
873 first value is 1 and the number is inclusive, ie 1-2 are the
|
|
|
874 first two bases of the sequence. Start cannot be larger than
|
|
|
875 end but can be etrace_index_.
|
|
|
876 Returns : A scalar.
|
|
|
877 Args : A position.
|
|
|
878
|
|
|
879 =cut
|
|
|
880
|
|
|
881 sub trace_index_at {
|
|
|
882 my ($self,$val) = @_;
|
|
|
883 return $self->{qual_ref}->trace_index_at($val);
|
|
|
884 }
|
|
|
885
|
|
|
886 =head2 to_string()
|
|
|
887
|
|
|
888 Title : to_string()
|
|
|
889 Usage : $quality = $obj->to_string();
|
|
|
890 Function: Return a textual representation of what the object contains.
|
|
|
891 For this module, this function will return:
|
|
|
892 qual
|
|
|
893 seq
|
|
|
894 display_id
|
|
|
895 accession_number
|
|
|
896 primary_id
|
|
|
897 desc
|
|
|
898 id
|
|
|
899 length_sequence
|
|
|
900 length_quality
|
|
|
901 Returns : A scalar.
|
|
|
902 Args : None.
|
|
|
903
|
|
|
904 =cut
|
|
|
905
|
|
|
906 sub to_string {
|
|
|
907 my ($self,$out,$result) = shift;
|
|
|
908 $out = "qual: ".join(',',@{$self->qual()})."\n";
|
|
|
909 foreach (qw(seq display_id accession_number primary_id desc id)) {
|
|
|
910 $result = $self->$_();
|
|
|
911 if (!$result) { $result = "<unset>"; }
|
|
|
912 $out .= "$_: $result\n";
|
|
|
913 }
|
|
|
914 return $out;
|
|
|
915 }
|
|
|
916 1;
|
