Mercurial > repos > mahtabm > ensembl
comparison variant_effect_predictor/Bio/Seq/SeqWithQuality.pm @ 0:1f6dce3d34e0
Uploaded
author | mahtabm |
---|---|
date | Thu, 11 Apr 2013 02:01:53 -0400 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
-1:000000000000 | 0:1f6dce3d34e0 |
---|---|
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; |