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