t/Maq.t

   1 #-*-perl-*-
   2 #$Id$
   3
   4 use strict;
   5 use warnings;
   6 no warnings qw(once);
   7 our $home;
   8 BEGIN {
   9     $home = '.'; # set to '.' for Build use,
  10                       # '..' for debugging from .t file
  11     unshift @INC, $home;
  12     use Bio::Root::Test;
  13     test_begin(-tests => 52,
  14                -requires_modules => [qw(IPC::Run Bio::Tools::Run::Maq)]);
  15 }
  16
  17 use File::Temp qw(tempfile tempdir);
  18 use Bio::Tools::Run::WrapperBase;
  19
  20 # test command functionality
  21
  22 ok my $maqfac = Bio::Tools::Run::Maq->new(
  23     -command            => 'assemble',
  24     -single_end_quality => 1,
  25     -het_fraction       => 0.005,
  26     -max_mismatches     => 4
  27     ), "make a factory using command 'assemble'";
  28 # ParameterBaseI compliance : really AssemblerBase tests...
  29 ok $maqfac->parameters_changed, "parameters changed on construction";
  30 ok $maqfac->het_fraction, "access parameter";
  31 ok !$maqfac->parameters_changed, "parameters_changed cleared on read";
  32 ok $maqfac->set_parameters( -error_dep_coeff => 0.5 ), "set a param not set in constructor";
  33 ok $maqfac->parameters_changed, "parameters_changed set";
  34 is ($maqfac->error_dep_coeff, 0.5, "parameter really set");
  35 is ($maqfac->het_fraction, 0.005, "original parameter unchanged");
  36 ok !$maqfac->parameters_changed, "parameters_changed cleared on read";
  37 ok $maqfac->set_parameters( -het_fraction => 0.01 ), "change an original parameter";
  38 is ($maqfac->het_fraction, 0.01, "parameter really changed");
  39 ok $maqfac->reset_parameters( -het_fraction => 0.05 ), "reset parameters with arg";
  40 ok !$maqfac->max_mismatches, "original parameters undefined";
  41 is ($maqfac->het_fraction, 0.05, "parameter really reset via arg");
  42 #back to beginning
  43 $maqfac->set_parameters(
  44     -command            => 'assemble',
  45     -single_end_quality => 1,
  46     -het_fraction       => 0.005,
  47     -max_mismatches     => 4
  48     );
  49 ok $maqfac->parameters_changed, "parameters changed";
  50
  51 is( scalar $maqfac->available_parameters, 9, "all available options");
  52 is( scalar $maqfac->available_parameters('params'), 7, "available parameters" );
  53 is( scalar $maqfac->available_parameters('switches'), 2, "available switches" );
  54 my %pms = $maqfac->get_parameters;
  55 is_deeply( \%pms,
  56            { command            => 'assemble',
  57              het_fraction       => 0.005,
  58              max_mismatches     => 4,
  59              single_end_quality => 1}, "get_parameters correct");
  60 is( $maqfac->command, 'assemble', "command attribute set");
  61
  62 is_deeply( $maqfac->{_options}->{_commands},
  63            [@Bio::Tools::Run::Maq::program_commands],
  64            "internal command array set" );
  65
  66 is_deeply( $maqfac->{_options}->{_prefixes},
  67            {%Bio::Tools::Run::Maq::command_prefixes},
  68            "internal prefix hash set");
  69
  70 is_deeply( $maqfac->{_options}->{_params},
  71            [qw( command error_dep_coeff het_fraction max_mismatches max_quality_sum min_map_quality num_haplotypes)],
  72            "commands filtered by prefix");
  73
  74
  75 my @a = @{$maqfac->_translate_params};
  76 is shift @a, 'assemble', 'translate_params: command correct';
  77 my ($k, %h);
  78 for (@a) {
  79     (/^-/) ? ( $h{$k = $_} = undef ) : ( $h{$k} = $_ );
  80 }
  81 is_deeply( \%h, { '-m' => 4, '-r' => 0.005, '-s' => undef }, 'translate_params: options correct');
  82
  83 # test run_maq filearg parsing
  84 # a pipeline...
  85
  86 SKIP : {
  87     test_skip( -requires_executable => $maqfac,
  88                -tests => 27 );
  89     my $rd1 = test_input_file('r1.fq');
  90     my $rd2 = test_input_file('r2.fq');
  91     my $refseq = test_input_file('campycoli.fas');
  92
  93     my $tdir = tempdir( "maqXXXX", CLEANUP => 1);
  94     my ($r1h, $r1f) = tempfile( "rd1XXXX", DIR => $tdir );
  95     $r1h->close;
  96     my ($r2h, $r2f) = tempfile( "rd2XXXX", DIR => $tdir );
  97     $r2h->close;
  98     my ($refh, $reff) = tempfile( "refXXXX", DIR => $tdir );
  99     $refh->close;
 100     my ($map1h, $map1f) = tempfile( "mapXXXX", DIR => $tdir );
 101     $map1h->close;
 102     my ($map2h, $map2f) = tempfile( "mapXXXX", DIR => $tdir );
 103     $map2h->close;
 104     my ($mmaph, $mmapf) = tempfile( "mapXXXX", DIR => $tdir );
 105     $mmaph->close;
 106     my ($cnsh, $cnsf) = tempfile( "cnsXXXX", DIR => $tdir );
 107     $cnsh->close;
 108     my ($maqh, $maqf) = tempfile( "maqXXXX", DIR => $tdir );
 109     $maqh->close;
 110     my ($fqh, $fqf) = tempfile( "faqXXXX", DIR => $tdir );
 111     $fqh->close;
 112
 113     ok $maqfac = Bio::Tools::Run::Maq->new(
 114         -command            => 'fasta2bfa'
 115         ), "make fasta2bfa conversion factory";
 116
 117     ok $maqfac->run_maq( -fas => $refseq,
 118                          -bfa => $reff ), "convert refseq to bfa";
 119
 120     like($maqfac->stderr, qr/1 sequence/, "maq success");
 121
 122     ok $maqfac = Bio::Tools::Run::Maq->new(
 123         -command            => 'fastq2bfq'
 124         ), "make fastq2bfq conversion factory";
 125
 126     ok $maqfac->run_maq( -faq => $rd1,
 127                          -bfq => $r1f ), "convert r1.fq to bfa";
 128     like($maqfac->stderr, qr/125 sequences were loaded/, "maq success");
 129     ok $maqfac->run_maq( -faq => $rd2,
 130                          -bfq => $r2f ), "convert r2.fq to bfa";
 131     like($maqfac->stderr, qr/125 sequences were loaded/, "maq success");
 132
 133     ok $maqfac = Bio::Tools::Run::Maq->new(
 134         -command            => 'map',
 135         ), "make map factory";
 136
 137     ok $maqfac->run_maq( -map => $map1f,
 138                          -bfa => $reff,
 139                          -bfq1 => $r1f ), "map single-end reads";
 140     ok $maqfac->run_maq( -map => $map2f,
 141                          -bfa => $reff,
 142                          -bfq2 => $r2f,
 143                          -bfq1 => $r1f ), "map paired-end reads";
 144
 145     ok $maqfac = Bio::Tools::Run::Maq->new(
 146         -command            => 'mapmerge'
 147         ), "make mapmerge factory";
 148
 149     ok $maqfac->run_maq( -out_map => $mmapf,
 150                          -in_map => [$map1f, $map2f] ), "merge maps";
 151
 152     ok $maqfac = Bio::Tools::Run::Maq->new(
 153         -command            => 'assemble'
 154         ), "make assemble factory";
 155
 156     ok $maqfac->run_maq( -cns => $cnsf,
 157                          -bfa => $reff,
 158                          -map => $mmapf ), "assemble consensus";
 159
 160     ok $maqfac = Bio::Tools::Run::Maq->new(
 161         -command            => 'mapview'
 162         ), "make mapview converter";
 163
 164
 165     ok $maqfac->run_maq( -map => $mmapf,
 166                          -txt => $maqf ), "convert mmap";
 167
 168     ok $maqfac = Bio::Tools::Run::Maq->new(
 169         -command            => 'cns2fq'
 170         ), "make consensus->fastq converter";
 171     ok $maqfac->run_maq( -cns => $cnsf,
 172                          -faq => $fqf ), "convert consensus -> fastq";
 173
 174     # test run (assembly pipeline)
 175     # these parms are the maq defaults for the respective programs
 176     ok $maqfac = Bio::Tools::Run::Maq->new(
 177         -map_max_mismatches => 2,
 178         -asm_max_mismatches => 7,
 179         -c2q_min_map_quality => 40
 180         ), "make an assembly factory";
 181
 182     is( $maqfac->command, 'run', "command attribute set");
 183     is( $maqfac->map_max_mismatches, 2, "map param set");
 184     is( $maqfac->asm_max_mismatches, 7, "asm param set");
 185     is( $maqfac->c2q_min_map_quality, 40, "c2q param set");
 186     ok my $assy = $maqfac->run($rd1,$refseq,$rd2), "make full assy";
 187     #some fuzziness in these: sometimes maq gives 41+4, sometimes 42+6.
 188     cmp_ok( $assy->get_nof_contigs, '>=', 37, "number of contigs");
 189     cmp_ok( $assy->get_nof_singlets,'>=',4, "number of singlets");
 190
 191 }
 192 1;
 193
 194 #  sub test_input_file {
 195 #      return "./data/".shift;
 196 #  }