t/Bowtie.t

   1 #-*-perl-*-
   2
   3
   4 use strict;
   5 use warnings;
   6 no warnings qw(once);
   7 our $home;
   8 our $ulimit;
   9 BEGIN {
  10     $home = '.';    # set to '.' for Build use,
  11                     # '..' for debugging from .t file
  12     unshift @INC, $home;
  13     use Bio::Root::Test;
  14
  15     if ($^O =~ m/mswin/i) {
  16         # "ulimit" does not exists and "ulimit -n 3>&-" does not work
  17         # as expected in Windows, so set $ulimit to a fixed reasonable number
  18         $ulimit = 1000;
  19     }
  20     else {
  21         eval { $ulimit = qx/ulimit -n 3>&-/ };
  22         if ($@ || !defined($ulimit)) {
  23             # skip all run tests, we can't ensure the ulimit is high enough for
  24             # these tests (needs ulimit -n of ~1000)
  25             $ulimit = 0;
  26         } else {
  27             chomp $ulimit
  28         }
  29     }
  30     my $DEBUG = test_debug();
  31     print STDERR $ulimit if $DEBUG == 1;
  32
  33     test_begin(-tests => 73,
  34                -requires_modules => [qw(IPC::Run Bio::Tools::Run::Bowtie)]);
  35
  36 }
  37
  38 use Bio::Tools::Run::WrapperBase;
  39 use Bio::SeqIO;
  40
  41 # test command functionality
  42
  43 ok my $bowtiefac = Bio::Tools::Run::Bowtie->new(
  44     -command            => 'paired',
  45     -try_hard           => 1,
  46     -min_insert_size    => 300,
  47     -solexa             => 1,
  48     -max_mismatches     => 4
  49     ), "make a factory using command 'assemble'";
  50 # ParameterBaseI compliance : really AssemblerBase tests...
  51 ok $bowtiefac->parameters_changed, "parameters changed on construction";
  52 ok $bowtiefac->min_insert_size, "access parameter";
  53 ok !$bowtiefac->parameters_changed, "parameters_changed cleared on read";
  54 ok $bowtiefac->set_parameters( -trim5 => 10 ), "set a param not set in constructor";
  55 ok $bowtiefac->parameters_changed, "parameters_changed set";
  56 is ($bowtiefac->trim5, 10, "parameter really set");
  57 is ($bowtiefac->min_insert_size, 300, "original parameter unchanged");
  58 ok !$bowtiefac->parameters_changed, "parameters_changed cleared on read";
  59 ok $bowtiefac->set_parameters( -min_insert_size => 100 ), "change an original parameter";
  60 is ($bowtiefac->min_insert_size, 100, "parameter really changed");
  61 ok $bowtiefac->reset_parameters( -min_insert_size => 200 ), "reset parameters with arg";
  62 ok !$bowtiefac->max_mismatches, "original parameters undefined";
  63 is ($bowtiefac->min_insert_size, 200, "parameter really reset via arg");
  64 ok $bowtiefac->set_parameters( -raw => 1 ), "set an exclusive group parameter";
  65 ok $bowtiefac->raw, "parameter really set";
  66 ok $bowtiefac->set_parameters( -fastq => 1 ), "set an incompatible parameter";
  67 ok $bowtiefac->fastq, "parameter really set";
  68 ok !$bowtiefac->raw, "original exclusive parameter really unset";
  69
  70 $bowtiefac->set_parameters(
  71     -command            => 'paired',
  72     -try_hard           => 1,
  73     -suppress           => 1000,
  74     -max_mismatches     => 4
  75     );
  76 ok $bowtiefac->parameters_changed, "parameters changed";
  77
  78 is( scalar $bowtiefac->available_parameters, 68, "all available options");
  79 is( scalar $bowtiefac->available_parameters('params'), 33, "available parameters" );
  80 is( scalar $bowtiefac->available_parameters('switches'), 35, "available switches" );
  81 #back to beginning - but with single
  82 $bowtiefac = Bio::Tools::Run::Bowtie->new(
  83     -command            => 'single',
  84     -try_hard           => 1,
  85     -min_insert_size    => 300,
  86     -solexa             => 1,
  87     -max_mismatches     => 4
  88     );
  89 ok $bowtiefac->parameters_changed, "parameters changed";
  90
  91 is( scalar $bowtiefac->available_parameters, 60, "all available options");
  92 is( scalar $bowtiefac->available_parameters('params'), 28, "available parameters" );
  93 is( scalar $bowtiefac->available_parameters('switches'), 32, "available switches" );
  94 my %pms = $bowtiefac->get_parameters;
  95 is_deeply( \%pms,
  96                 { command            => 'single',
  97                   max_mismatches     => 4,
  98                   solexa             => 1,
  99                   try_hard           => 1,
 100                   sam_format         => 1}, "get_parameters correct"); # we are single so filter 300 out
 101                                                                        # and again, we default to SAM
 102 is( $bowtiefac->command, 'single', "command attribute set");
 103
 104 is_deeply( $bowtiefac->{_options}->{_commands},
 105            [@Bio::Tools::Run::Bowtie::program_commands],
 106            "internal command array set" );
 107
 108 is_deeply( $bowtiefac->{_options}->{_prefixes},
 109            {%Bio::Tools::Run::Bowtie::command_prefixes},
 110            "internal prefix hash set");
 111
 112 is_deeply( $bowtiefac->{_options}->{_params},
 113            [qw( command qualities skip upto trim5 trim3 max_seed_mismatches
 114                 max_qual_mismatch max_quality_sum snp_penalty snp_frac
 115                 seed_length max_mismatches max_backtracks max_search_ram
 116                 report_n_alignments supress supress_random offset_base
 117                 defaul_mapq sam_rg suppress_columns alignmed_file
 118                 unaligned_file excess_file threads offrate random_seed )],
 119            "commands filtered by prefix");
 120
 121 my @a = @{$bowtiefac->_translate_params};
 122 my ($k, %h);
 123 for (@a) {
 124     (/^-/) ? ( $h{$k = $_} = undef ) : ( $h{$k} = $_ );
 125 }
 126 is_deeply( \%h, { '-v' => 4, '--solexa-quals' => undef, '-y' => undef, '-S' => undef }, 'translate_params: options correct'); # we default to SAM so '-S' appears
 127
 128
 129 # test run_bowtie filearg parsing
 130
 131 SKIP : {
 132     test_skip( -requires_executable => $bowtiefac,
 133                -tests => 40 ); # three tests not included due to absence of SAM functionality at this stage
 134     skip("; set 'ulimit -n' to 1000 for bowtie tests", 40) unless $ulimit >= 1000;
 135     my $rdr = test_input_file('bowtie', 'reads', 'e_coli_1000.raw');
 136     my $rda = test_input_file('bowtie', 'reads', 'e_coli_1000.fa');
 137     my $rdq = test_input_file('bowtie', 'reads', 'e_coli_1000.fq');
 138     my $rdc = test_input_file('bowtie', 'reads', 'e_coli.cb');
 139     my $rda1 = test_input_file('bowtie', 'reads', 'e_coli_1000_1.fa');
 140     my $rda2 = test_input_file('bowtie', 'reads', 'e_coli_1000_2.fa');
 141     my $rdq1 = test_input_file('bowtie', 'reads', 'e_coli_1000_1.fq');
 142     my $rdq2 = test_input_file('bowtie', 'reads', 'e_coli_1000_2.fq');
 143     my $refseq = test_input_file('bowtie', 'indexes', 'e_coli');
 144     my @inlstr;
 145     my @inlobj;
 146     my $in = Bio::SeqIO->new( -file => $rda, -format => 'Fasta' );
 147     while ( my $seq = $in->next_seq() ) {
 148         push @inlstr,$seq->seq();
 149         push @inlobj,$seq;
 150     }
 151     my $inlstr = "'".join(',',@inlstr)."'";
 152
 153     # unpaired reads
 154     ok $bowtiefac = Bio::Tools::Run::Bowtie->new(
 155         -command            => 'single'
 156         ), "make unpaired reads bowtie factory";
 157
 158     $bowtiefac->set_parameters( -raw => 1 );
 159     ok $bowtiefac->_run( -ind => $refseq,
 160                          -seq => $rdr ), "read raw sequence";
 161
 162     like($bowtiefac->stderr, qr/reads processed: 1000/, "bowtie success");
 163
 164     $bowtiefac->reset_parameters( -raw => 0 );
 165     $bowtiefac->set_parameters( -fasta => 1 );
 166     ok $bowtiefac->_run( -ind => $refseq,
 167                          -seq => $rda ), "read fasta sequence";
 168
 169     like($bowtiefac->stderr, qr/reads processed: 1000/, "bowtie success");
 170
 171     $bowtiefac->reset_parameters( -fasta => 0 );
 172     $bowtiefac->set_parameters( -fastq => 1 );
 173     ok $bowtiefac->_run( -ind => $refseq,
 174                          -seq => $rdq ), "read fastq sequence";
 175
 176     like($bowtiefac->stderr, qr/reads processed: 1000/, "bowtie success");
 177
 178
 179     # paired reads
 180     ok $bowtiefac = Bio::Tools::Run::Bowtie->new(
 181         -command            => 'paired'
 182         ), "make paired reads bowtie factory";
 183
 184     $bowtiefac->set_parameters( -fasta => 1 );
 185     ok $bowtiefac->_run( -ind => $refseq,
 186                          -seq => $rda1,  -seq2 => $rda2 ), "read paired fasta sequence";
 187
 188     like($bowtiefac->stderr, qr/reads processed: 1000|Reported 1000/, "bowtie success");
 189
 190     $bowtiefac->reset_parameters( -fasta => 0 );
 191     $bowtiefac->set_parameters( -fastq => 1 );
 192     ok $bowtiefac->_run( -ind => $refseq,
 193                          -seq => $rdq1,  -seq2 => $rdq2 ), "read paired fastq sequence";
 194
 195     like($bowtiefac->stderr, qr/reads processed: 1000|Reported 1000/, "bowtie success");
 196
 197
 198     # test single
 199     # these parms are the bowtie defaults - getting raw is not default for module
 200     ok $bowtiefac = Bio::Tools::Run::Bowtie->new(
 201         -command             => 'single',
 202         -max_seed_mismatches => 2,
 203         -seed_length         => 28,
 204         -max_qual_mismatch   => 70,
 205         -want                => 'raw'
 206         ), "make a single alignment factory";
 207
 208     is( $bowtiefac->command, 'single', "command attribute set");
 209     is( $bowtiefac->max_seed_mismatches, 2, "seed mismatch param set");
 210     is( $bowtiefac->seed_length, 28, "seed length param set");
 211     is( $bowtiefac->max_qual_mismatch, 70, "quality mismatch param set");
 212     is( $bowtiefac->want, 'raw', "return type set");
 213     my $sam;
 214     $bowtiefac->set_parameters( -fastq => 1 );
 215     ok $sam = $bowtiefac->run($rdq, $refseq), "make file based alignment";
 216     ok eval { (-e $sam)&&(-r _) }, "make readable output";
 217     open (FILE, $sam);
 218     my $lines =()= <FILE>;
 219     close FILE;
 220     is( $lines, 1003, "number of alignments");
 221     SKIP: {
 222         test_skip( -tests => 4,
 223                -requires_module => 'Bio::DB::Sam');
 224         is($bowtiefac->want( 'Bio::Assembly::Scaffold' ), 'Bio::Assembly::Scaffold', "change mode");
 225         ok my $assy = $bowtiefac->run($rdq, $refseq), "make alignment";
 226         is( $assy->get_nof_contigs, 4, "number of contigs");
 227         is( $assy->get_nof_singlets, 691, "number of singlets");
 228     }
 229     # tests from here may fail due to insufficient memory - works with >=2GB
 230     # test crossbow
 231     # these parms are again the bowtie defaults - getting raw is not default for module
 232     SKIP: {
 233         skip("Crossbow functionality broken", 5);
 234         ok $bowtiefac = Bio::Tools::Run::Bowtie->new(
 235             -command             => 'crossbow',
 236             -max_seed_mismatches => 2,
 237             -seed_length         => 28,
 238             -max_qual_mismatch   => 70
 239         ), "make a crossbow alignment factory";
 240
 241         is( $bowtiefac->command, 'crossbow', "command attribute set");
 242         ok $sam = $bowtiefac->run($rdc, $refseq), "make file based alignment";
 243         ok eval { (-e $sam)&&(-r _) }, "make readable output";
 244         open (FILE, $sam);
 245         $lines =()= <FILE>;
 246         close FILE;
 247         is( $lines, 6, "number of alignments"); # 3 alignments and 3 SAM header lines
 248     }
 249
 250     SKIP: {
 251         skip("Error with bowtie v1.2.1: Reads file contained a pattern with more than 1024 quality values.", 10);
 252         ok $bowtiefac = Bio::Tools::Run::Bowtie->new(
 253             -command             => 'single',
 254             -max_seed_mismatches => 2,
 255             -seed_length         => 28,
 256             -max_qual_mismatch   => 70
 257             ), "make a single alignment factory";
 258
 259         ok $bowtiefac->set_parameters( -inline => 1 );
 260         ok $bowtiefac->_run( -ind => $refseq,
 261                              -seq => $inlstr ), "read sequence as strings in memory";
 262
 263         like($bowtiefac->stderr, qr/reads processed: 1000/, "bowtie success");
 264
 265         ok $bowtiefac->run( \@inlstr, $refseq ), "read sequence as seq objects";
 266
 267         like($bowtiefac->stderr, qr/reads processed: 1000/, "bowtie success");
 268
 269         ok $bowtiefac->run( \@inlobj, $refseq ), "read sequence as seq objects";
 270
 271         like($bowtiefac->stderr, qr/reads processed: 1000/, "bowtie success");
 272
 273         $bowtiefac->set_parameters( -inline => 0 );
 274         ok $sam = $bowtiefac->run($rdr,$refseq), "make variable based alignment";
 275
 276         like($bowtiefac->stderr, qr/reads processed: 1000/, "bowtie success");
 277     }
 278 }
 279
 280 1;