t/Assembly/ContigSpectrum.t

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