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 => 51,
  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 is( join(' ', @{$maqfac->_translate_params}),
  74     "assemble -m 4 -r 0.005 -s", "translate params" );
  75
  76 # test run_maq filearg parsing
  77 # a pipeline...
  78
  79 SKIP : {
  80     test_skip( -requires_executable => $maqfac,
  81                -tests => 27 );
  82     my $rd1 = test_input_file('r1.fq');
  83     my $rd2 = test_input_file('r2.fq');
  84     my $refseq = test_input_file('campycoli.fas');
  85
  86     my $tdir = tempdir( "maqXXXX", CLEANUP => 1);
  87     my ($r1h, $r1f) = tempfile( "rd1XXXX", DIR => $tdir );
  88     $r1h->close;
  89     my ($r2h, $r2f) = tempfile( "rd2XXXX", DIR => $tdir );
  90     $r2h->close;
  91     my ($refh, $reff) = tempfile( "refXXXX", DIR => $tdir );
  92     $refh->close;
  93     my ($map1h, $map1f) = tempfile( "mapXXXX", DIR => $tdir );
  94     $map1h->close;
  95     my ($map2h, $map2f) = tempfile( "mapXXXX", DIR => $tdir );
  96     $map2h->close;
  97     my ($mmaph, $mmapf) = tempfile( "mapXXXX", DIR => $tdir );
  98     $mmaph->close;
  99     my ($cnsh, $cnsf) = tempfile( "cnsXXXX", DIR => $tdir );
 100     $cnsh->close;
 101     my ($maqh, $maqf) = tempfile( "maqXXXX", DIR => $tdir );
 102     $maqh->close;
 103     my ($fqh, $fqf) = tempfile( "faqXXXX", DIR => $tdir );
 104     $fqh->close;
 105
 106     ok $maqfac = Bio::Tools::Run::Maq->new(
 107         -command            => 'fasta2bfa'
 108         ), "make fasta2bfa conversion factory";
 109
 110     ok $maqfac->run_maq( -fas => $refseq,
 111                          -bfa => $reff ), "convert refseq to bfa";
 112
 113     like($maqfac->stderr, qr/1 sequence/, "maq success");
 114
 115     ok $maqfac = Bio::Tools::Run::Maq->new(
 116         -command            => 'fastq2bfq'
 117         ), "make fastq2bfq conversion factory";
 118
 119     ok $maqfac->run_maq( -faq => $rd1,
 120                          -bfq => $r1f ), "convert r1.fq to bfa";
 121     like($maqfac->stderr, qr/125 sequences were loaded/, "maq success");
 122     ok $maqfac->run_maq( -faq => $rd2,
 123                          -bfq => $r2f ), "convert r2.fq to bfa";
 124     like($maqfac->stderr, qr/125 sequences were loaded/, "maq success");
 125
 126     ok $maqfac = Bio::Tools::Run::Maq->new(
 127         -command            => 'map',
 128         ), "make map factory";
 129
 130     ok $maqfac->run_maq( -map => $map1f,
 131                          -bfa => $reff,
 132                          -bfq1 => $r1f ), "map single-end reads";
 133     ok $maqfac->run_maq( -map => $map2f,
 134                          -bfa => $reff,
 135                          -bfq2 => $r2f,
 136                          -bfq1 => $r1f ), "map paired-end reads";
 137
 138     ok $maqfac = Bio::Tools::Run::Maq->new(
 139         -command            => 'mapmerge'
 140         ), "make mapmerge factory";
 141
 142     ok $maqfac->run_maq( -out_map => $mmapf,
 143                          -in_map => [$map1f, $map2f] ), "merge maps";
 144
 145     ok $maqfac = Bio::Tools::Run::Maq->new(
 146         -command            => 'assemble'
 147         ), "make assemble factory";
 148
 149     ok $maqfac->run_maq( -cns => $cnsf,
 150                          -bfa => $reff,
 151                          -map => $mmapf ), "assemble consensus";
 152
 153     ok $maqfac = Bio::Tools::Run::Maq->new(
 154         -command            => 'mapview'
 155         ), "make mapview converter";
 156
 157
 158     ok $maqfac->run_maq( -map => $mmapf,
 159                          -txt => $maqf ), "convert mmap";
 160
 161     ok $maqfac = Bio::Tools::Run::Maq->new(
 162         -command            => 'cns2fq'
 163         ), "make consensus->fastq converter";
 164     ok $maqfac->run_maq( -cns => $cnsf,
 165                          -faq => $fqf ), "convert consensus -> fastq";
 166
 167     # test run (assembly pipeline)
 168     # these parms are the maq defaults for the respective programs
 169     ok $maqfac = Bio::Tools::Run::Maq->new(
 170         -map_max_mismatches => 2,
 171         -asm_max_mismatches => 7,
 172         -c2q_min_map_quality => 40
 173         ), "make an assembly factory";
 174
 175     is( $maqfac->command, 'run', "command attribute set");
 176     is( $maqfac->map_max_mismatches, 2, "map param set");
 177     is( $maqfac->asm_max_mismatches, 7, "asm param set");
 178     is( $maqfac->c2q_min_map_quality, 40, "c2q param set");
 179     ok my $assy = $maqfac->run($rd1,$refseq,$rd2), "make full assy";
 180     #some fuzziness in these: sometimes maq gives 41+4, sometimes 42+6.
 181     cmp_ok( $assy->get_nof_contigs, '>=', 37, "number of contigs");
 182     cmp_ok( $assy->get_nof_singlets,'>=',4, "number of singlets");
 183
 184 }
 185 1;
 186
 187 #  sub test_input_file {
 188 #      return "./data/".shift;
 189 #  }