3 our $API_VERSION = qv
('1.1.3');
6 use File
::Path
qw(mkpath rmtree);
10 use lib
'.'; # for core package test scripts only
15 -requires_modules
=> [q
(Bio
::SeqIO
::mbsout
)],
16 -requires_networking
=> 0
19 use_ok
('Bio::SeqIO::mbsout');
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" );
27 test_file_1
( 0, "mbsout/mbsout_infile1" );
28 test_file_2
( 0, "mbsout/mbsout_infile2" );
29 test_file_3
( 0, "mbsout/mbsout_infile3" );
35 $dir = test_input_file
($dir);
43 ##############################################################################
45 ##############################################################################
49 $infile = test_input_file
($infile);
51 #print_file1( $infile, $gzip );
53 my $file_sequence = $infile;
55 $file_sequence = "gunzip -c <$file_sequence |";
57 my $mbsout = Bio
::SeqIO
->new(
58 -file
=> "$file_sequence",
62 isa_ok
( $mbsout, 'Bio::SeqIO::mbsout' );
64 my $rh_base_conversion_table = $mbsout->get_base_conversion_table;
66 isa_ok
( $mbsout, 'Bio::SeqIO::mbsout' );
72 'command: mbs 6 -t 0.001 -r 0.00025 -s 5000 2500 -f 3 1 traj ',
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,
84 TRAJ_FILENAME
=> 'traj'
87 foreach my $attribute ( keys %attributes ) {
88 my $func = lc($attribute);
90 if ( $attribute =~ m/POPS|SEEDS|POSITIONS/ ) {
91 $func = ucfirst($func);
94 $func = 'get_' . $func;
95 my @returns = $mbsout->$func();
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;
103 else { $got = \
@returns }
105 my $expected = $attributes{$attribute};
107 if ( defined $got && defined $expected ) {
108 is_deeply
( $got, $expected, "Get $attribute" );
110 else { is_deeply
( $got, $expected, "Get $attribute" ) }
113 # Testing next_hap at beginning of run
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" );
120 # Testing next_hap after beginning of run
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" );
127 # Testing next_pop after beginning of pop
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" );
135 @data_got = $mbsout->get_next_hap;
136 @data_expected = qw(4444444);
137 is_deeply
( \
@data_got, \
@data_expected, "Get next_hap" );
140 @data_got = $mbsout->get_next_hap;
141 @data_expected = qw(5555555);
142 is_deeply
( \
@data_got, \
@data_expected, "Get next_hap" );
144 # Testing next_run after beginning of run
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" );
151 # Testing next_run at beginning of run
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" );
158 # Testing next_run at beginning of run
160 convert_bases_to_nums
( $rh_base_conversion_table, $mbsout->get_next_run );
168 is_deeply
( \
@data_got, \
@data_expected,
169 "Get next_run at beginning of run" );
171 is
( $mbsout->get_next_run_num, undef, 'have all lines been read?' );
175 ##############################################################################
177 ##############################################################################
181 $infile = test_input_file
($infile);
183 #print_file2( $infile, $gzip );
185 my $file_sequence = $infile;
187 $file_sequence = "gunzip -c <$file_sequence |";
190 my $mbsout = Bio
::SeqIO
->new(
191 -file
=> "$file_sequence",
195 isa_ok
( $mbsout, 'Bio::SeqIO::mbsout' );
201 'command: mbs 6 -t 0.001 -r 0.00025 -s 5000 2500 -f 5 1 traj ',
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,
213 TRAJ_FILENAME
=> 'traj'
216 foreach my $attribute ( keys %attributes ) {
217 my $func = lc($attribute);
218 if ( $attribute =~ m/POSITIONS/ ) {
219 $func = ucfirst($func);
221 elsif ( $attribute =~ m/\_file/ ) {
225 $func = 'get_' . $func;
226 my @returns = $mbsout->$func();
227 my ( $return, $got );
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;
234 else { $got = \
@returns }
236 my $expected = $attributes{$attribute};
238 if ( defined $got && defined $expected ) {
239 is_deeply
( $got, $expected, "Get $attribute" );
241 else { is_deeply
( $got, $expected, "Get $attribute" ) }
244 my $rh_base_conversion_table = $mbsout->get_base_conversion_table;
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" );
252 # Testing next_hap after beginning of run
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" );
259 # Surprise test! testing mbsout::outgroup
260 my $outgroup = $mbsout->outgroup;
261 is
( $outgroup, 0, "Testing mbsout::outgroup" );
263 # Testing next_run after beginning of run
265 convert_bases_to_nums
( $rh_base_conversion_table, $mbsout->get_next_run );
271 is_deeply
( \
@data_got, \
@data_expected,
272 "Get next_run after beginning of run" );
274 # Testing next_run after beginning of run
276 convert_bases_to_nums
( $rh_base_conversion_table, $mbsout->get_next_run );
284 is_deeply
( \
@data_got, \
@data_expected,
285 "Get next_run after beginning of run" );
287 # Testing next_run at beginning of run
289 convert_bases_to_nums
( $rh_base_conversion_table, $mbsout->get_next_run );
297 is_deeply
( \
@data_got, \
@data_expected,
298 "Get next_run at beginning of run" );
300 # Testing next_hap at beginning of run 2
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" );
307 # Testing next_run after hap
309 convert_bases_to_nums
( $rh_base_conversion_table, $mbsout->get_next_run );
316 is_deeply
( \
@data_got, \
@data_expected, "Get next_run after hap" );
318 is
( $mbsout->get_next_run_num, 5, 'next run should be 5.' );
323 ##############################################################################
325 ##############################################################################
329 $infile = test_input_file
($infile);
331 #print_file3( $infile, $gzip );
333 my $file_sequence = $infile;
335 $file_sequence = "gunzip -c <$file_sequence |";
337 my $mbsout = Bio
::SeqIO
->new(
338 -file
=> "$file_sequence",
342 isa_ok
( $mbsout, 'Bio::SeqIO::mbsout' );
344 my $rh_base_conversion_table = $mbsout->get_base_conversion_table;
350 'command: mbs 3 -t 0.001 -r 0.00025 -s 5000 2500 -f 1 1 traj ',
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,
362 TRAJ_FILENAME
=> 'traj'
365 foreach my $attribute ( keys %attributes ) {
366 my $func = lc($attribute);
368 if ( $attribute =~ m/POPS|SEEDS|POSITIONS/ ) {
369 $func = ucfirst($func);
372 $func = 'get_' . $func;
373 my @returns = $mbsout->$func();
374 my ( $return, $got );
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;
381 else { $got = \
@returns }
383 my $expected = $attributes{$attribute};
385 if ( defined $got && defined $expected ) {
386 is_deeply
( $got, $expected, "Get $attribute" );
388 else { is_deeply
( $got, $expected, "Get $attribute" ) }
391 # Testing next_run at beginning of run
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" );
398 is
( $mbsout->get_next_run_num, undef, 'have all lines been read?' );
400 # Testing what happens when we read from empty stream
401 @data_got = $mbsout->get_next_run;
403 is_deeply
( \
@data_got, \
@data_expected, "Get next_run at eof" );
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" );
410 # Testing what happens when we read from empty stream
411 @data_got = $mbsout->get_next_seq;
413 is_deeply
( \
@data_got, \
@data_expected, "Get next_seq at eof" );
419 my $destination = shift;
423 command: mbs 6 -t 0.001 -r 0.00025 -s 5000 2500 -f 3 1 traj
425 //0-1 allele: a a a d d d d
427 positions: 79.1001 80.1001 81.101 82.101 83.10001 84.801 85
434 //1-1 allele: a a a a d d d
436 positions: 79.1001 80.1001 81.101 82.101 83.10001 84.801 85
443 //2-1 allele: a d d d d d d
445 positions: 79.1001 80.1001 81.101 82.101 83.10001 84.801 85
458 else { $gzip = ' '; }
459 open my $OUT, "$gzip >$destination" or croak
"Could not write file '$destination': $!\n";
467 my $destination = shift;
471 command: mbs 6 -t 0.001 -r 0.00025 -s 5000 2500 -f 5 1 traj
473 //0-1 allele: a a a d d d d
475 positions: 79.1001 80.1001 81.101 82.101 83.10001 84.801 85
482 //1-1 allele: a a a d d d d
484 positions: 79.1001 80.1001 81.101 82.101 83.10001 84.801 85
491 //2-1 allele: a a a d d d d
494 positions: 79.1001 80.1001 81.101 82.101 83.10001 84.801 85
501 //3-1 allele: a a a a a d d
504 positions: 79.1001 80.1001 81.101 82.101 83.10001 84.801 85
511 //4-1 allele: a a d d d d d
514 positions: 79.1001 80.1001 81.101 82.101 83.10001 84.801 85
522 else { $gzip = ' '; }
524 open my $OUT, "$gzip >$destination" or croak
"Could not write file '$destination': $!\n";
532 my $destination = shift;
536 command: mbs 3 -t 0.001 -r 0.00025 -s 5000 2500 -f 1 1 traj
538 //0-1 allele: a a a d d d d
540 positions: 79.1001 80.1001 81.101 82.101 83.10001 84.801 85
549 else { $gzip = ' '; }
551 open my $OUT, "$gzip >$destination" or croak
"Could not write file '$destination': $!\n";
558 my ( $ra_in, $out ) = @_;
559 open my $OUT, '>', $out or croak
"\nCould not write outfile '$out': $!\n";
560 print $OUT ("@$ra_in");
564 sub convert_bases_to_nums
{
565 my ( $rh_base_conversion_table, @seqs ) = @_;
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;
573 push @out_seqstrings, $seqstring;
575 return @out_seqstrings;