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 => 60);
  11
  12     use_ok('Bio::Seq::Quality');
  13 }
  14
  15 my $DEBUG = test_debug();
  16
  17 # create some random sequence object with no id
  18 my $seqobj_broken = Bio::Seq::Quality->new( -seq => "ATCGATCGA",
  19                                           );
  20
  21 my $seqobj;
  22 lives_ok {
  23         $seqobj = Bio::Seq::Quality->new( -seq => "ATCGATCGA",
  24                                      -id  => 'QualityFragment-12',
  25                                      -accession_number => 'X78121',
  26                                    );
  27 };
  28
  29
  30 # create some random quality object with the same number of qualities and the same identifiers
  31 my $string_quals = "10 20 30 40 50 40 30 20 10";
  32 my $qualobj;
  33 lives_ok {
  34         $qualobj = Bio::Seq::Quality->new( -qual => $string_quals,
  35                                       -id  => 'QualityFragment-12',
  36                                       -accession_number => 'X78121',
  37                                         );
  38 };
  39
  40 # check to see what happens when you construct the Quality object
  41 ok my $swq1 = Bio::Seq::Quality->new( -seq => "ATCGATCGA",
  42                                       -id  => 'QualityFragment-12',
  43                                       -accession_number => 'X78121',
  44                                       -qual     =>      $string_quals);
  45
  46
  47 print("Testing various weird constructors...\n") if $DEBUG;
  48 print("\ta) No ids, Sequence object, no quality...\n") if $DEBUG;
  49 # w for weird
  50 my $wswq1;
  51 lives_ok {
  52         $wswq1 = Bio::Seq::Quality->new( -seq  =>  "ATCGATCGA",
  53                                          -qual  =>      "");
  54 };
  55 print $@ if $DEBUG;
  56
  57
  58 print("\tb) No ids, no sequence, quality object...\n") if $DEBUG;
  59         # note that you must provide a alphabet for this one.
  60 $wswq1 = Bio::Seq::Quality->new( -seq => "",
  61                                         -qual => $string_quals,
  62                                         -alphabet => 'dna'
  63 );
  64 print("\tc) Absolutely nothing. (HAHAHAHA)...\n") if $DEBUG;
  65 lives_ok {
  66         $wswq1 = Bio::Seq::Quality->new( -seq => "",
  67                                                 -qual => "",
  68                                                 -alphabet => 'dna'
  69         );
  70 };
  71
  72
  73 print("\td) Absolutely nothing but an ID\n") if $DEBUG;
  74 lives_ok {
  75     $wswq1 = Bio::Seq::Quality->new( -seq => "",
  76                                             -qual => "",
  77                                             -alphabet => 'dna',
  78                                             -id => 'an object with no sequence and no quality but with an id'
  79         );
  80 };
  81
  82 print("\td) No sequence, No quality, No ID...\n") if $DEBUG;
  83 warnings_like {
  84         $wswq1 = Bio::Seq::Quality->new( -seq  =>       "",
  85                                     -qual =>    "",
  86                                     -verbose => 0);
  87 } qr/Got a sequence with no letters in it cannot guess alphabet/;
  88
  89 print("Testing various methods and behaviors...\n") if $DEBUG;
  90
  91 print("1. Testing the seq() method...\n") if $DEBUG;
  92         print("\t1a) get\n") if $DEBUG;
  93         my $original_seq = $swq1->seq();
  94         is ($original_seq, "ATCGATCGA");
  95         print("\t1b) set\n") if $DEBUG;
  96         ok ($swq1->seq("AAAAAAAAAAAA"));
  97         print("\t1c) get (again, to make sure the set was done.)\n") if $DEBUG;
  98         is($swq1->seq(), "AAAAAAAAAAAA");
  99         print("\tSetting the sequence back to the original value...\n") if $DEBUG;
 100         $swq1->seq($original_seq);
 101
 102
 103 print("2. Testing the qual() method...\n") if $DEBUG;
 104         print("\t2a) get\n") if $DEBUG;
 105         my @qual = @{$swq1->qual()};
 106         my $str_qual = join(' ',@qual);
 107         is $str_qual, "10 20 30 40 50 40 30 20 10";
 108         print("\t2b) set\n") if $DEBUG;
 109         ok $swq1->qual("10 10 10 10 10");
 110         print("\t2c) get (again, to make sure the set was done.)\n") if $DEBUG;
 111         my @qual2 = @{$swq1->qual()};
 112         my $str_qual2 = join(' ',@qual2);
 113         is($str_qual2, "10 10 10 10 10 0 0 0 0"); ###!
 114         print("\tSetting the quality back to the original value...\n") if $DEBUG;
 115         $swq1->qual($str_qual);
 116
 117 print("3. Testing the length() method...\n") if $DEBUG;
 118         print("\t3a) When lengths are equal...\n") if $DEBUG;
 119         is($swq1->length(), 9);
 120         print("\t3b) When lengths are different\n") if $DEBUG;
 121         $swq1->qual("10 10 10 10 10");
 122         isnt ($swq1->length(), "DIFFERENT");
 123
 124
 125 print("6. Testing the subqual() method...\n") if $DEBUG;
 126      my $t_subqual = "10 20 30 40 50 60 70 80 90";
 127      $swq1->qual($t_subqual);
 128      print("\t6d) Testing the subqual at the start (border condition)\n") if $DEBUG;
 129           # ok ('1 2 3' eq join(' ',@{$swq1->subqual(1,3)}));
 130      print("\t6d) Testing the subqual at the end (border condition)\n") if $DEBUG;
 131           # ok ('7 8 9' eq join(' ',@{$swq1->subqual(7,9)}));
 132      print("\t6d) Testing the subqual in the middle\n") if $DEBUG;
 133           # ok ('4 5 6' eq join(' ',@{$swq1->subqual(4,6)}));
 134
 135 print("7. Testing cases where quality is zero...\n") if $DEBUG;
 136 $swq1 = Bio::Seq::Quality->new(-seq =>  'G',
 137                                -qual => '0',
 138                                      );
 139 my $swq2 = Bio::Seq::Quality->new(-seq =>  'G',
 140                                   -qual => '65',
 141                                      );
 142 is $swq1->length, $swq2->length;
 143
 144 $swq1 = Bio::Seq::Quality->new(-seq =>  'GC',
 145                                -qual => '0 0',
 146                                      );
 147 $swq2 = Bio::Seq::Quality->new(-seq =>  'GT',
 148                                -qual => '65 0',
 149                                      );
 150 is $swq1->length, $swq2->length;
 151
 152
 153 #
 154 # end of test inherited from seqwithquality.t
 155 #
 156 #################################################################
 157 #
 158 # testing new functionality
 159 #
 160
 161 my $qual = '0 1 2 3 4 5 6 7 8 9 11 12';
 162 my $trace = '0 5 10 15 20 25 30 35 40 45 50 55';
 163
 164 ok my $seq = Bio::Seq::Quality->new
 165     ( -qual => $qual,
 166       -trace_indices => $trace,
 167       -seq =>  'atcgatcgatcg',
 168       -id  => 'human_id',
 169       -accession_number => 'S000012',
 170       -verbose => $DEBUG >= 0 ? $DEBUG : 0
 171 );
 172
 173 is_deeply $seq->qual, [split / /, $qual];
 174 is_deeply $seq->trace, [split / /, $trace];
 175 is_deeply $seq->trace_indices, [split / /, $trace]; #deprecated
 176
 177 is $seq->qual_text, $qual;
 178 is $seq->trace_text, $trace;
 179
 180 is join (' ', @{$seq->subqual(2, 3)}), '1 2';
 181 is $seq->subqual_text(2, 3), '1 2';
 182 is join (' ', @{$seq->subqual(2, 3, "9 9")}), '9 9';
 183 is $seq->subqual_text(2, 3, "8 8"), '8 8';
 184
 185 is join (' ', @{$seq->subtrace(2, 3)}), '5 10';
 186 is $seq->subtrace_text(2, 3), '5 10';
 187 is join (' ', @{$seq->subtrace(2, 3, "9 9")}), '9 9';
 188 is $seq->subtrace_text(2, 3, "8 8"), '8 8';
 189
 190
 191 is $seq->trace_index_at(5), 20;
 192 is join(' ', @{$seq->sub_trace_index(5,6)}), "20 25";
 193
 194 is $seq->baseat(2), 't';
 195
 196
 197 #############################################
 198 #
 199 # same tests using Seq::Meta::Array methods follow ...
 200 #
 201
 202 my $meta = '0 1 2 3 4 5 6 7 8 9 11 12';
 203 $trace = '0 5 10 15 20 25 30 35 40 45 50 55';
 204 my @trace_array = qw(0 5 10 15 20 25 30 35 40 45 50 55);
 205
 206 ok $seq = Bio::Seq::Quality->new
 207     ( -meta => $meta,
 208       -seq =>  'atcgatcgatcg',
 209       -id  => 'human_id',
 210       -accession_number => 'S000012',
 211       -verbose => $DEBUG >= 0 ? $DEBUG : 0
 212 );
 213
 214 $seq->named_meta('trace', \@trace_array);
 215
 216 is_deeply $seq->meta, [split / /, $meta];
 217 is_deeply $seq->named_meta('trace'), [split / /, $trace];
 218
 219 is $seq->meta_text, $meta;
 220 is $seq->named_meta_text('trace'), $trace;
 221
 222 is join (' ', @{$seq->submeta(2, 3)}), '1 2';
 223 is $seq->submeta_text(2, 3), '1 2';
 224 is join (' ', @{$seq->submeta(2, 3, "9 9")}), '9 9';
 225 is $seq->submeta_text(2, 3, "8 8"), '8 8';
 226
 227 is join (' ', @{$seq->named_submeta('trace', 2, 3)}), '5 10';
 228 is $seq->named_submeta_text('trace', 2, 3), '5 10';
 229 is join (' ', @{$seq->named_submeta('trace', 2, 3, "9 9")}), '9 9';
 230 is $seq->named_submeta_text('trace', 2, 3, "8 8"), '8 8';
 231
 232
 233 ok $seq = Bio::Seq::Quality->new(
 234     -seq => "ATGGGGGTGGTGGTACCCTATGGGGGTGGTGGTACCCT",
 235     -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",
 236     -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"
 237 );
 238
 239 my $rev;
 240 ok $rev = $seq->revcom;
 241 is $rev->seq, 'AGGGTACCACCACCCCCATAGGGTACCACCACCCCCAT';
 242 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";
 243
 244
 245 # selecting ranges based on quality
 246
 247 # test seq with three high quality regions (13, 12 and 3), one very short (3)
 248 ok $seq = Bio::Seq::Quality->new(
 249     -seq => "ATGGGGGTGGTGGTACCCTATGGGGGTGGTGGTACCCT",
 250     -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"
 251     );
 252
 253
 254 is $seq->threshold, undef;
 255 is $seq->threshold(10), 10;
 256 is $seq->threshold(13), 13;
 257
 258 is $seq->count_clear_ranges, 3;
 259
 260 my $newseq = $seq->get_clear_range;
 261 is $newseq->length, 12;
 262
 263
 264 my @ranges = $seq->get_all_clean_ranges;
 265 is scalar @ranges, 3;
 266 my $min_length = 10;
 267 @ranges = $seq->get_all_clean_ranges($min_length);
 268 is scalar @ranges, 2;