~ubuntu-branches/ubuntu/trusty/bioperl/trusty

« back to all changes in this revision

Viewing changes to Bio/Seq/PrimaryQual.pm

  • Committer: Package Import Robot
  • Author(s): Charles Plessy
  • Date: 2013-09-22 13:39:48 UTC
  • mfrom: (3.1.11 sid)
  • Revision ID: package-import@ubuntu.com-20130922133948-c6z62zegjyp7ztou
Tags: 1.6.922-1
* New upstream release.
* Replaces and Breaks grinder (<< 0.5.3-3~) because of overlaping contents.
  Closes: #722910
* Stop Replacing and Breaking bioperl ( << 1.6.9 ): not needed anymore. 

Show diffs side-by-side

added added

removed removed

Lines of Context:
22
22
 # you can use either a space-delimited string for quality
23
23
 
24
24
 my $string_quals = "10 20 30 40 50 40 30 20 10";
25
 
 my $qualobj = Bio::Seq::PrimaryQual->new
26
 
              ( '-qual' => $string_quals,
27
 
                '-id'  => 'QualityFragment-12',
28
 
                '-accession_number' => 'X78121',
29
 
                );
 
25
 my $qualobj = Bio::Seq::PrimaryQual->new(
 
26
     -qual             => $string_quals,
 
27
     -id               => 'QualityFragment-12',
 
28
     -accession_number => 'X78121',
 
29
 );
30
30
 
31
31
 # _or_ you can use an array of quality values
32
32
 
33
33
 my @q2 = split/ /,$string_quals;
34
 
 $qualobj = Bio::Seq::PrimaryQual->new( '-qual' => \@q2,
35
 
       '-primary_id'     =>      'chads primary_id',
36
 
       '-desc'           =>      'chads desc',
37
 
       '-accession_number' => 'chads accession_number',
38
 
      '-id'             =>      'chads id'
39
 
      );
 
34
 $qualobj = Bio::Seq::PrimaryQual->new(
 
35
     -qual              => \@q2,
 
36
     -primary_id        => 'chads primary_id',
 
37
     -desc              => 'chads desc',
 
38
     -accession_number  => 'chads accession_number',
 
39
     -id                => 'chads id'
 
40
 );
40
41
 
41
42
 # to get the quality values out:
42
43
 
96
97
=cut
97
98
 
98
99
 
99
 
# Let the code begin...
100
 
 
101
100
package Bio::Seq::PrimaryQual;
102
 
use vars qw(%valid_type);
 
101
 
103
102
use strict;
104
103
 
105
 
 
106
104
use base qw(Bio::Root::Root Bio::Seq::QualI);
107
105
 
 
106
our $MATCHPATTERN = '0-9eE\.\s+-';
 
107
 
108
108
 
109
109
=head2 new()
110
110
 
111
111
 Title   : new()
112
 
 Usage   : $qual = Bio::Seq::PrimaryQual->new
113
 
        ( -qual => '10 20 30 40 50 50 20 10',
114
 
          -id  => 'human_id',
115
 
          -accession_number => 'AL000012',
116
 
        );
 
112
 Usage   : $qual = Bio::Seq::PrimaryQual->new(
 
113
               -qual             => '10 20 30 40 50 50 20 10',
 
114
               -id               => 'human_id',
 
115
               -accession_number => 'AL000012',
 
116
           );
117
117
 
118
118
 Function: Returns a new Bio::Seq::PrimaryQual object from basic 
119
 
        constructors, being a string _or_ a reference to an array for the
120
 
        sequence and strings for id and accession_number. Note that you
121
 
        can provide an empty quality string.
 
119
           constructors, being a string _or_ a reference to an array for the
 
120
           sequence and strings for id and accession_number. Note that you
 
121
           can provide an empty quality string.
122
122
 Returns : a new Bio::Seq::PrimaryQual object
123
123
 
124
124
=cut
125
125
 
