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