t/BWA.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     unshift @INC, '../..';
  13     use Bio::Root::Test;
  14     test_begin(-tests => 37,
  15                -requires_modules => [qw(IPC::Run Bio::Tools::Run::BWA Bio::Tools::Run::Samtools Bio::DB::Sam)]);
  16 }
  17
  18 use File::Copy;
  19
  20 use Bio::Tools::Run::WrapperBase;
  21 use Bio::Assembly::IO::sam;
  22
  23 # test command functionality
  24
  25 ok my $bwafac = Bio::Tools::Run::BWA->new(
  26     -command              => 'aln',
  27     -n_threads            => 1,
  28     -subopt_hit_threshold => 35
  29     ), "make a factory using command 'aln'";
  30
  31 # ParameterBaseI compliance : really AssemblerBase tests...
  32 ok $bwafac->parameters_changed, "parameters changed on construction";
  33 ok $bwafac->subopt_hit_threshold, "access parameter";
  34 ok !$bwafac->parameters_changed, "parameters_changed cleared on read";
  35 ok $bwafac->set_parameters( -reverse_no_comp => 1 ), "set a param not set in constructor";
  36 ok $bwafac->parameters_changed, "parameters_changed set";
  37 ok ($bwafac->reverse_no_comp, "parameter really set");
  38 is ($bwafac->subopt_hit_threshold, 35, "original parameter unchanged");
  39 ok !$bwafac->parameters_changed, "parameters_changed cleared on read";
  40 ok $bwafac->set_parameters( -subopt_hit_threshold => 33 ), "change an original parameter";
  41 is ($bwafac->subopt_hit_threshold, 33, "parameter really changed");
  42 ok $bwafac->reset_parameters( -subopt_hit_threshold => 34 ), "reset parameters with arg";
  43 ok !$bwafac->n_threads, "original parameters undefined";
  44 is ($bwafac->subopt_hit_threshold, 34, "parameter really reset via arg");
  45 #back to beginning
  46 $bwafac->set_parameters(
  47     -command            => 'aln',
  48     -n_threads          => 1,
  49     -subopt_hit_threshold => 35
  50     );
  51 ok $bwafac->parameters_changed, "parameters changed";
  52
  53 is( scalar $bwafac->available_parameters, 16, "all available options");
  54 is( scalar $bwafac->available_parameters('params'), 14, "available parameters" );
  55 is( scalar $bwafac->available_parameters('switches'), 2, "available switches" );
  56 my %pms = $bwafac->get_parameters;
  57 is_deeply( \%pms,
  58            { command            => 'aln',
  59              subopt_hit_threshold => 35,
  60              n_threads     => 1}, "get_parameters correct");
  61 is( $bwafac->command, 'aln', "command attribute set");
  62
  63 is_deeply( $bwafac->{_options}->{_commands},
  64            [@Bio::Tools::Run::BWA::program_commands],
  65            "internal command array set" );
  66
  67 is_deeply( $bwafac->{_options}->{_prefixes},
  68            {%Bio::Tools::Run::BWA::command_prefixes},
  69            "internal prefix hash set");
  70
  71 is_deeply( $bwafac->{_options}->{_params},
  72            [qw( command max_edit_dist max_gap_opens max_gap_extns deln_protect_3p deln_protect_ends subseq_seed max_edit_dist_seed n_threads mm_penalty gap_open_penalty gap_extn_penalty subopt_hit_threshold trim_parameter )],
  73            "commands filtered by prefix");
  74
  75 my @a = @{$bwafac->_translate_params};
  76 is shift @a, 'aln', 'translate_params: command correct';
  77 my ($k, %h);
  78 for (@a) {
  79     (/^-/) ? ( $h{$k = $_} = undef ) : ( $h{$k} = $_ );
  80 }
  81 is_deeply( \%h, { '-R' => 35, '-t' => 1 }, 'translate_params: options correct');
  82
  83
  84 # test run_bwa filearg parsing
  85 # a pipeline...
  86
  87 SKIP : {
  88     test_skip( -requires_executable => $bwafac,
  89                -tests => 12 );
  90
  91     my %tmpfiles;
  92     for (qw(rd1 rd2 refseq sai1f sai2f samf bamf)) {
  93         $tmpfiles{$_} = test_output_file();
  94     }
  95
  96     my $rd1 = test_input_file('r1bwa.fq');
  97     my $rd2 = test_input_file('r2bwa.fq');
  98     my $refseq = test_input_file('Ft.frag.fas');
  99
 100     ok my $bwa = Bio::Tools::Run::BWA->new( -command => 'index' ), "make refseq index factory";
 101     ok $bwa->run_bwa( -fas => $refseq ), "index refseq";
 102     ok $bwa = Bio::Tools::Run::BWA->new( -command => 'aln' ), "make aln factory";
 103     ok $bwa->run_bwa( -fas => $refseq, -faq => $rd1, -sai => $tmpfiles{sai1f} ), "map read1 to refseq";
 104     ok $bwa->run_bwa( -fas => $refseq, -faq => $rd2, -sai => $tmpfiles{sai2f}), "map read 2 to refseq";
 105     ok $bwa = Bio::Tools::Run::BWA->new( -command => 'sampe' ), "paired read assembly factory";
 106     ok $bwa->run_bwa( -fas => $refseq, -sai1 => $tmpfiles{sai1f}, -faq1 => $rd1,
 107                       -sai2 => $tmpfiles{sai2f}, -faq2 => $rd2, -sam => $tmpfiles{samf}), "assemble paired reads";
 108
 109     #test run (assembly pipeline)
 110     ok $bwa = Bio::Tools::Run::BWA->new(), "make a full assembly factory";
 111     is ($bwa->command, 'run', "command attribute set");
 112     ok my $assy = $bwa->run($rd1, $refseq, $rd2), "make full assy";
 113
 114     TODO: {
 115         local $TODO = "latest bwa doesn't work with assembly pipeline";
 116         is ($assy->get_nof_contigs, 204, "number of contigs");
 117         is ($assy->get_nof_singlets, 220, "number of singlets");
 118     }
 119 }