126
 
 
127
 
 
128
 
 
129
126
sub new {
130
127
    my ($class, @args) = @_;
131
128
    my $self = $class->SUPER::new(@args);
159
156
    return $self;
160
157
}
161
158
 
 
159
 
162
160
=head2 qual()
163
161
 
164
162
 Title   : qual()
165
163
 Usage   : @quality_values  = @{$obj->qual()};
166
 
 Function: Returns the quality as a reference to an array containing the
167
 
           quality values. The individual elements of the quality array are
168
 
           not validated and can be any numeric value.
 
164
 Function: Get or set the quality as a reference to an array containing the
 
165
           quality values. An error is generated if the quality scores are
 
166
           invalid, see validate_qual().
169
167
 Returns : A reference to an array.
170
168
 
171
169
=cut
174
172
    my ($self,$value) = @_;
175
173
 
176
174
    if( ! defined $value || length($value) == 0 ) { 
177
 
        $self->{'qual'} ||= [];
 
175
        $self->{'qual'} ||= [];
178
176
    } elsif( ref($value) =~ /ARRAY/i ) {
179
 
        # if the user passed in a reference to an array
180
 
        $self->{'qual'} = $value;
181
 
    } elsif(! $self->validate_qual($value)){
182
 
        $self->throw("Attempting to set the quality to [$value] which does not look healthy");      
 
177
        # if the user passed in a reference to an array
 
178
        $self->{'qual'} = $value;
183
179
    } else {
184
 
        $value =~ s/^\s+//;
185
 
        $self->{'qual'} = [split(/\s+/,$value)];
 
180
        $self->validate_qual($value, 1);
 
181
        $value =~ s/^\s+//;
 
182
        $self->{'qual'} = [split(/\s+/,$value)];
186
183
    }
187
184
    
188
185
    return $self->{'qual'};
200
197
=cut
201
198
 
202
199
sub seq {
203
 
    join ' ', @{ shift->qual };
 
200
    return join ' ', @{ shift->qual };
204
201
}
205
202
 
206
203
 
207
204
=head2 validate_qual($qualstring)
208
205
 
209
 
 Title   : validate_qual($qualstring)
210
 
 Usage   : print("Valid.") if { &validate_qual($self,$qualities); }
211
 
 Function: Make sure that the quality, if it has length > 0, contains at
212
 
        least one digit. Note that quality strings are parsed into arrays
213
 
        using split/\d+/,$quality_string, so make sure that your quality
214
 
        scalar looks like this if you want it to be parsed properly.
