t/Seq/Quality.t

   1 # -*-Perl-*- Test Harness script for Bioperl
   2 # $Id$
   3
   4 use strict;
   5
   6 BEGIN {
   7     use lib '.';
   8     use Bio::Root::Test;
   9
  10     test_begin(-tests => 85);
  11
  12     use_ok('Bio::Seq::Quality');
  13 }
  14
  15 use Bio::SeqIO;
  16
  17 my $DEBUG = test_debug();
  18
  19 # create some random sequence object with no id
  20 my $seqobj_broken = Bio::Seq::Quality->
  21     new( -seq => "ATCGATCGA",
  22          );
  23
  24 my $seqobj;
  25 lives_ok {
  26     $seqobj = Bio::Seq::Quality->
  27         new( -seq => "ATCGATCGA",
  28              -id  => 'QualityFragment-12',
  29              -accession_number => 'X78121',
  30              );
  31 };
  32
  33
  34 # create some random quality object with the same number of qualities
  35 # and the same identifiers
  36 my $string_quals = "10 20 30 40 50 40 30 20 10";
  37 my $qualobj;
  38 lives_ok {
  39     $qualobj = Bio::Seq::Quality->
  40         new( -qual => $string_quals,
  41              -id  => 'QualityFragment-12',
  42              -accession_number => 'X78121',
  43              );
  44 };
  45
  46 # check to see what happens when you construct the Quality object
  47 ok my $swq1 = Bio::Seq::Quality->
  48     new( -seq => "ATCGATCGA",
  49          -id  => 'QualityFragment-12',
  50          -accession_number => 'X78121',
  51          -qual  =>      $string_quals);
  52
  53
  54 print("Testing various weird constructors...\n") if $DEBUG;
  55 print("\ta) No ids, Sequence object, no quality...\n") if $DEBUG;
  56
  57 # w for weird
  58 my $wswq1;
  59 lives_ok {
  60     $wswq1 = Bio::Seq::Quality->
  61         new( -seq  => "ATCGATCGA",
  62              -qual => "");
  63 };
  64 print $@ if $DEBUG;
  65
  66
  67 print("\tb) No ids, no sequence, quality object...\n") if $DEBUG;
  68         # note that you must provide a alphabet for this one.
  69 $wswq1 = Bio::Seq::Quality->
  70     new( -seq => "",
  71          -qual => $string_quals,
  72          -alphabet => 'dna'
  73          );
  74
  75 print("\tc) Absolutely nothing. (HAHAHAHA)...\n") if $DEBUG;
  76 lives_ok {
  77     $wswq1 = Bio::Seq::Quality->new( -seq => "",
  78                                      -qual => "",
  79                                      -alphabet => 'dna'
  80                                      );
  81 };
  82
  83
  84 print("\td) Absolutely nothing but an ID\n") if $DEBUG;
  85 lives_ok {
  86     $wswq1 = Bio::Seq::Quality->
  87         new( -seq => "",
  88              -qual => "",
  89              -alphabet => 'dna',
  90              -id => 'an object with no sequence and no quality but with an id'
  91              );
  92 };
  93
  94 print("\td) No sequence, no quality, no ID...\n") if $DEBUG;
  95 warnings_like {
  96     $wswq1 = Bio::Seq::Quality->
  97         new( -seq  =>   "",
  98              -qual =>   "",
  99              -verbose => 0);
 100 } qr/not guess alphabet/i;
 101
 102 print("Testing various methods and behaviors...\n") if $DEBUG;
 103
 104 print("1. Testing the seq() method...\n") if $DEBUG;
 105         print("\t1a) get\n") if $DEBUG;
 106         my $original_seq = $swq1->seq();
 107         is ($original_seq, "ATCGATCGA");
 108         print("\t1b) set\n") if $DEBUG;
 109         ok ($swq1->seq("AAAAAAAAAAAA"));
 110         print("\t1c) get (again, to make sure the set was done.)\n") if $DEBUG;
 111         is($swq1->seq(), "AAAAAAAAAAAA");
 112         print("\tSetting the sequence back to the original value...\n") if $DEBUG;
 113         $swq1->seq($original_seq);
 114
 115
 116 print("2. Testing the qual() method...\n") if $DEBUG;
 117         print("\t2a) get\n") if $DEBUG;
 118         my @qual = @{$swq1->qual()};
 119         my $str_qual = join(' ',@qual);
 120         is $str_qual, "10 20 30 40 50 40 30 20 10";
 121         print("\t2b) set\n") if $DEBUG;
 122         ok $swq1->qual("10 10 10 10 10");
 123         print("\t2c) get (again, to make sure the set was done.)\n") if $DEBUG;
 124         my @qual2 = @{$swq1->qual()};
 125         my $str_qual2 = join(' ',@qual2);
 126         is($str_qual2, "10 10 10 10 10 0 0 0 0"); ###!
 127         print("\tSetting the quality back to the original value...\n") if $DEBUG;
 128         $swq1->qual($str_qual);
 129
 130 print("3. Testing the length() method...\n") if $DEBUG;
 131         print("\t3a) When lengths are equal...\n") if $DEBUG;
 132         is($swq1->length(), 9);
 133         print("\t3b) When lengths are different\n") if $DEBUG;
 134         $swq1->qual("10 10 10 10 10");
 135         isnt ($swq1->length(), "DIFFERENT");
 136
 137
 138 print("6. Testing the subqual() method...\n") if $DEBUG;
 139      my $t_subqual = "10 20 30 40 50 60 70 80 90";
 140      $swq1->qual($t_subqual);
 141      print("\t6d) Testing the subqual at the start (border condition)\n") if $DEBUG;
 142      ok ('10 20 30' eq join(' ',@{$swq1->subqual(1,3)}));
 143      print("\t6d) Testing the subqual at the end (border condition)\n") if $DEBUG;
 144      ok ('70 80 90' eq join(' ',@{$swq1->subqual(7,9)}));
 145      print("\t6d) Testing the subqual in the middle\n") if $DEBUG;
 146      ok ('40 50 60' eq join(' ',@{$swq1->subqual(4,6)}));
 147
 148 print("7. Testing cases where quality is zero...\n") if $DEBUG;
 149 $swq1 = Bio::Seq::Quality->new(-seq =>  'G',
 150                                -qual => '0',
 151                                );
 152 my $swq2 = Bio::Seq::Quality->new(-seq =>  'G',
 153                                   -qual => '65',
 154                                   );
 155 is $swq1->length, $swq2->length;
 156
 157 $swq1 = Bio::Seq::Quality->new(-seq =>  'GC',
 158                                -qual => '0 0',
 159                                );
 160 $swq2 = Bio::Seq::Quality->new(-seq =>  'GT',
 161                                -qual => '65 0',
 162                                );
 163 my $swq3 = Bio::Seq::Quality->new(-seq =>  'AG',
 164                                -qual => '0 60',
 165                                );
 166 is $swq1->length, $swq2->length;
 167 is $swq1->length, $swq3->length;
 168
 169
 170 #
 171 # end of test inherited from seqwithquality.t
 172 #
 173 #################################################################
 174 #
 175 # testing new functionality
 176 #
 177
 178 my $qual = '0 1 2 3 4 5 6 7 8 9 11 12 13';
 179 my $trace = '0 5 10 15 20 25 30 35 40 45 50 55 60';
 180
 181 ok my $seq = Bio::Seq::Quality->new
 182     ( -qual => $qual,
 183       -trace_indices => $trace,
 184       -seq =>  'atcgatcgatcgt',
 185       -id  => 'human_id',
 186       -accession_number => 'S000012',
 187       -verbose => $DEBUG >= 0 ? $DEBUG : 0
 188 );
 189
 190
 191 print("2. Testing the trace() method...\n") if $DEBUG;
 192         print("\t2a) get\n") if $DEBUG;
 193         my @trace = @{$seq->trace()};
 194         my $str_trace = join(' ',@trace);
 195         is $str_trace, $trace;
 196         print("\t2b) set\n") if $DEBUG;
 197         ok $seq->trace("10 10 10 10 10");
 198         print("\t2c) get (again, to make sure the set was done.)\n") if $DEBUG;
 199         my @trace2 = @{$seq->trace()};
 200         my $str_trace2 = join(' ',@trace2);
 201         is($str_trace2, "10 10 10 10 10 0 0 0 0 0 0 0 0"); ###!
 202         print("\tSetting the trace back to the original value...\n") if $DEBUG;
 203         $seq->trace($trace);
 204
 205
 206
 207 is_deeply $seq->qual, [split / /, $qual];
 208 is_deeply $seq->trace, [split / /, $trace];
 209 is_deeply $seq->trace_indices, [split / /, $trace]; #deprecated
 210
 211 is $seq->qual_text, $qual;
 212 is $seq->trace_text, $trace;
 213
 214 is join (' ', @{$seq->subqual(2, 3)}), '1 2';
 215 is $seq->subqual_text(2, 3), '1 2';
 216 is join (' ', @{$seq->subqual(2, 3, "9 9")}), '9 9';
 217 is $seq->subqual_text(2, 3, "8 8"), '8 8';
 218
 219 is join (' ', @{$seq->subtrace(2, 3)}), '5 10';
 220 is $seq->subtrace_text(2, 3), '5 10';
 221 is join (' ', @{$seq->subtrace(2, 3, "9 9")}), '9 9';
 222 is $seq->subtrace_text(2, 3, "8 8"), '8 8';
 223
 224
 225 is $seq->trace_index_at(5), 20;
 226 is join(' ', @{$seq->sub_trace_index(5,6)}), "20 25";
 227
 228 is $seq->baseat(2), 't';
 229 is $seq->baseat(3), 'c';
 230 is $seq->baseat(4), 'g';
 231 is $seq->baseat(5), 'a';
 232
 233
 234 #############################################
 235 #
 236 # same tests using Seq::Meta::Array methods follow ...
 237 #
 238
 239 my $meta = '0 1 2 3 4 5 6 7 8 9 11 12';
 240 $trace = '0 5 10 15 20 25 30 35 40 45 50 55';
 241 my @trace_array = qw(0 5 10 15 20 25 30 35 40 45 50 55);
 242
 243 ok $seq = Bio::Seq::Quality->new
 244     ( -meta => $meta,
 245       -seq =>  'atcgatcgatcg',
 246       -id  => 'human_id',
 247       -accession_number => 'S000012',
 248       -verbose => $DEBUG >= 0 ? $DEBUG : 0
 249 );
 250
 251 $seq->named_meta('trace', \@trace_array);
 252
 253 is_deeply $seq->meta, [split / /, $meta];
 254 is_deeply $seq->named_meta('trace'), [split / /, $trace];
 255
 256 is $seq->meta_text, $meta;
 257 is $seq->named_meta_text('trace'), $trace;
 258
 259 is join (' ', @{$seq->submeta(2, 3)}), '1 2';
 260 is $seq->submeta_text(2, 3), '1 2';
 261 is join (' ', @{$seq->submeta(2, 3, "9 9")}), '9 9';
 262 is $seq->submeta_text(2, 3, "8 8"), '8 8';
 263
 264 is join (' ', @{$seq->named_submeta('trace', 2, 3)}), '5 10';
 265 is $seq->named_submeta_text('trace', 2, 3), '5 10';
 266 is join (' ', @{$seq->named_submeta('trace', 2, 3, "9 9")}), '9 9';
 267 is $seq->named_submeta_text('trace', 2, 3, "8 8"), '8 8';
 268
 269
 270 ok $seq = Bio::Seq::Quality->new(
 271     -seq => "ATGGGGGTGGTGGTACCCTATGGGGGTGGTGGTACCCT",
 272     -qual => "10 59 12 75 63 76 84 36 42 10 35 97 81 50 81 53 93 13 38 10 59 12 75 63 76 84 36 42 10 35 97 81 50 81 53 93 13 38",
 273     -trace_indices => "1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38"
 274 );
 275
 276 my $rev;
 277 ok $rev = $seq->revcom;
 278 is $rev->seq, 'AGGGTACCACCACCCCCATAGGGTACCACCACCCCCAT';
 279 is $rev->qual_text, "38 13 93 53 81 50 81 97 35 10 42 36 84 76 63 75 12 59 10 38 13 93 53 81 50 81 97 35 10 42 36 84 76 63 75 12 59 10";
 280
 281
 282 # selecting ranges based on quality
 283
 284 # test seq with three high quality regions (13, 12 and 3), one very short (3)
 285 ok $seq = Bio::Seq::Quality->new(
 286     -seq => "ATGGGGGTGGTGGTACCCTATGGGGGTGGTGGTACCCT",
 287     -qual => "0 5 10 20 30 40 40 50 50 50 50 50 40 10 10 10 5 5 20 20 30 40 50 44 44 50 50 50 50 50 5 5 40 40 40 40 50 50"
 288     );
 289
 290
 291 is $seq->threshold, undef;
 292 is $seq->threshold(10), 10;
 293 is $seq->threshold(13), 13;
 294
 295 is $seq->count_clear_ranges, 3;
 296
 297 my $newseq = $seq->get_clear_range;
 298 is $newseq->length, 12;
 299
 300
 301 my @ranges = $seq->get_all_clean_ranges;
 302 is scalar @ranges, 3;
 303 my $min_length = 10;
 304 @ranges = $seq->get_all_clean_ranges($min_length);
 305 is scalar @ranges, 2;
 306
 307 my $seqio = Bio::SeqIO->new(
 308     -file   => test_input_file('test_clear_range.fastq'),
 309     -format => 'fastq'
 310 );
 311
 312 while ( my $seq = $seqio->next_seq() ) {
 313     my $newqualobj;
 314     lives_ok { $newqualobj = $seq->get_clear_range(0) };
 315     if ($newqualobj) {
 316         is($newqualobj->id, $seq->id, 'Bug 2845');
 317     } else {
 318         ok(0, "No object returned via get_clear_range()");
 319     }
 320 }
 321
 322
 323
 324
 325 #############################################
 326 #
 327 # try testing some 'meta morphic relations'
 328 #
 329
 330 ## belief; As the threshold is increased, the number of clear ranges
 331 ## (ncr) should not decrease.
 332
 333 ## belief; As the thrshold is increased, the length of the clear
 334 ## ranges (lcr) should not decrease.
 335
 336 ## belief; As the threshold is incrazed, the clear range length (clr)
 337 ## should not increase. Sorry for the terribe var names.
 338
 339 ## belief; The number of clear ranges should vary between zero and
 340 ## half the sequence length.
 341
 342 ## belief; The length of the clear ranges should vary between zero and
 343 ## the sequence length.
 344
 345 ## belief; The length of the clear range should vary between zero and
 346 ## the sequence length.
 347
 348 ## belief; The lenght of the clear range should not be larger than the
 349 ## length of hte clear ranges.
 350
 351
 352 my @bases = qw (A T C G a t c g);
 353 my @qualities = 0..65;
 354
 355
 356 ## See beliefs above:
 357 my $ncr_thresh_sanity = 0;
 358 my $lcr_thresh_sanity = 0;
 359 my $clr_thresh_sanity = 0;
 360
 361 my $ncr_range_sanity = 0;
 362 my $lcr_range_sanity = 0;
 363 my $clr_range_sanity = 0;
 364
 365 my $final_loss_of_sanity = 0;
 366
 367
 368
 369 ## Go time:
 370
 371 for (1..100){
 372     $seq = join("", map {$bases[rand(@bases)]} 1..1000  );
 373     $qual = join(" ", map {$qualities[rand(@qualities)]} 1..1000  );
 374
 375     $seq = Bio::Seq::Quality->
 376         new(
 377             -seq => $seq,
 378             -qual => $qual,
 379             );
 380
 381     $seq->threshold(10);
 382     my $a_ncr = $seq->count_clear_ranges;
 383     my $a_lcr = $seq->clear_ranges_length;
 384     my $a_clr = scalar(@{$seq->get_clear_range->qual});
 385
 386     $ncr_range_sanity ++ if $a_ncr >= 0 && $a_ncr <=  500;
 387     $lcr_range_sanity ++ if $a_lcr >= 0 && $a_lcr <= 1000;
 388     $clr_range_sanity ++ if $a_clr >= 0 && $a_clr <= 1000;
 389     $final_loss_of_sanity ++ if $a_lcr >= $a_clr;
 390
 391     $seq->threshold(20);
 392     my $b_ncr = $seq->count_clear_ranges;
 393     my $b_lcr = $seq->clear_ranges_length;
 394     my $b_clr = scalar(@{$seq->get_clear_range->qual});
 395
 396     $ncr_range_sanity ++ if $b_ncr >= 0 && $b_ncr <=  500;
 397     $lcr_range_sanity ++ if $b_lcr >= 0 && $b_lcr <= 1000;
 398     $clr_range_sanity ++ if $b_clr >= 0 && $b_clr <= 1000;
 399     $final_loss_of_sanity ++ if $b_lcr >= $b_clr;
 400
 401
 402     $seq->threshold(30);
 403     my $c_ncr = $seq->count_clear_ranges;
 404     my $c_lcr = $seq->clear_ranges_length;
 405     my $c_clr = scalar(@{$seq->get_clear_range->qual});
 406
 407     $ncr_range_sanity ++ if $c_ncr >= 0 && $c_ncr <=  500;
 408     $lcr_range_sanity ++ if $c_lcr >= 0 && $c_lcr <= 1000;
 409     $clr_range_sanity ++ if $c_clr >= 0 && $c_clr <= 1000;
 410     $final_loss_of_sanity ++ if $c_lcr >= $c_clr;
 411
 412
 413
 414     $ncr_thresh_sanity ++ if
 415         $a_ncr <= $b_ncr &&
 416         $b_ncr <= $c_ncr;
 417
 418     $lcr_thresh_sanity ++ if
 419         $a_ncr <= $b_ncr &&
 420         $b_ncr <= $c_ncr;
 421
 422     $clr_thresh_sanity ++ if
 423         $a_clr >= $b_clr &&
 424         $b_clr >= $c_clr;
 425
 426 }
 427
 428 is $ncr_thresh_sanity, 100;
 429 is $lcr_thresh_sanity, 100;
 430 is $clr_thresh_sanity, 100;
 431
 432 is $ncr_range_sanity, 300;
 433 is $lcr_range_sanity, 300;
 434 is $clr_range_sanity, 300;
 435
 436 is $final_loss_of_sanity, 300;
 437
 438
 439
 440
 441 ## Test the mask sequence function ...
 442
 443 ## Ideally we'd at least test each function with each permutation of constructors.
 444
 445 my $x = Bio::Seq::Quality->
 446     new( -seq => "aaaattttccccgggg",
 447          -qual =>"1 1 1 1 2 2 2 2 1 1 1 1 3 3 3 3");
 448
 449 $x->threshold(1);
 450
 451 is $x->mask_below_threshold, "aaaattttccccgggg";
 452
 453 $x->threshold(2);
 454
 455 is $x->mask_below_threshold, "XXXXttttXXXXgggg";
 456
 457 $x->threshold(3);
 458
 459 is $x->mask_below_threshold, "XXXXXXXXXXXXgggg";
 460
 461 $x->threshold(4);
 462
 463 is $x->mask_below_threshold, "XXXXXXXXXXXXXXXX";
 464