t/Assembly/ContigSpectrum.t

   1 use strict;
   2
   3 BEGIN {
   4     use lib '.';
   5     use Bio::Root::Test;
   6
   7     test_begin( -tests            => 236,
   8                 -requires_modules => [ qw(Bio::Assembly::Tools::ContigSpectrum)] );
   9     use_ok('Bio::Assembly::IO');
  10     use_ok('Bio::Assembly::Tools::ContigSpectrum');
  11 }
  12
  13 my $verbose = test_debug();
  14
  15 my $in = Bio::Assembly::IO->new(
  16   -file => test_input_file('contigspectrumtest.tigr'),
  17   -format => 'tigr'
  18 );
  19 isa_ok($in, 'Bio::Assembly::IO');
  20 my $sc = $in->next_assembly;
  21 isa_ok($sc, 'Bio::Assembly::Scaffold');
  22
  23 # Try all the get/set methods
  24 ok(my $csp = Bio::Assembly::Tools::ContigSpectrum->new, 'get/set methods');
  25 isa_ok($csp, 'Bio::Assembly::Tools::ContigSpectrum');
  26 ok($csp->id('asdf'));
  27 is($csp->id, 'asdf');
  28 ok($csp->nof_seq(123));
  29 is($csp->nof_seq, 123);
  30 ok($csp->nof_rep(456));
  31 is($csp->nof_rep, 456);
  32 ok($csp->max_size(789));
  33 is($csp->max_size, 789);
  34 ok($csp->nof_overlaps(111));
  35 is($csp->nof_overlaps, 111);
  36 ok($csp->min_overlap(50));
  37 is($csp->min_overlap, 50);
  38 ok($csp->avg_overlap(54.3));
  39 is($csp->avg_overlap, 54.3);
  40 ok($csp->min_identity(89.1));
  41 is($csp->min_identity, 89.1);
  42 ok($csp->avg_identity(98.7));
  43 is($csp->avg_identity, 98.7);
  44 ok($csp->avg_seq_len(123.456));
  45 is($csp->avg_seq_len, 123.456);
  46 ok($csp->eff_asm_params(1));
  47 is($csp->eff_asm_params, 1);
  48
  49 # contig spectrum based on simple spectrum
  50 ok(my $spectrum_csp = Bio::Assembly::Tools::ContigSpectrum->new, 'simple spectrum');
  51 ok($spectrum_csp->spectrum({1=>1, 2=>2, 3=>3}));
  52 is($spectrum_csp->eff_asm_params, 0);
  53 is($spectrum_csp->nof_seq, 14);
  54 is($spectrum_csp->max_size, 3);
  55 is($spectrum_csp->nof_rep, 1);
  56 is($spectrum_csp->nof_overlaps, 0);
  57 is($spectrum_csp->min_overlap, undef);
  58 is($spectrum_csp->avg_overlap, 0);
  59 is($spectrum_csp->min_identity, undef);
  60 is($spectrum_csp->avg_identity, 0);
  61 is($spectrum_csp->avg_seq_len, 0);
  62 is(scalar $spectrum_csp->assembly, 0);
  63
  64 ok(my $string = $spectrum_csp->to_string(1));
  65 is($string, '1 2 3');
  66 ok($string = $spectrum_csp->to_string(2));
  67 is($string, "1\t2\t3");
  68 ok($string = $spectrum_csp->to_string(3));
  69 is($string, "1\n2\n3");
  70
  71 # mixed contig spectrum imported from assembly
  72 ok(my $mixed_csp = Bio::Assembly::Tools::ContigSpectrum->new(
  73   -assembly       => $sc,
  74   -eff_asm_params => 1 ), 'mixed contig spectrum');
  75 is_deeply($mixed_csp->spectrum, {1=>0, 2=>3, 6=>1, 9=>1}); # [0 3 0 0 0 1 0 0 1]
  76 is($mixed_csp->eff_asm_params, 1);
  77 is($mixed_csp->max_size, 9);
  78 is($mixed_csp->nof_rep, 1);
  79 is($mixed_csp->nof_seq, 21);
  80 float_is($mixed_csp->avg_seq_len, 303.81);
  81 is($mixed_csp->nof_overlaps, 16);
  82 is($mixed_csp->min_overlap, 35);
  83 is($mixed_csp->avg_overlap, 155.875);
  84 float_is($mixed_csp->min_identity, 96.8421);
  85 float_is($mixed_csp->avg_identity, 98.8826);
  86 is(scalar $mixed_csp->assembly, 1);
  87
  88 # dissolved contig spectrum
  89 ok(my $dissolved_csp = Bio::Assembly::Tools::ContigSpectrum->new(
  90   -dissolve => [$mixed_csp, 'ZZZ'] ), 'dissolved contig spectrum');
  91 is_deeply($dissolved_csp->spectrum, {1=>2, 2=>1}); # [2 1]
  92 is($dissolved_csp->eff_asm_params, 0);
  93 is($dissolved_csp->max_size, 2);
  94 is($dissolved_csp->nof_rep, 1);
  95 is($dissolved_csp->nof_seq, 4);
  96 float_is($dissolved_csp->avg_seq_len, 321);
  97 # eff_asm_params haven't been requested
  98 is($dissolved_csp->nof_overlaps, 0);
  99 is($dissolved_csp->min_overlap, undef);
 100 is($dissolved_csp->avg_overlap, 0);
 101 is($dissolved_csp->min_identity, undef);
 102 is($dissolved_csp->avg_identity, 0);
 103
 104 ok($dissolved_csp = Bio::Assembly::Tools::ContigSpectrum->new(
 105   -dissolve => [$mixed_csp, 'sdsu'] ));
 106 is_deeply($dissolved_csp->spectrum, {1=>3, 6=>1}); # [3 0 0 0 0 1]
 107 is($dissolved_csp->eff_asm_params, 0);
 108 is($dissolved_csp->max_size, 6);
 109 is($dissolved_csp->nof_rep, 1);
 110 is($dissolved_csp->nof_seq, 9);
 111 float_is($dissolved_csp->avg_seq_len, 441.222222222222);
 112 # eff_asm_params haven't been requested
 113 is($dissolved_csp->nof_overlaps, 0);
 114 is($dissolved_csp->min_overlap, undef);
 115 is($dissolved_csp->avg_overlap, 0);
 116 is($dissolved_csp->min_identity, undef);
 117 is($dissolved_csp->avg_identity, 0);
 118
 119 ok($dissolved_csp = Bio::Assembly::Tools::ContigSpectrum->new(
 120   -dissolve => [$mixed_csp, 'ABC'] ));
 121 is_deeply($dissolved_csp->spectrum, {1=>2, 6=>1}); # [2 0 0 0 0 1]
 122 is($dissolved_csp->eff_asm_params, 0);
 123 is($dissolved_csp->max_size, 6);
 124 is($dissolved_csp->nof_rep, 1);
 125 is($dissolved_csp->nof_seq, 8);
 126 is($dissolved_csp->avg_seq_len, 140.625);
 127 # eff_asm_params haven't been requested
 128 is($dissolved_csp->nof_overlaps, 0);
 129 is($dissolved_csp->min_overlap, undef);
 130 is($dissolved_csp->avg_overlap, 0);
 131 is($dissolved_csp->min_identity, undef);
 132 is($dissolved_csp->avg_identity, 0);
 133
 134 ok($dissolved_csp = Bio::Assembly::Tools::ContigSpectrum->new(
 135   -min_overlap  => 62,
 136   -min_identity => 1,
 137   -dissolve     => [$mixed_csp, 'ABC'] ));
 138 is_deeply($dissolved_csp->spectrum, {1=>2, 6=>1}); # [2 0 0 0 0 1]
 139
 140 ok($dissolved_csp = Bio::Assembly::Tools::ContigSpectrum->new(
 141    -min_overlap  => 63,
 142    -min_identity => 1,
 143    -dissolve     => [$mixed_csp, 'ABC'] ));
 144 is_deeply($dissolved_csp->spectrum, {1=>3, 5=>1}); # [3 0 0 0 1]
 145
 146 # after dissolving, the remaining assembly objects should be 3 singlets and 1 6-contig
 147 my @contigs = ($dissolved_csp->assembly);
 148 is(scalar @contigs, 4);
 149 my @contig_ids = sort qw( 144 652_1 652_2 652_3 );
 150 is_deeply( [sort map($_->id, @contigs)], \@contig_ids );
 151 my @contig_sizes = sort qw( 1 1 1 5 );
 152 is_deeply( [sort map($_->num_sequences, @contigs)], \@contig_sizes );
 153 my @contig_isas = sort qw( Bio::Assembly::Singlet Bio::Assembly::Singlet
 154 Bio::Assembly::Singlet Bio::Assembly::Contig );
 155 is_deeply( [sort map(ref $_, @contigs)], \@contig_isas );
 156 my @reads = ($contigs[1])->each_seq;
 157 my @read_ids = sort qw(ABC|9980040 ABC|9937790 ABC|9956706 ABC|9960711 ABC|9976538);
 158 is_deeply( [sort map($_->id, @reads)], \@read_ids );
 159
 160 ok($dissolved_csp = Bio::Assembly::Tools::ContigSpectrum->new(
 161   -min_overlap  => 62,
 162   -min_identity => 97,
 163   -dissolve     => [$mixed_csp, 'ABC'] ));
 164 is_deeply($dissolved_csp->spectrum, {1=>2, 6=>1}); # [2 0 0 0 0 1]
 165
 166 ok($dissolved_csp = Bio::Assembly::Tools::ContigSpectrum->new(
 167    -min_overlap  => 62,
 168    -min_identity => 98,
 169    -dissolve     => [$mixed_csp, 'ABC'] ));
 170 is_deeply($dissolved_csp->spectrum, {1=>2, 6=>1}); # [2 0 0 0 0 1]
 171
 172 ok($dissolved_csp = Bio::Assembly::Tools::ContigSpectrum->new(
 173   -dissolve       => [$mixed_csp, 'ABC'],
 174   -eff_asm_params => 1 ));
 175 is_deeply($dissolved_csp->spectrum, {1=>2, 6=>1}); # [2 0 0 0 0 1]
 176 is($dissolved_csp->eff_asm_params, 1);
 177 is($dissolved_csp->max_size, 6);
 178 is($dissolved_csp->nof_rep, 1);
 179 is($dissolved_csp->nof_seq, 8);
 180 float_is($dissolved_csp->avg_seq_len, 140.625);
 181 is($dissolved_csp->nof_overlaps, 5);
 182 float_is($dissolved_csp->avg_overlap, 76.8);
 183 float_is($dissolved_csp->avg_identity, 100.0);
 184 # min_overlap and min_identity not explicitely specified for the dissolved csp
 185 # min_overlap and min_identity are thus taken from the mixed csp
 186 is($dissolved_csp->min_overlap, 35);
 187 float_is($dissolved_csp->min_identity, 96.8421);
 188
 189 # cross contig spectrum
 190 ok(my $cross_csp = Bio::Assembly::Tools::ContigSpectrum->new(
 191   -cross => $mixed_csp), 'cross-contig spectrum');
 192 is_deeply($cross_csp->spectrum, {1=>7, 2=>2, 9=>1}); # [7 2 0 0 0 0 0 0 1]
 193
 194 # assembly should have 2 singlets and 1 9-contig
 195 @contigs = $cross_csp->assembly;
 196 is(scalar @contigs, 3);
 197 @contig_sizes = sort qw( 2 2 9 );
 198 is_deeply( [sort map($_->num_sequences, @contigs)], \@contig_sizes );
 199 @contig_isas = sort qw( Bio::Assembly::Contig Bio::Assembly::Contig Bio::Assembly::Contig);
 200 is_deeply( [sort map(ref $_, @contigs)], \@contig_isas );
 201 @read_ids = sort qw(sdsu|SDSU_RFPERU_006_E04.x01.phd.1 ZZZ|SDSU_RFPERU_010_B05.x01.phd.1);
 202 is_deeply( [sort map($_->id, $contigs[0]->each_seq)], \@read_ids );
 203 @read_ids = sort qw(sdsu|SDSU_RFPERU_013_H05.x01.phd.1 ABC|SDSU_RFPERU_005_F02.x01.phd.1);
 204 is_deeply( [sort map($_->id, $contigs[1]->each_seq)], \@read_ids );
 205 @read_ids = sort qw( ZZZ|9962187 ABC|9937790 ABC|9944760 ABC|9956706
 206         sdsu|9986984 ABC|9960711 ABC|9970175 ABC|9976538 ABC|9980040);
 207 is_deeply( [sort map($_->id, $contigs[2]->each_seq)], \@read_ids );
 208
 209 # effective assembly params
 210 ok($cross_csp = Bio::Assembly::Tools::ContigSpectrum->new(
 211   -cross          => $mixed_csp,
 212   -eff_asm_params => 1 ), 'cross-contig spectrum');
 213 is_deeply($cross_csp->spectrum, {1=>7, 2=>2, 9=>1}); # [7 2 0 0 0 0 0 0 1]
 214 is($cross_csp->nof_rep, 1);
 215 is($cross_csp->eff_asm_params, 1);
 216 is($cross_csp->max_size, 9);
 217 is($cross_csp->nof_seq, 13);
 218 float_is($cross_csp->avg_seq_len, 206.308);
 219 is($cross_csp->nof_overlaps, 10);
 220 float_is($cross_csp->avg_overlap, 76.9);
 221 float_is($cross_csp->avg_identity, 99.2357);
 222 ## min_overlap and min_identity not explicitely specified for the cross csp
 223 ## min_overlap and min_identity are thus taken from the mixed csp
 224 is($cross_csp->min_overlap, 35);
 225 float_is($cross_csp->min_identity, 96.8421);
 226
 227 # with a specified minimum overlap and identity
 228 ok($cross_csp = Bio::Assembly::Tools::ContigSpectrum->new(
 229   -cross          => $mixed_csp,
 230   -min_overlap    => 50,
 231   -min_identity   => 98 ), 'cross-contig spectrum');
 232 is_deeply($cross_csp->spectrum, {1=>3, 2=>1, 7=>1}); # [3 1 0 0 0 0 1]
 233 is($cross_csp->nof_rep, 1);
 234 is($cross_csp->eff_asm_params, 0);
 235 is($cross_csp->max_size, 7);
 236 is($cross_csp->nof_seq, 9);
 237 float_is($cross_csp->avg_seq_len, 191.222);
 238 is($cross_csp->min_overlap, 50);
 239 float_is($cross_csp->min_identity, 98);
 240
 241 # sum of contig spectra
 242 ok(my $sum_csp = Bio::Assembly::Tools::ContigSpectrum->new(-eff_asm_params=>1), 'contig spectrum sum');
 243 ok($sum_csp->add($dissolved_csp));
 244 ok($sum_csp->add($mixed_csp));
 245 is_deeply($sum_csp->spectrum, {1=>2, 2=>3, 6=>2, 9=>1}); # [2 3 0 0 0 2 0 0 1]
 246 is($sum_csp->eff_asm_params, 1);
 247 is($sum_csp->max_size, 9);
 248 is($sum_csp->nof_rep, 2);
 249 is($sum_csp->nof_seq, 29);
 250 float_is($sum_csp->avg_seq_len, 258.7934);
 251 is($sum_csp->nof_overlaps, 21);
 252 is($sum_csp->min_overlap, 35);
 253 float_is($sum_csp->avg_overlap, 137.0476);
 254 float_is($sum_csp->min_identity, 96.8421);
 255 float_is($sum_csp->avg_identity, 99.1487);
 256 is(scalar $sum_csp->assembly, 4);
 257
 258 # average of contig spectra
 259 ok(my $avg_csp = Bio::Assembly::Tools::ContigSpectrum->new(-eff_asm_params=>1), 'average contig spectrum');
 260 ok($avg_csp = $avg_csp->average([$dissolved_csp, $mixed_csp]));
 261 is_deeply($avg_csp->spectrum, {1=>1, 2=>1.5, 6=>1, 9=>0.5}); # [1 1 0 0 0 1 0 0 0.5]
 262 is($avg_csp->eff_asm_params, 1);
 263 is($avg_csp->max_size, 9);
 264 is($avg_csp->nof_rep, 2);
 265 is($avg_csp->nof_seq, 14.5);
 266 float_is($avg_csp->avg_seq_len, 258.7934);
 267 is($avg_csp->nof_overlaps, 10.5);
 268 is($avg_csp->min_overlap, 35);
 269 float_is($avg_csp->avg_overlap, 137.0476);
 270 float_is($avg_csp->min_identity, 96.8421);
 271 float_is($avg_csp->avg_identity, 99.1487);
 272 is(scalar $avg_csp->assembly, 4);
 273
 274 # drop assembly info from contig spectrum
 275 ok($mixed_csp->drop_assembly(), 'drop assembly');
 276 is(scalar $mixed_csp->assembly(), 0);
 277
 278 # score
 279 my $test_csp;
 280 my $spectrum;
 281 ok($test_csp = Bio::Assembly::Tools::ContigSpectrum->new(-spectrum=>$spectrum), 'contig spectrum score');
 282 is($test_csp->score, undef);
 283 $spectrum = {1=>120};
 284 ok($test_csp = Bio::Assembly::Tools::ContigSpectrum->new(-spectrum=>$spectrum));
 285 is($test_csp->score, 0);
 286 $spectrum = {120=>1};
 287 ok($test_csp = Bio::Assembly::Tools::ContigSpectrum->new(-spectrum=>$spectrum));
 288 is($test_csp->score, 1);
 289 is($test_csp->score(240), 0.248953974895397);
 290 $spectrum = {1=>120, 120=>1};
 291 ok($test_csp = Bio::Assembly::Tools::ContigSpectrum->new(-spectrum=>$spectrum));
 292 is($test_csp->score, 0.248953974895397);
 293
 294 # large contig (27 reads)
 295 $in = Bio::Assembly::IO->new(
 296   -file   => test_input_file('27-contig_Newbler.ace'),
 297   -format => 'ace-454'
 298 );
 299 isa_ok($in, 'Bio::Assembly::IO');
 300 $sc = $in->next_assembly;
 301 isa_ok($sc, 'Bio::Assembly::Scaffold');
 302 ok(my $large_csp = Bio::Assembly::Tools::ContigSpectrum->new(
 303   -assembly       => $sc,
 304   -eff_asm_params => 1 ), 'large contig spectrum');
 305 is(scalar $large_csp->assembly(), 1);
 306 is_deeply($large_csp->spectrum, {1=>0, 27=>1});
 307 is($large_csp->eff_asm_params, 1);
 308 is($large_csp->max_size, 27);
 309 is($large_csp->nof_rep, 1);
 310 is($large_csp->nof_seq, 27);
 311 float_is($large_csp->avg_seq_len, 100);
 312 is($large_csp->nof_overlaps, 26);
 313 is($large_csp->min_overlap, 54);
 314 ok( $large_csp->avg_overlap >= 88.7692307692307 );
 315 float_is($large_csp->min_identity, 33.3333);
 316 ok( $large_csp->avg_identity >= 74 );
 317
 318 ok(my $large_xcsp = Bio::Assembly::Tools::ContigSpectrum->new(
 319   -cross          => $large_csp,
 320   -eff_asm_params => 1           ), 'large cross-contig spectrum');
 321 is($large_xcsp->nof_overlaps, 26);
 322 # operation returns sometimes 88.7692307692308 and sometimes 88.8076923076923...
 323 ok( $large_xcsp->avg_overlap >= 88.7692307692307 );
 324 ok( $large_xcsp->avg_overlap <= 88.8076923076924 );
 325 is_deeply($large_xcsp->spectrum, {1=>21, 27=>1});
 326
 327 ok( $large_xcsp = Bio::Assembly::Tools::ContigSpectrum->new(
 328   -cross          => $large_csp,
 329   -min_overlap    => 100) );
 330 is_deeply($large_xcsp->spectrum, {1=>18, 2=>5, 3=>1, 7=>1});
 331 my @xcontigs = sort {$a->id cmp $b->id} $large_xcsp->assembly;
 332 is(scalar @xcontigs, 7); # the cross-1-contigs are not included
 333 my @xcontig_ids = sort qw( contig00001_1 contig00001_2 contig00001_3 contig00001_4
 334 contig00001_5 contig00001_6 contig00001_7 );
 335 is_deeply( [map($_->id, @xcontigs)], \@xcontig_ids );
 336 my @xcontig_sizes = sort qw( 2 2 2 2 2 3 7 );
 337 is_deeply( [sort map($_->num_sequences, @xcontigs)], \@xcontig_sizes );
 338 my $xcontig = $xcontigs[5];
 339
 340 if ($verbose) {
 341     for my $ctg (@xcontigs) {
 342         print STDERR $ctg->id.":".join(',',$ctg->get_seq_ids)."\n";
 343     }
 344 }
 345
 346 TODO: {
 347     local $TODO= 'hash randomization problems lead to random switching of the assemblies';
 348     lives_ok { $xcontig->get_seq_coord($xcontig->get_seq_by_name('species1635|5973'))->start };
 349     lives_ok { $xcontig->get_seq_coord($xcontig->get_seq_by_name('species158|7890'))->start };
 350     lives_ok { $xcontig->get_seq_coord($xcontig->get_seq_by_name('species2742|48'))->end };
 351
 352     # ORIGINAL TESTS FOLLOWING:
 353     #is( $xcontig->get_seq_coord($xcontig->get_seq_by_name('species1635|5973'))->start, 1);
 354     #is( $xcontig->get_seq_coord($xcontig->get_seq_by_name('species158|7890'))->start, 1);
 355     #is( $xcontig->get_seq_coord($xcontig->get_seq_by_name('species2742|48'))->end, 140);
 356 }
 357
 358 # one contig at a time
 359 $in = Bio::Assembly::IO->new(
 360   -file => test_input_file('contigspectrumtest.tigr'),
 361   -format => 'tigr'
 362 );
 363 $sc = $in->next_assembly;
 364 ok($csp = Bio::Assembly::Tools::ContigSpectrum->new(
 365   -eff_asm_params => 1 ), 'one contig at a time');
 366 for my $contig ($sc->all_contigs) {
 367   ok($csp->assembly($contig));
 368 }
 369
 370 is(scalar $csp->assembly(), 5);
 371 is_deeply($csp->spectrum, {1=>0, 2=>3, 6=>1, 9=>1}); # [0 3 0 0 0 1 0 0 1]
 372 is($csp->eff_asm_params, 1);
 373 is($csp->max_size, 9);
 374 is($csp->nof_rep, 5);
 375 is($csp->nof_seq, 21);
 376 float_is($csp->avg_seq_len, 303.81);
 377 is($csp->nof_overlaps, 16);
 378 is($csp->min_overlap, 35);
 379 is($csp->avg_overlap, 155.875);
 380 float_is($csp->min_identity, 96.8421);
 381 float_is($csp->avg_identity, 98.8826);