215
 
 Returns : 1 for a valid sequence (WHY? Shouldn\'t it return 0? <boggle>)
216
 
 Args    : a scalar (any scalar, why PrimarySeq author?) and a scalar
217
 
        containing the string to validate.
 
206
 Title   : validate_qual($qualstring)
 
207
 Usage   : print("Valid.") if { &validate_qual($self, $quality_string); }
 
208
 Function: Test that the given quality string is valid. It is expected to
 
209
           contain space-delimited numbers that can be parsed using split /\d+/.
 
210
           However, this validation takes shortcuts and only tests that the
 
211
           string contains characters valid in numbers: 0-9 . eE +-
 
212
           Note that empty quality strings are valid too.
 
213
 Returns : 1 for a valid sequence, 0 otherwise
 
214
 Args    : - Scalar containing the quality string to validate.
 
215
           - Boolean to optionally throw an error if validation failed
218
216
 
219
217
=cut
220
218
 
221
219
sub validate_qual {
222
 
    # how do I validate quality values?
223
 
    # \d+\s+\d+..., I suppose
224
 
    my ($self,$qualstr) = @_;
225
 
    # why the CORE?? -- (Because Bio::PrimarySeqI namespace has a 
226
 
    #                    length method, you have to qualify 
227
 
    #                    which length to use)
228
 
    return 0 if (!defined $qualstr || CORE::length($qualstr) <= 0);   
229
 
    return 1 if( $qualstr =~ /\d/);
230
 
    
231
 
    return 0;
 
220
    my ($self, $qualstr, $throw) = @_;
 
221
    if ( (defined $qualstr                ) &&
 
222
         ($qualstr !~ /^[$MATCHPATTERN]*$/) ) {
 
223
        if ($throw) {
 
224
            $self->throw("Failed validation of quality score from  '".
 
225
               (defined($self->id)||'[unidentified sequence]')."'. No numeric ".
 
226
               "value found.\n");
 
227
        }
 
228
        return 0;
 
229
    }
 
230
    return 1;
232
231
}
233
232
 
 
233
 
234
234
=head2 subqual($start,$end)
235
235
 
236
236
 Title   : subqual($start,$end)
237
237
 Usage   : @subset_of_quality_values = @{$obj->subqual(10,40)};
238
238
 Function: returns the quality values from $start to $end, where the
239
 
        first value is 1 and the number is inclusive, ie 1-2 are the
240
 
        first two bases of the sequence. Start cannot be larger than
241
 
        end but can be equal.
 
239
           first value is 1 and the number is inclusive, ie 1-2 are the
 
240
           first two bases of the sequence. Start cannot be larger than
 
241
           end but can be equal.
242
242
 Returns : A reference to an array.
243
243
 Args    : a start position and an end position
244
244
 
245
245
=cut
246
246
 
247
 
 
248
247
sub subqual {
249
 
   my ($self,$start,$end) = @_;
250
 
 
251
 
   if( $start > $end ){
252
 
       $self->throw("in subqual, start [$start] has to be greater than end [$end]");
253
 
   }
254
 
 
255
 
   if( $start <= 0 || $end > $self->length ) {
256
 
       $self->throw("You have to have start positive and length less than the total length of sequence [$start:$end] Total ".$self->length."");
257
 
   }
258
 
 
259
 
   # remove one from start, and then length is end-start
260
 
 
261
 
   $start--;
262
 
        $end--;
263
 
        my @sub_qual_array = @{$self->{qual}}[$start..$end];
264
 
 
265
 
        #   return substr $self->seq(), $start, ($end-$start);
266
 
        return \@sub_qual_array;
 
248
    my ($self,$start,$end) = @_;
 
249
 
 
250
    if( $start > $end ){
 
251
        $self->throw("in subqual, start [$start] has to be greater than end [$end]");
 
252
    }
 
253
 
 
254
    if( $start <= 0 || $end > $self->length ) {
 
255
        $self->throw("You have to have start positive and length less than the total length of sequence [$start:$end] Total ".$self->length."");
 
256
    }
 
257
 
 
258
    # remove one from start, and then length is end-start
 
259
 
 
260
    $start--;
 
261
    $end--;
 
262
    my @sub_qual_array = @{$self->{qual}}[$start..$end];
 
263
 
 
264
    #   return substr $self->seq(), $start, ($end-$start);
 
265
    return \@sub_qual_array;
267
266
 
268
267
}
269
268
 
 
269
 
270
270
=head2 display_id()
271
271
 
272
272
 Title   : display_id()
273
273
 Usage   : $id_string = $obj->display_id();
274
274
 Function: returns the display id, aka the common name of the Quality
275
 
        object.
276
 
        The semantics of this is that it is the most likely string to be
277
 
        used as an identifier of the quality sequence, and likely to have
278
 
        "human" readability.  The id is equivalent to the ID field of the
279
 
        GenBank/EMBL databanks and the id field of the Swissprot/sptrembl
280
 
        database. In fasta format, the >(\S+) is presumed to be the id,
281
 
        though some people overload the id to embed other information.
282
 
        Bioperl does not use any embedded information in the ID field,
283
 
        and people are encouraged to use other mechanisms (accession
284
 
        field for example, or extending the sequence object) to solve
285
 
        this. Notice that $seq->id() maps to this function, mainly for
286
 
        legacy/convience issues
 
275
           object.
 
276
           The semantics of this is that it is the most likely string to be
 
277
           used as an identifier of the quality sequence, and likely to have
 
278
           "human" readability.  The id is equivalent to the ID field of the
 
279
           GenBank/EMBL databanks and the id field of the Swissprot/sptrembl
 
280
           database. In fasta format, the >(\S+) is presumed to be the id,
 
281
           though some people overload the id to embed other information.
 
282
           Bioperl does not use any embedded information in the ID field,
 
283
           and people are encouraged to use other mechanisms (accession
 
284
           field for example, or extending the sequence object) to solve
 
285
           this. Notice that $seq->id() maps to this function, mainly for
 
286
           legacy/convience issues
287
287
 Returns : A string
288
288
 Args    : None
289
289
 
295
295
      $obj->{'display_id'} = $value;
296
296
    }
