t/Seq/SimulatedRead.t

   1 # -*-Perl-*- Test Harness script for Bioperl
   2 # $Id$
   3
   4 use strict;
   5 use warnings;
   6
   7 BEGIN {
   8     use Bio::Root::Test;
   9     test_begin(-tests => 194);
  10
  11     use_ok('Bio::Seq');
  12     use_ok('Bio::Seq::Quality');
  13     use_ok('Bio::PrimarySeq');
  14     use_ok('Bio::LocatableSeq');
  15     use_ok('Bio::Seq::SimulatedRead');
  16 }
  17
  18 my $VERBOSE = test_debug();
  19
  20 my ($ref, $ref2, $ref3, $ref4, $ref5, $read, $errors);
  21
  22 $ref = Bio::Seq::Quality->new(-id    => 'human_id',
  23                                -seq   => 'TAAAAAAACCCC',
  24                                -qual  => '1 2 3 4 5 6 7 8 9 10 11 12',
  25                                -trace => '0 5 10 15 20 25 30 35 40 45 50 55',
  26                                -desc  => 'The human genome' );
  27
  28 $ref2 = Bio::Seq->new(-id   => 'other_genome',
  29                        -seq  => 'ACGTACGT',
  30                        -desc => '"Secret" sequence');
  31
  32 $ref3 = Bio::PrimarySeq->new(-seq => 'ACACTGATCTAGCGTCGTGCTAGCTGACGTAGCTGAT' );
  33
  34 $ref4 = Bio::LocatableSeq->new(-id  => 'a_thaliana',
  35                                 -seq => 'CGTATTCTGAGGAGAGCTCT' );
  36
  37
  38 # Basic object
  39
  40 ok $read = Bio::Seq::SimulatedRead->new();
  41 isa_ok $read, 'Bio::Seq::SimulatedRead';
  42 isa_ok $read, 'Bio::LocatableSeq';
  43 isa_ok $read, 'Bio::Seq::Quality';
  44
  45 $errors->{'1'}->{'+'} = 'T';
  46 ok $read = Bio::Seq::SimulatedRead->new(-reference => $ref, -errors => $errors );
  47 is $read->reference, $ref;
  48 ok $read->errors;
  49 is $read->errors->{'1'}->{'+'}->[0], 'T';
  50
  51 ok $read = Bio::Seq::SimulatedRead->new(-reference => $ref, -track => 0 );
  52 is $read->start, 1;
  53 is $read->end, 12;
  54 is $read->seq, 'TAAAAAAACCCC';
  55 is $read->track, 0;
  56 is $read->desc, undef;
  57 is $read->revcom->seq, 'GGGGTTTTTTTA';
  58
  59 ok $read = Bio::Seq::SimulatedRead->new(-reference => $ref, -track => 1 );
  60 is $read->start, 1;
  61 is $read->end, 12;
  62 is $read->seq, 'TAAAAAAACCCC';
  63 is join(' ',@{$read->qual}), '';
  64 is $read->track, 1;
  65 is $read->desc, 'reference=human_id start=1 end=12 strand=+1 description="The human genome"';
  66
  67 ok $read = Bio::Seq::SimulatedRead->new(-reference => $ref, -track => 1, -coord_style => 'bioperl' );
  68 is $read->desc, 'reference=human_id start=1 end=12 strand=+1 description="The human genome"';
  69
  70 ok $read = Bio::Seq::SimulatedRead->new(-reference => $ref, -track => 1, -coord_style => 'genbank' );
  71 is $read->desc, 'reference=human_id position=1..12 description="The human genome"';
  72
  73 ok $read = Bio::Seq::SimulatedRead->new(-reference => $ref, -qual_levels => [30, 10]);
  74 is $read->start, 1;
  75 is $read->end, 12;
  76 is $read->seq, 'TAAAAAAACCCC';
  77 is join(' ', @{$read->qual}), '30 30 30 30 30 30 30 30 30 30 30 30';
  78 is $read->track, 1;
  79 is $read->desc, 'reference=human_id start=1 end=12 strand=+1 description="The human genome"';
  80 is $read->revcom->seq, 'GGGGTTTTTTTA';
  81
  82 ok $read = Bio::Seq::SimulatedRead->new(-reference => $ref2 );
  83 is $read->start, 1;
  84 is $read->end, 8;
  85 is $read->seq, 'ACGTACGT';
  86 is join(' ',@{$read->qual}), '';
  87 is $read->desc, 'reference=other_genome start=1 end=8 strand=+1 description="\"Secret\" sequence"';
  88
  89 ok $read = Bio::Seq::SimulatedRead->new(-reference => $ref3 );
  90 is $read->start, 1;
  91 is $read->end, 37;
  92 is $read->seq, 'ACACTGATCTAGCGTCGTGCTAGCTGACGTAGCTGAT';
  93 is join(' ',@{$read->qual}), '';
  94 is $read->desc, 'start=1 end=37 strand=+1';
  95
  96 ok $read = Bio::Seq::SimulatedRead->new(-reference => $ref, -strand => -1, -qual_levels => [30, 10]);
  97 is $read->seq, 'GGGGTTTTTTTA';
  98 is join(' ', @{$read->qual}), '30 30 30 30 30 30 30 30 30 30 30 30';
  99 is $read->desc, 'reference=human_id start=1 end=12 strand=-1 description="The human genome"';
 100
 101 ok $read = Bio::Seq::SimulatedRead->new(-reference => $ref, -strand => -1, -qual_levels => [30, 10], -coord_style => 'genbank' );
 102 is $read->desc, 'reference=human_id position=complement(1..12) description="The human genome"';
 103
 104 ok $read = Bio::Seq::SimulatedRead->new(-reference => $ref, -start => 2, -end => 8, -qual_levels => [30, 10]);
 105 is $read->seq, 'AAAAAAA';
 106 is join(' ', @{$read->qual}), '30 30 30 30 30 30 30';
 107 is $read->desc, 'reference=human_id start=2 end=8 strand=+1 description="The human genome"';
 108
 109 ok $read = Bio::Seq::SimulatedRead->new(-reference => $ref, -strand => -1, -start => 2, -end => 8, -qual_levels => [30, 10]);
 110 is $read->seq, 'TTTTTTT';
 111 is join(' ', @{$read->qual}), '30 30 30 30 30 30 30';
 112 is $read->desc, 'reference=human_id start=2 end=8 strand=-1 description="The human genome"';
 113
 114 $errors = {};
 115 $errors->{'6'}->{'+'} = 'GG';
 116 ok $read = Bio::Seq::SimulatedRead->new(-reference => $ref, -strand => -1, -start => 2, -end => 8, -errors => $errors, -qual_levels => [30, 10]);
 117 is $read->start, 2;
 118 is $read->end, 8;
 119 is $read->seq, 'TTTTTTGGT';
 120 is join(' ', @{$read->qual}), '30 30 30 30 30 30 10 10 30';
 121 is $read->desc, 'reference=human_id start=2 end=8 strand=-1 errors=6+G,6+G description="The human genome"';
 122
 123 $errors = {};
 124 $errors->{'6'}->{'+'} = 'GG';
 125 $errors->{'1'}->{'%'} = 'T';
 126 $errors->{'3'}->{'-'} = undef;
 127 ok $read = Bio::Seq::SimulatedRead->new(-reference => $ref, -strand => 1, -start => 2, -end => 8, -errors => $errors, -qual_levels => [30, 10]);
 128 is $read->start, 2;
 129 is $read->end, 8;
 130 is $read->seq, 'TAAAAGGA';
 131 is join(' ', @{$read->qual}), '10 30 30 30 30 10 10 30';
 132 is $read->desc, 'reference=human_id start=2 end=8 strand=+1 errors=1%T,3-,6+G,6+G description="The human genome"';
 133
 134 $errors = {};
 135 $errors->{'6'}->{'+'} = 'GG';
 136 ok $read = Bio::Seq::SimulatedRead->new(-reference => $ref, -errors => $errors, -qual_levels => [30, 10]);
 137 is $read->start, 1;
 138 is $read->end, 12;
 139 is $read->seq, 'TAAAAAGGAACCCC';
 140 is join(' ', @{$read->qual}), '30 30 30 30 30 30 10 10 30 30 30 30 30 30';
 141 is $read->desc, 'reference=human_id start=1 end=12 strand=+1 errors=6+G,6+G description="The human genome"';
 142
 143 ok $read = Bio::Seq::SimulatedRead->new(-reference => $ref, -errors => $errors, -mid => 'ACGT', -errors => $errors, -qual_levels => [30, 10]);
 144 is $read->start, 1;
 145 is $read->end, 12;
 146 is $read->seq, 'ACGTTAGGAAAAAACCCC';
 147 is join(' ', @{$read->qual}), '30 30 30 30 30 30 10 10 30 30 30 30 30 30 30 30 30 30';
 148 is $read->desc, 'reference=human_id start=1 end=12 strand=+1 mid=ACGT errors=6+G,6+G description="The human genome"';
 149
 150 ok $read = Bio::Seq::SimulatedRead->new(-reference => $ref, -mid => 'TTTAAA', -qual_levels => [30, 10]);
 151 is $read->start, 1;
 152 is $read->end, 12;
 153 is $read->seq, 'TTTAAATAAAAAAACCCC';
 154 is join(' ', @{$read->qual}), '30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30';
 155 is $read->desc, 'reference=human_id start=1 end=12 strand=+1 mid=TTTAAA description="The human genome"';
 156
 157 ok $read = Bio::Seq::SimulatedRead->new(-reference => $ref, -mid => '', -qual_levels => []);
 158 is $read->start, 1;
 159 is $read->end, 12;
 160 is $read->seq, 'TAAAAAAACCCC';
 161 is join(' ', @{$read->qual}), '';
 162 is $read->desc, 'reference=human_id start=1 end=12 strand=+1 description="The human genome"';
 163
 164
 165 # Redundant errors
 166
 167 $errors = {};
 168 $errors->{'6'}->{'+'} = ['G', 'G'];
 169 $errors->{'1'}->{'%'} = ['A', 'G', 'T'];
 170 $errors->{'3'}->{'-'} = [undef, undef];
 171 ok $read = Bio::Seq::SimulatedRead->new(-reference => $ref, -strand => 1, -start => 2, -end => 8, -errors => $errors, -qual_levels => [30, 10]), 'redundant errors';
 172 is $read->start, 2;
 173 is $read->end, 8;
 174 is $read->seq, 'TAAAAGGA';
 175 is join(' ', @{$read->qual}), '10 30 30 30 30 10 10 30';
 176 is $read->desc, 'reference=human_id start=2 end=8 strand=+1 errors=1%A,1%G,1%T,3-,3-,6+G,6+G description="The human genome"';
 177
 178
 179 # Specifying errors() after new()
 180
 181 ok $read = Bio::Seq::SimulatedRead->new(-reference => $ref, -mid => '', -qual_levels => []);
 182 is $read->start, 1;
 183 is $read->end, 12;
 184 is $read->seq, 'TAAAAAAACCCC';
 185 is join(' ', @{$read->qual}), '';
 186 is $read->desc, 'reference=human_id start=1 end=12 strand=+1 description="The human genome"';
 187
 188 $errors = {};
 189 ok $read->errors($errors), 'errors()';
 190 is $read->start, 1;
 191 is $read->end, 12;
 192 is $read->seq, 'TAAAAAAACCCC';
 193 is $read->desc, 'reference=human_id start=1 end=12 strand=+1 description="The human genome"';
 194
 195 $errors = {};
 196 $errors->{'6'}->{'+'} = 'GG';
 197 ok $read->errors($errors);
 198 is $read->seq, 'TAAAAAGGAACCCC';
 199 is $read->start, 1;
 200 is $read->end, 12;
 201 is $read->desc, 'reference=human_id start=1 end=12 strand=+1 errors=6+G,6+G description="The human genome"';
 202
 203
 204 # More tracking tests
 205
 206 ok $read = Bio::Seq::SimulatedRead->new(-reference => $ref, -mid => 'ACGT', -qual_levels => [], -coord_style => 'genbank' );
 207 is $read->desc, 'reference=human_id position=1..12 mid=ACGT description="The human genome"';
 208
 209 ok $read->mid('AAAA');
 210 is $read->desc, 'reference=human_id position=1..12 mid=AAAA description="The human genome"';
 211
 212 $errors = {};
 213 $errors->{'6'}->{'+'} = 'GG';
 214 ok $read->errors($errors);
 215 is $read->desc, 'reference=human_id position=1..12 mid=AAAA errors=6+G,6+G description="The human genome"';
 216
 217 ok not($read->track(0)), 'track()';
 218 is $read->track, 0;
 219 is $read->desc, undef;
 220 ok $read->track(1);
 221 is $read->track, 1;
 222 is $read->desc, 'reference=human_id position=1..12 mid=AAAA errors=6+G,6+G description="The human genome"';
 223
 224
 225 # qual_levels() method
 226
 227 ok $read = Bio::Seq::SimulatedRead->new(-verbose => $VERBOSE, );
 228 ok $read->qual_levels([30, 10]), 'qual_levels()';
 229 is join(' ', @{$read->qual_levels}), '30 10';
 230
 231 # reference() method
 232
 233 ok $read->reference($ref), 'reference()';
 234 is $read->reference(), $ref;
 235
 236 # mid() method
 237
 238 ok $read = Bio::Seq::SimulatedRead->new(-verbose => $VERBOSE, ), 'mid()';
 239 ok $read->qual_levels([30, 10]);
 240 ok $read->reference($ref);
 241 ok $read->mid('ACGT');
 242 ok $read->mid, 'ACGT';
 243
 244 is $read->seq, 'ACGTTAAAAAAACCCC';
 245 is join(' ', @{$read->qual}), '30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30';
 246 is $read->desc, 'reference=human_id start=1 end=12 strand=+1 mid=ACGT description="The human genome"';
 247
 248 ok $read->mid('TTTAAA');
 249 ok $read->mid, 'TTTAAA';
 250 is $read->seq, 'TTTAAATAAAAAAACCCC';
 251 is join(' ', @{$read->qual}), '30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30';
 252 is $read->desc, 'reference=human_id start=1 end=12 strand=+1 mid=TTTAAA description="The human genome"';
 253
 254
 255 # Edge case... mutation of the last bases of a simulated read with MID
 256
 257 $errors = {};
 258 $errors->{'18'}->{'%'} = 'T';
 259 $read->errors($errors);
 260 is $read->seq, 'TTTAAATAAAAAAACCCT';
 261
 262
 263 # Try different BioPerl object types
 264
 265 ok $read = Bio::Seq::SimulatedRead->new(-reference => $ref  ), 'Bio::Seq::Quality';
 266 ok $read = Bio::Seq::SimulatedRead->new(-reference => $ref2 ), 'Bio::Seq';
 267 ok $read = Bio::Seq::SimulatedRead->new(-reference => $ref3 ), 'Bio::PrimarySeq';
 268 ok $read = Bio::Seq::SimulatedRead->new(-reference => $ref4 ), 'Bio::LocatableSeq';
 269
 270
 271 # More detailed tests of the error specifications
 272
 273 $errors = {};
 274 $errors->{'0'}->{'-'} = undef;
 275 warning_like {$read = Bio::Seq::SimulatedRead->new(-reference => $ref, -errors => $errors )}
 276     qr/Positions of substitutions and deletions have to be strictly positive but got 0/;
 277 is $read->seq, 'TAAAAAAACCCC';
 278
 279 $errors = {};
 280 $errors->{'1'}->{'-'} = undef;
 281 ok $read = Bio::Seq::SimulatedRead->new(-reference => $ref, -errors => $errors );
 282 is $read->seq, 'AAAAAAACCCC';
 283
 284 $errors = {};
 285 $errors->{'12'}->{'-'} = undef;
 286 ok $read = Bio::Seq::SimulatedRead->new(-reference => $ref, -errors => $errors );
 287 is $read->seq, 'TAAAAAAACCC';
 288
 289 $errors = {};
 290 $errors->{'13'}->{'-'} = undef;
 291 warning_like {$read = Bio::Seq::SimulatedRead->new(-reference => $ref, -errors => $errors )}
 292     qr/Position 13 is beyond end of read \(12 residues\)/; # there should be a warning too
 293 is $read->seq, 'TAAAAAAACCCC';
 294
 295 $errors = {};
 296 $errors->{'0'}->{'%'} = 'G';
 297 warning_like {$read = Bio::Seq::SimulatedRead->new(-reference => $ref, -errors => $errors )}
 298     qr/Positions of substitutions and deletions have to be strictly positive/; # there should be a warning too
 299 is $read->seq, 'TAAAAAAACCCC';
 300
 301 $errors = {};
 302 $errors->{'1'}->{'%'} = 'G';
 303 ok $read = Bio::Seq::SimulatedRead->new(-reference => $ref, -errors => $errors );
 304 is $read->seq, 'GAAAAAAACCCC';
 305
 306 $errors = {};
 307 $errors->{'12'}->{'%'} = 'G';
 308 ok $read = Bio::Seq::SimulatedRead->new(-reference => $ref, -errors => $errors );
 309 is $read->seq, 'TAAAAAAACCCG';
 310
 311 $errors = {};
 312 $errors->{'13'}->{'%'} = 'G';
 313 warning_like { $read = Bio::Seq::SimulatedRead->new(-reference => $ref, -errors => $errors ) }
 314     qr/Position 13 is beyond end of read \(12 residues\)/; # there should be a warning too
 315 is $read->seq, 'TAAAAAAACCCC';
 316
 317 $errors = {};
 318 $errors->{'0'}->{'+'} = 'A';
 319 ok $read = Bio::Seq::SimulatedRead->new(-reference => $ref, -errors => $errors );
 320 is $read->seq, 'ATAAAAAAACCCC';
 321
 322 $errors = {};
 323 $errors->{'1'}->{'+'} = 'A';
 324 ok $read = Bio::Seq::SimulatedRead->new(-reference => $ref, -errors => $errors );
 325 is $read->seq, 'TAAAAAAAACCCC';
 326
 327 $errors = {};
 328 $errors->{'12'}->{'+'} = 'A';
 329 ok $read = Bio::Seq::SimulatedRead->new(-reference => $ref, -errors => $errors );
 330 is $read->seq, 'TAAAAAAACCCCA';
 331
 332 $errors = {};
 333 $errors->{'13'}->{'+'} = 'A';
 334 warning_like {$read = Bio::Seq::SimulatedRead->new(-reference => $ref, -errors => $errors ) }
 335     qr/Position 13 is beyond end of read \(12 residues\)/; # there should be a warning too
 336 is $read->seq, 'TAAAAAAACCCC';
 337
 338 $errors = {};
 339 $errors->{'1'}->{'%'} = 'G';
 340 $errors->{'2'}->{'%'} = 'G';
 341 $errors->{'3'}->{'%'} = 'G';
 342 ok $read = Bio::Seq::SimulatedRead->new(-reference => $ref, -errors => $errors );
 343 is $read->seq, 'GGGAAAAACCCC';
 344
 345 $errors = {};
 346 $errors->{'1'}->{'+'} = 'G';
 347 $errors->{'2'}->{'+'} = 'G';
 348 $errors->{'3'}->{'+'} = 'G';
 349 ok $read = Bio::Seq::SimulatedRead->new(-reference => $ref, -errors => $errors );
 350 is $read->seq, 'TGAGAGAAAAACCCC';
 351
 352 $errors = {};
 353 $errors->{'1'}->{'-'} = undef;
 354 $errors->{'2'}->{'-'} = undef;
 355 $errors->{'3'}->{'-'} = undef;
 356 ok $read = Bio::Seq::SimulatedRead->new(-reference => $ref, -errors => $errors );
 357 is $read->seq, 'AAAAACCCC';
 358
 359 $errors = {};
 360 $errors->{'1'}->{'+'} = 'GGG';
 361 $errors->{'2'}->{'-'} = undef;
 362 ok $read = Bio::Seq::SimulatedRead->new(-reference => $ref, -errors => $errors );
 363 is $read->seq, 'TGGGAAAAAACCCC';
 364
 365 $errors = {};
 366 $errors->{'2'}->{'+'} = 'CC';
 367 $errors->{'2'}->{'-'} = undef;
 368 ok $read = Bio::Seq::SimulatedRead->new(-reference => $ref, -errors => $errors );
 369 is $read->seq, 'TCCAAAAAACCCC';
 370
 371 $errors = {};
 372 $errors->{'2'}->{'%'} = 'C';
 373 $errors->{'2'}->{'-'} = undef;
 374 ok $read = Bio::Seq::SimulatedRead->new(-reference => $ref, -errors => $errors );
 375 is $read->seq, 'TAAAAAACCCC';