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; |