297
297
    return $obj->{'display_id'};
298
 
 
299
298
}
300
299
 
 
300
 
301
301
=head2 header()
302
302
 
303
303
 Title   : header()
318
318
 
319
319
}
320
320
 
 
321
 
321
322
=head2 accession_number()
322
323
 
323
324
 Title   : accession_number()
346
347
    return $acc;
347
348
}
348
349
 
 
350
 
349
351
=head2 primary_id()
350
352
 
351
353
 Title   : primary_id()
361
363
=cut
362
364
 
363
365
sub primary_id {
364
 
   my ($obj,$value) = @_;
365
 
   if( defined $value) {
366
 
      $obj->{'primary_id'} = $value;
 
366
    my ($obj,$value) = @_;
 
367
    if( defined $value) {
 
368
        $obj->{'primary_id'} = $value;
367
369
    }
368
 
   return $obj->{'primary_id'};
369
 
 
 
370
    return $obj->{'primary_id'};
370
371
}
371
372
 
 
373
 
372
374
=head2 desc()
373
375
 
374
376
 Title   : desc()
382
384
=cut
383
385
 
384
386
sub desc {
385
 
   my ($obj,$value) = @_;
386
 
   if( defined $value) {
387
 
      $obj->{'desc'} = $value;
 
387
    my ($obj,$value) = @_;
 
388
    if( defined $value) {
 
389
        $obj->{'desc'} = $value;
388
390
    }
389
391
    return $obj->{'desc'};
390
 
 
391
392
}
392
393
 
 
394
 
393
395
=head2 id()
394
396
 
395
397
 Title   : id()
396
398
 Usage   : $id = $qual->id();
397
399
 Function: Return the ID of the quality. This should normally be (and
398
 
        actually is in the implementation provided here) just a synonym
399
 
        for display_id().
 
400
           actually is in the implementation provided here) just a synonym
 
401
           for display_id().
400
402
 Returns : A string.
401
403
 Args    : None.
402
404
 
403
405
=cut
404
406
 
405
407
sub id {
406
 
   my ($self,$value) = @_;
407
 
   if( defined $value ) {
 
408
    my ($self,$value) = @_;
 
409
    if( defined $value ) {
408
410
        return $self->display_id($value);
409
 
   }
410
 
   return $self->display_id();
 
411
    }
 
412
    return $self->display_id();
411
413
}
412
414
 
 
415
 
413
416
=head2 length()
414
417
 
415
 
 Title   : length()
416
 
 Usage   : $length = $qual->length();
 
418
 Title   : length()
 
419
 Usage   : $length = $qual->length();
417
420
 Function: Return the length of the array holding the quality values.
418
 
        Under most circumstances, this should match the number of quality
419
 
        values but no validation is done when the PrimaryQual object is
420
 
        constructed and non-digits could be put into this array. Is this
421
 
        a bug? Just enough rope...
 
421
           Under most circumstances, this should match the number of quality
 
422
           values but no validation is done when the PrimaryQual object is
 
423
           constructed and non-digits could be put into this array. Is this
 
424
           a bug? Just enough rope...
422
425
 Returns : A scalar (the number of elements in the quality array).
423
426
 Args    : None.
424
427
 
427
430
sub length {
428
431
    my $self = shift;
429
432
    if (ref($self->{qual}) ne "ARRAY") {
430
 
        $self->warn("{qual} is not an array here. Why? It appears to be ".ref($self->{qual})."(".$self->{qual}."). Good thing this can _never_ happen.");
 
433
        $self->warn("{qual} is not an array here. Why? It appears to be ".ref($self->{qual})."(".$self->{qual}."). Good thing this can _never_ happen.");
431
434
    }
432
435
    return scalar(@{$self->{qual}});
433
436
}
434
437
 
435
 
=head2 qualat($position)
436
 
 
437
 
 Title   : qualat($position)
 
438
 
 
439
=head2 qualat()
 
440
 
 
441
 Title   : qualat
438
442
 Usage   : $quality = $obj->qualat(10);
439
443
 Function: Return the quality value at the given location, where the
440
 
        first value is 1 and the number is inclusive, ie 1-2 are the first
441
 
        two bases of the sequence. Start cannot be larger than end but can
442
 
        be equal.
 
444
           first value is 1 and the number is inclusive, ie 1-2 are the first
 
445
           two bases of the sequence. Start cannot be larger than end but can
 
446
           be equal.
443
447
 Returns : A scalar.
444
448
 Args    : A position.
445
449
 
449
453
    my ($self,$val) = @_;
450
454
    my @qualat = @{$self->subqual($val,$val)};
451
455
    if (scalar(@qualat) == 1) {
452
 
        return $qualat[0];
453
 
    }
454
 
    else {
455
 
        $self->throw("AAAH! qualat provided more then one quality.");
 
456
        return $qualat[0];
 
457
    } else {
 
458
        $self->throw("qualat() provided more than one quality.");
456
459
    }
457
460
458
461
 
461
464
 Title   : to_string()
462
465
 Usage   : $quality = $obj->to_string();
463
466
 Function: Return a textual representation of what the object contains.
464
 
        For this module, this function will return:
465
 
                qual
466
 
                display_id
467
 
                accession_number
468
 
                primary_id
469
 
                desc
470
 
                id
471
 
                length
 
467
           For this module, this function will return:
 
468
                qual
 
469
                display_id
 
470
                accession_number
 
471
                primary_id
 
472
                desc
 
473
                id
 
474
                length
472
475
 Returns : A scalar.
473
476
 Args    : None.
474
477
 
478
481
    my ($self,$out,$result) = shift;
479
482
    $out = "qual: ".join(',',@{$self->qual()});
480
483
    foreach (qw(display_id accession_number primary_id desc id length)) {
481
 
        $result = $self->$_();
482
 
        if (!$result) { $result = "<unset>"; }
483
 
        $out .= "$_: $result\n";
 
484
        $result = $self->$_();
 
485
        if (!$result) { $result = "<unset>"; }
 
486
        $out .= "$_: $result\n";
484
487
    }
485
488
    return $out;
486
489
}
489
492
sub to_string_automatic {
490
493
    my ($self,$sub_result,$out) = shift;
491
494
    foreach (sort keys %$self) {
492
 
        print("Working on $_\n");
493
 
        eval { $self->$_(); };
494
 
        if ($@) { $sub_result = ref($_); }
495
 
        elsif (!($sub_result = $self->$_())) {
496
 
            $sub_result = "<unset>";
497
 
        }
498
 
        if (ref($sub_result) eq "ARRAY") {
499
 
            print("This thing ($_) is an array!\n");
500
 
            $sub_result = join(',',@$sub_result);       
501
 
        }
502
 
        $out .= "$_: ".$sub_result."\n";
 
495
        print("Working on $_\n");
 
496
        eval { $self->$_(); };
 
497
        if ($@) { $sub_result = ref($_); }
 
498
        elsif (!($sub_result = $self->$_())) {
 
499
            $sub_result = "<unset>";
 
500
        }
 
501
        if (ref($sub_result) eq "ARRAY") {
 
502
            print("This thing ($_) is an array!\n");
 
503
            $sub_result = join(',',@$sub_result);
 
504
        }
 
505
        $out .= "$_: ".$sub_result."\n";
503
506
    }
504
507
    return $out;
505
508