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