t/SeqIO/mbsout.t

   1 #!/usr/bin/perl
   2 use version;
   3 our $API_VERSION = qv('1.1.3');
   4
   5 use strict;
   6 use File::Path qw(mkpath rmtree);
   7 use Carp;
   8
   9 BEGIN {
  10     use Bio::Root::Test;
  11
  12     test_begin(
  13         -tests               => 74,
  14         -requires_modules    => [q(Bio::SeqIO::mbsout)],
  15         -requires_networking => 0
  16     );
  17
  18     use_ok('Bio::SeqIO::mbsout');
  19 }
  20
  21 # skip tests if the msout.pm module is too old.
  22 cmp_ok( $Bio::SeqIO::mbsout::API_VERSION,
  23     '>=', qv('1.1.3'), "Bio::SeqIO::mbsout is at least api version 1.1.3" );
  24
  25 create_dir("mbsout");
  26 test_file_1( 0, "mbsout/mbsout_infile1" );
  27 test_file_2( 0, "mbsout/mbsout_infile2" );
  28 test_file_3( 0, "mbsout/mbsout_infile3" );
  29
  30 sub create_dir {
  31
  32     my $dir = shift;
  33
  34     $dir = test_input_file($dir);
  35
  36     unless ( -d $dir ) {
  37         mkpath($dir);
  38     }
  39 }
  40
  41 sub test_file_1 {
  42 ##############################################################################
  43 ## Test file 1
  44 ##############################################################################
  45
  46     my $gzip   = shift;
  47     my $infile = shift;
  48     $infile = test_input_file($infile);
  49
  50     #print_file1( $infile, $gzip );
  51
  52     my $file_sequence = $infile;
  53     if ($gzip) {
  54         $file_sequence = "gunzip -c <$file_sequence |";
  55     }
  56     my $mbsout = Bio::SeqIO->new(
  57         -file   => "$file_sequence",
  58         -format => 'mbsout',
  59     );
  60
  61     isa_ok( $mbsout, 'Bio::SeqIO::mbsout' );
  62
  63     my $rh_base_conversion_table = $mbsout->get_base_conversion_table;
  64
  65     isa_ok( $mbsout, 'Bio::SeqIO::mbsout' );
  66
  67     my %attributes = (
  68         RUNS     => 3,
  69         SEGSITES => 7,
  70         MBS_INFO_LINE =>
  71           'command:        mbs 6 -t 0.001 -r 0.00025 -s 5000 2500 -f 3 1 traj ',
  72         TOT_RUN_HAPS      => 6,
  73         NEXT_RUN_NUM      => 1,
  74         LAST_READ_HAP_NUM => 0,
  75         POSITIONS => [qw(79.1001 80.1001 81.101 82.101 83.10001 84.801 85)],
  76         CURRENT_RUN_SEGSITES      => 7,
  77         POP_MUT_PARAM_PER_SITE    => 0.001,
  78         POP_RECOMB_PARAM_PER_SITE => 0.00025,
  79         NSITES                    => 5000,
  80         SELPOS                    => 2500,
  81         NFILES                    => 3,
  82         NREPS                     => 1,
  83         TRAJ_FILENAME             => 'traj'
  84     );
  85
  86     foreach my $attribute ( keys %attributes ) {
  87         my $func = lc($attribute);
  88
  89         if ( $attribute =~ m/POPS|SEEDS|POSITIONS/ ) {
  90             $func = ucfirst($func);
  91         }
  92
  93         $func = 'get_' . $func;
  94         my @returns = $mbsout->$func();
  95         my ( $return, $got );
  96
  97         # If there were more than one return value, then compare references to
  98         # arrays instead of scalars
  99         unless ( @returns > 1 ) {
 100             $got = shift @returns;
 101         }
 102         else { $got = \@returns }
 103
 104         my $expected = $attributes{$attribute};
 105
 106         if ( defined $got && defined $expected ) {
 107             is_deeply( $got, $expected, "Get $attribute" );
 108         }
 109         else { is_deeply( $got, $expected, "Get $attribute" ) }
 110     }
 111
 112     # Testing next_hap at beginning of run
 113     my @data_got =
 114       convert_bases_to_nums( $rh_base_conversion_table, $mbsout->get_next_seq );
 115     my @data_expected = qw(1111111);
 116     is_deeply( \@data_got, \@data_expected,
 117         "Get next_hap at beginning of run" );
 118
 119     # Testing next_hap after beginning of run
 120     @data_got =
 121       convert_bases_to_nums( $rh_base_conversion_table, $mbsout->get_next_seq );
 122     @data_expected = qw(5555555);
 123     is_deeply( \@data_got, \@data_expected,
 124         "Get next_hap after beginning of run" );
 125
 126     # Testing next_pop after beginning of pop
 127     @data_got =
 128       convert_bases_to_nums( $rh_base_conversion_table, $mbsout->get_next_seq );
 129     @data_expected = qw(4444444);
 130     is_deeply( \@data_got, \@data_expected,
 131         "Get next_pop after beginning of pop" );
 132
 133     # Testing next_hap
 134     @data_got      = $mbsout->get_next_hap;
 135     @data_expected = qw(4444444);
 136     is_deeply( \@data_got, \@data_expected, "Get next_hap" );
 137
 138     # Testing next_hap
 139     @data_got      = $mbsout->get_next_hap;
 140     @data_expected = qw(5555555);
 141     is_deeply( \@data_got, \@data_expected, "Get next_hap" );
 142
 143     # Testing next_run after beginning of run
 144     @data_got =
 145       convert_bases_to_nums( $rh_base_conversion_table, $mbsout->get_next_run );
 146     @data_expected = qw(4444444);
 147     is_deeply( \@data_got, \@data_expected,
 148         "Get next_run after beginning of run" );
 149
 150     # Testing next_run at beginning of run
 151     @data_got =
 152       convert_bases_to_nums( $rh_base_conversion_table, $mbsout->get_next_run );
 153     @data_expected = qw(5555555 5555555 5555555 1010101 1111111 1515151);
 154     is_deeply( \@data_got, \@data_expected,
 155         "Get next_run at beginning of run" );
 156
 157     # Testing next_run at beginning of run
 158     @data_got =
 159       convert_bases_to_nums( $rh_base_conversion_table, $mbsout->get_next_run );
 160     @data_expected = qw(
 161       1414141
 162       1414141
 163       1515151
 164       1414141
 165       1515151
 166       1515151);
 167     is_deeply( \@data_got, \@data_expected,
 168         "Get next_run at beginning of run" );
 169
 170     is( $mbsout->get_next_run_num, undef, 'have all lines been read?' );
 171 }
 172
 173 sub test_file_2 {
 174 ##############################################################################
 175 ## Test file 2
 176 ##############################################################################
 177
 178     my $gzip   = shift;
 179     my $infile = shift;
 180     $infile = test_input_file($infile);
 181
 182     #print_file2( $infile, $gzip );
 183
 184     my $file_sequence = $infile;
 185     if ($gzip) {
 186         $file_sequence = "gunzip -c <$file_sequence |";
 187     }
 188
 189     my $mbsout = Bio::SeqIO->new(
 190         -file   => "$file_sequence",
 191         -format => 'mbsout',
 192     );
 193
 194     isa_ok( $mbsout, 'Bio::SeqIO::mbsout' );
 195
 196     my %attributes = (
 197         RUNS     => 5,
 198         SEGSITES => 7,
 199         MBS_INFO_LINE =>
 200           'command:        mbs 6 -t 0.001 -r 0.00025 -s 5000 2500 -f 5 1 traj ',
 201         TOT_RUN_HAPS      => 6,
 202         NEXT_RUN_NUM      => 1,
 203         LAST_READ_HAP_NUM => 0,
 204         POSITIONS => [qw(79.1001 80.1001 81.101 82.101 83.10001 84.801 85)],
 205         CURRENT_RUN_SEGSITES      => 7,
 206         POP_MUT_PARAM_PER_SITE    => 0.001,
 207         POP_RECOMB_PARAM_PER_SITE => 0.00025,
 208         NSITES                    => 5000,
 209         SELPOS                    => 2500,
 210         NFILES                    => 5,
 211         NREPS                     => 1,
 212         TRAJ_FILENAME             => 'traj'
 213     );
 214
 215     foreach my $attribute ( keys %attributes ) {
 216         my $func = lc($attribute);
 217         if ( $attribute =~ m/POSITIONS/ ) {
 218             $func = ucfirst($func);
 219         }
 220         elsif ( $attribute =~ m/\_file/ ) {
 221             $func = q(infile);
 222         }
 223
 224         $func = 'get_' . $func;
 225         my @returns = $mbsout->$func();
 226         my ( $return, $got );
 227
 228         # If there were more than one return value, then compare references to
 229         # arrays instead of scalars
 230         unless ( @returns > 1 ) {
 231             $got = shift @returns;
 232         }
 233         else { $got = \@returns }
 234
 235         my $expected = $attributes{$attribute};
 236
 237         if ( defined $got && defined $expected ) {
 238             is_deeply( $got, $expected, "Get $attribute" );
 239         }
 240         else { is_deeply( $got, $expected, "Get $attribute" ) }
 241     }
 242
 243     my $rh_base_conversion_table = $mbsout->get_base_conversion_table;
 244
 245     # Testing next_hap at beginning of run
 246     my @data_got      = $mbsout->get_next_hap;
 247     my @data_expected = qw(1111111);
 248     is_deeply( \@data_got, \@data_expected,
 249         "Get next_hap at beginning of run" );
 250
 251     # Testing next_hap after beginning of run
 252     @data_got =
 253       convert_bases_to_nums( $rh_base_conversion_table, $mbsout->get_next_seq );
 254     @data_expected = '5555555';
 255     is_deeply( \@data_got, \@data_expected,
 256         "Get next_hap after beginning of run" );
 257
 258     # Surprise test! testing mbsout::outgroup
 259     my $outgroup = $mbsout->outgroup;
 260     is( $outgroup, 0, "Testing mbsout::outgroup" );
 261
 262     # Testing next_run after beginning of run
 263     @data_got =
 264       convert_bases_to_nums( $rh_base_conversion_table, $mbsout->get_next_run );
 265     @data_expected = qw(
 266       4444444
 267       4444444
 268       5555555
 269       4444444);
 270     is_deeply( \@data_got, \@data_expected,
 271         "Get next_run after beginning of run" );
 272
 273     # Testing next_run after beginning of run
 274     @data_got =
 275       convert_bases_to_nums( $rh_base_conversion_table, $mbsout->get_next_run );
 276     @data_expected = qw(
 277       5555555
 278       5555555
 279       5555555
 280       1010101
 281       1111111
 282       1515151);
 283     is_deeply( \@data_got, \@data_expected,
 284         "Get next_run after beginning of run" );
 285
 286     # Testing next_run at beginning of run
 287     @data_got =
 288       convert_bases_to_nums( $rh_base_conversion_table, $mbsout->get_next_run );
 289     @data_expected = qw(
 290       1414141
 291       1414141
 292       1515151
 293       1414141
 294       1515151
 295       1515151);
 296     is_deeply( \@data_got, \@data_expected,
 297         "Get next_run at beginning of run" );
 298
 299     # Testing next_hap at beginning of run 2
 300     @data_got =
 301       convert_bases_to_nums( $rh_base_conversion_table, $mbsout->get_next_seq );
 302     @data_expected = '1515151';
 303     is_deeply( \@data_got, \@data_expected,
 304         "Get next_hap at beginning of run 2" );
 305
 306     # Testing next_run after hap
 307     @data_got =
 308       convert_bases_to_nums( $rh_base_conversion_table, $mbsout->get_next_run );
 309     @data_expected = qw(
 310       5050505
 311       5151515
 312       5555555
 313       5454545
 314       5454545);
 315     is_deeply( \@data_got, \@data_expected, "Get next_run after hap" );
 316
 317     is( $mbsout->get_next_run_num, 5, 'next run should be 5.' );
 318
 319 }
 320
 321 sub test_file_3 {
 322 ##############################################################################
 323 ## Test file 3
 324 ##############################################################################
 325
 326     my $gzip   = shift;
 327     my $infile = shift;
 328     $infile = test_input_file($infile);
 329
 330     #print_file3( $infile, $gzip );
 331
 332     my $file_sequence = $infile;
 333     if ($gzip) {
 334         $file_sequence = "gunzip -c <$file_sequence |";
 335     }
 336     my $mbsout = Bio::SeqIO->new(
 337         -file   => "$file_sequence",
 338         -format => 'mbsout',
 339     );
 340
 341     isa_ok( $mbsout, 'Bio::SeqIO::mbsout' );
 342
 343     my $rh_base_conversion_table = $mbsout->get_base_conversion_table;
 344
 345     my %attributes = (
 346         RUNS     => 1,
 347         SEGSITES => 7,
 348         MBS_INFO_LINE =>
 349           'command:        mbs 3 -t 0.001 -r 0.00025 -s 5000 2500 -f 1 1 traj ',
 350         TOT_RUN_HAPS      => 3,
 351         NEXT_RUN_NUM      => 1,
 352         LAST_READ_HAP_NUM => 0,
 353         POSITIONS => [qw(79.1001 80.1001 81.101 82.101 83.10001 84.801 85)],
 354         CURRENT_RUN_SEGSITES      => 7,
 355         POP_MUT_PARAM_PER_SITE    => 0.001,
 356         POP_RECOMB_PARAM_PER_SITE => 0.00025,
 357         NSITES                    => 5000,
 358         SELPOS                    => 2500,
 359         NFILES                    => 1,
 360         NREPS                     => 1,
 361         TRAJ_FILENAME             => 'traj'
 362     );
 363
 364     foreach my $attribute ( keys %attributes ) {
 365         my $func = lc($attribute);
 366
 367         if ( $attribute =~ m/POPS|SEEDS|POSITIONS/ ) {
 368             $func = ucfirst($func);
 369         }
 370
 371         $func = 'get_' . $func;
 372         my @returns = $mbsout->$func();
 373         my ( $return, $got );
 374
 375         # If there were more than one return value, then compare references to
 376         # arrays instead of scalars
 377         unless ( @returns > 1 ) {
 378             $got = shift @returns;
 379         }
 380         else { $got = \@returns }
 381
 382         my $expected = $attributes{$attribute};
 383
 384         if ( defined $got && defined $expected ) {
 385             is_deeply( $got, $expected, "Get $attribute" );
 386         }
 387         else { is_deeply( $got, $expected, "Get $attribute" ) }
 388     }
 389
 390     # Testing next_run at beginning of run
 391     my @data_got =
 392       convert_bases_to_nums( $rh_base_conversion_table, $mbsout->get_next_run );
 393     my @data_expected = qw(1111111 5555555 4444444);
 394     is_deeply( \@data_got, \@data_expected,
 395         "Get next_run at end/beginning of run" );
 396
 397     is( $mbsout->get_next_run_num, undef, 'have all lines been read?' );
 398
 399     # Testing what happens when we read from empty stream
 400     @data_got      = $mbsout->get_next_run;
 401     @data_expected = ();
 402     is_deeply( \@data_got, \@data_expected, "Get next_run at eof" );
 403
 404     # Testing what happens when we read from empty stream
 405     @data_got      = $mbsout->get_next_hap;
 406     @data_expected = undef;
 407     is_deeply( \@data_got, \@data_expected, "Get next_hap at eof" );
 408
 409     # Testing what happens when we read from empty stream
 410     @data_got      = $mbsout->get_next_seq;
 411     @data_expected = ();
 412     is_deeply( \@data_got, \@data_expected, "Get next_seq at eof" );
 413
 414 }
 415
 416 sub print_file1 {
 417
 418     my $destination = shift;
 419     my $gzip        = shift;
 420
 421     my $out = <<END
 422 command:        mbs 6 -t 0.001 -r 0.00025 -s 5000 2500 -f 3 1 traj
 423
 424 //0-1   allele: a a a d d d d
 425 segsites: 7
 426 positions: 79.1001 80.1001 81.101 82.101 83.10001 84.801 85
 427 1111111
 428 5555555
 429 4444444
 430 4444444
 431 5555555
 432 4444444
 433 //1-1   allele: a a a a d d d
 434 segsites: 7
 435 positions: 79.1001 80.1001 81.101 82.101 83.10001 84.801 85
 436 5555555
 437 5555555
 438 5555555
 439 1010101
 440 1111111
 441 1515151
 442 //2-1   allele: a d d d d d d
 443 segsites: 7
 444 positions: 79.1001 80.1001 81.101 82.101 83.10001 84.801 85
 445 1414141
 446 1414141
 447 1515151
 448 1414141
 449 1515151
 450 1515151
 451 END
 452       ;
 453
 454     if ($gzip) {
 455         $gzip = "| gzip";
 456     }
 457     else { $gzip = ' '; }
 458     open my $OUT, "$gzip >$destination" or croak "Could not write file '$destination': $!\n";
 459
 460     print $OUT $out;
 461     close $OUT;
 462 }
 463
 464 sub print_file2 {
 465
 466     my $destination = shift;
 467     my $gzip        = shift;
 468
 469     my $out = <<END
 470 command:        mbs 6 -t 0.001 -r 0.00025 -s 5000 2500 -f 5 1 traj
 471
 472 //0-1   allele: a a a d d d d
 473 segsites: 7
 474 positions: 79.1001 80.1001 81.101 82.101 83.10001 84.801 85
 475 1111111
 476 5555555
 477 4444444
 478 4444444
 479 5555555
 480 4444444
 481 //1-1   allele: a a a d d d d
 482 segsites: 7
 483 positions: 79.1001 80.1001 81.101 82.101 83.10001 84.801 85
 484 5555555
 485 5555555
 486 5555555
 487 1010101
 488 1111111
 489 1515151
 490 //2-1   allele: a a a d d d d
 491
 492 segsites: 7
 493 positions: 79.1001 80.1001 81.101 82.101 83.10001 84.801 85
 494 1414141
 495 1414141
 496 1515151
 497 1414141
 498 1515151
 499 1515151
 500 //3-1   allele: a a a a a d d
 501
 502 segsites: 7
 503 positions: 79.1001 80.1001 81.101 82.101 83.10001 84.801 85
 504 1515151
 505 5050505
 506 5151515
 507 5555555
 508 5454545
 509 5454545
 510 //4-1   allele: a a d d d d d
 511
 512 segsites: 7
 513 positions: 79.1001 80.1001 81.101 82.101 83.10001 84.801 85
 514 5555555
 515 END
 516       ;
 517
 518     if ($gzip) {
 519         $gzip = "| gzip";
 520     }
 521     else { $gzip = ' '; }
 522
 523     open my $OUT, "$gzip >$destination" or croak "Could not write file '$destination': $!\n";
 524
 525     print $OUT $out;
 526     close $OUT;
 527 }
 528
 529 sub print_file3 {
 530
 531     my $destination = shift;
 532     my $gzip        = shift;
 533
 534     my $out = <<END ;
 535 command:        mbs 3 -t 0.001 -r 0.00025 -s 5000 2500 -f 1 1 traj
 536
 537 //0-1   allele: a a a d d d d
 538 segsites: 7
 539 positions: 79.1001 80.1001 81.101 82.101 83.10001 84.801 85
 540 1111111
 541 5555555
 542 4444444
 543 END
 544
 545     if ($gzip) {
 546         $gzip = "| gzip";
 547     }
 548     else { $gzip = ' '; }
 549
 550     open my $OUT, "$gzip >$destination" or croak "Could not write file '$destination': $!\n";
 551
 552     print $OUT $out;
 553     close $OUT;
 554 }
 555
 556 sub print_to_file {
 557     my ( $ra_in, $out ) = @_;
 558     open my $OUT, '>', $out or croak "\nCould not write outfile '$out': $!\n";
 559     print $OUT ("@$ra_in");
 560     close $OUT;
 561 }
 562
 563 sub convert_bases_to_nums {
 564     my ( $rh_base_conversion_table, @seqs ) = @_;
 565
 566     my @out_seqstrings;
 567     foreach my $seq (@seqs) {
 568         my $seqstring = $seq->seq;
 569         foreach my $base ( keys %{$rh_base_conversion_table} ) {
 570             $seqstring =~ s/($base)/$rh_base_conversion_table->{$base}/g;
 571         }
 572         push @out_seqstrings, $seqstring;
 573     }
 574     return @out_seqstrings;
 575 }