lib/Bio/Seq/SeqWithQuality.pm

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