scripts/random_bowtie_tests.pl

   1 #!/usr/bin/perl -w
   2
   3 #
   4 # Generate and run a series of random (but non-trivial) test cases for
   5 # the Bowtie suite of tools, including bowtie, bowtie-build and
   6 # bowtie-inspect.  Most of the problems turned up this way are in the
   7 # form of assertions in the tools.  However, we also do some sanity-
   8 # checking of the results; e.g. we use the pe_verify.pl script to
   9 # verify that paired-end alignments are consistent with matched-up
  10 # single-end alignments.
  11 #
  12 # Usage: perl random_tester.pl [rand seed] \
  13 #                              [# outer iters] \
  14 #                              [# inner iters] \
  15 #                              [min # text bases] \
  16 #                              [max text bases to add to min] \
  17 #                              [min # read bases] \
  18 #                              [max read bases to add to min]
  19 #
  20 # Options:
  21 #   -n         don't attempt to compile binaries; use existing binaries
  22 #   -p "<pol>" use only the specified alignment policy (e.g. "-n 3")
  23 #
  24
  25 use List::Util qw[min max];
  26 use Getopt::Long;
  27 use FindBin qw($Bin);
  28
  29 my $setPolicy;
  30 my $help = 0;
  31 my $noCompile = 0;
  32 my $bowtie = "./bowtie";                 # path to version of 'bowtie' binary to use
  33 my $bowtie_build = "./bowtie-build";     # path to version of 'bowtie-build' binary to use
  34 my $bowtie_build_old = "./bowtie-build"; # path to old version of 'bowtie-build' binary to use (for inspect.pl)
  35 my $bowtie_inspect = "./bowtie-inspect"; # path to version of 'bowtie-build' binary to use
  36
  37 GetOptions(
  38         "bowtie=s"           => \$bowtie,
  39         "bowtie-build=s"     => \$bowtie_build,
  40         "bowtie-build-old=s" => \$bowtie_build_old,
  41         "bowtie-inspect=s"   => \$bowtie_inspect,
  42         "h|help"             => \$help,
  43         "n|no-compile"       => \$noCompile,
  44         "p|policy"           => \$setPolicy, # overrides pickPolicy()
  45         "e|pairs-only"       => \$pairedEndOnly) || die "Bad arguments";
  46
  47 if($help) {
  48         print "Usage: perl random_bowtie_tests.pl [options]* seed outer inner tbase trand pbase prand\n";
  49         exit 0;
  50 }
  51
  52 my $seed = 0;
  53 $seed = int $ARGV[0] if defined($ARGV[0]);
  54 srand $seed;
  55
  56 # make all the relevant binaries, unless we were asked not to
  57 unless($noCompile) {
  58         run("make bowtie bowtie-debug bowtie-build-debug ".
  59                "bowtie-inspect-debug bowtie-inspect bowtie-build") == 0 || die "Error building";
  60 }
  61
  62 # Alignment policies
  63 my @policies = (
  64         "-n 3", "-n 2", "-n 1", "-n 0",
  65         "-v 3", "-v 2", "-v 1", "-v 0"
  66 );
  67
  68 sub pickPolicy {
  69         my $pe = shift;
  70         my $r = int(rand($#policies + 1));
  71         my $pol = $policies[$r];
  72         $pol = $setPolicy if defined($setPolicy);
  73         if($pe && int(rand(2)) == 0) {
  74                 my $min = int(rand(200))+10;
  75                 $pol .= " -I $min";
  76                 $pol .= " -X ".int(rand(30)+$min);
  77         }
  78         if(!$pe && int(rand(2)) == 0) {
  79                 # Possibly ask for best alignments
  80                 $pol .= " --best";
  81         }
  82         if($pol =~ /-n/ && int(rand(2)) == 0) {
  83                 $pol .= " --nomaqround";
  84         }
  85         if($pol =~ /-n/ && int(rand(2)) == 0) {
  86                 $pol .= " -l ".int(rand(30)+8);
  87         }
  88         if($pol =~ /-n/ && int(rand(2)) == 0) {
  89                 $pol .= " -e ".int(rand(120)+40);
  90         }
  91         return $pol;
  92 }
  93
  94 my $outer = 5000;
  95 $outer = int $ARGV[1] if defined($ARGV[1]);
  96 my $limit = $outer;
  97
  98 my $inner = 100;
  99 $inner = int $ARGV[2] if defined($ARGV[2]);
 100
 101 my $tbase = 100;
 102 $tbase = int $ARGV[3] if defined($ARGV[3]);
 103 my $trand = 300;
 104 $trand = int $ARGV[4] if defined($ARGV[4]);
 105
 106 my $pbase = 10;
 107 $pbase = int $ARGV[5] if defined($ARGV[5]);
 108 my $prand = 30;
 109 $prand = int $ARGV[6] if defined($ARGV[6]);
 110
 111 my $ibase = 50;
 112 $ibase = int $ARGV[7] if defined($ARGV[7]);
 113 my $irand = 250;
 114 $irand = int $ARGV[8] if defined($ARGV[8]);
 115
 116 my $verbose = 0;
 117 my $exitOnFail = 1;
 118 my @dnaMap = ('A', 'T', 'C', 'G',
 119               'N',
 120               'M', 'R', 'W', 'S', 'Y', 'K', 'V', 'H', 'D', 'B');
 121
 122 sub nonACGTtoN {
 123         my $t = shift;
 124         $t =~ tr/-NnMmRrWwSsYyKkVvHhDdBbXx/N/;
 125         $t =~ /[ACGTN]+/ || die "Bad N-ized DNA string: $t";
 126         return $t;
 127 }
 128
 129 sub randGap() {
 130         my $or = int(rand(4));
 131         my $gap = "";
 132         if($or == 0) {
 133                 my $ir = int(rand(100))+1;
 134                 if(($ir & 1) == 1) {
 135                         for(my $i = 0; $i < $ir; $i++) {
 136                                 $gap .= 'N';
 137                         }
 138                 } else {
 139                         for(my $i = 0; $i < $ir; $i++) {
 140                                 $gap .= '-';
 141                         }
 142                 }
 143         }
 144         return $gap;
 145 }
 146
 147 # Generates a random DNA string of the given length
 148 sub randDna($) {
 149         my $num = shift;
 150         my $i;
 151         my $t = '';
 152         my $noAmbigs = int(rand(3)) == 0;
 153         for($i = 0; $i < $num; $i++) {
 154                 my $or = int(rand(50));
 155                 if($or == 0 && !$noAmbigs) {
 156                         # Add a random, possibly ambiguous character
 157                         $t .= $dnaMap[int(rand($#dnaMap+1))];
 158                 } elsif($or == 1 && !$noAmbigs) {
 159                         # Add a random-length streak of Ns (max: 20)
 160                         my $streak = int(rand(20))+1;
 161                         for(my $j = $i; $j < $num && $j < $streak; $j++) {
 162                                 $t .= 'N';
 163                         }
 164                 } else {
 165                         # Add a random non-ambiguous character
 166                         $t .= $dnaMap[int(rand(4))];
 167                 }
 168         }
 169         return $t;
 170 }
 171
 172 sub randColor($) {
 173         my %cmap = ("A" => "0", "C" => "1", "G" => "2", "T" => "3", "N" => ".");
 174         my $s = nonACGTtoN(randDna(shift));
 175         for(my $i = 0; $i < length($s); $i++) {
 176                 defined($cmap{substr($s, $i, 1)}) || die "No entry for \"".substr($s, $i, 1)."\"\n";
 177                 substr($s, $i, 1) = $cmap{substr($s, $i, 1)};
 178         }
 179         return $s;
 180 }
 181
 182 # Utility function that returns the reverse complement of its argument
 183 sub revcomp {
 184         my $r = shift;
 185         my $c = shift; # may be undef
 186         $r = reverse($r);
 187         $r =~ tr/aAcCgGtT/tTgGcCaA/ unless (defined($c) && $c);
 188         return $r;
 189 }
 190
 191 # Utility function that returns the complement of its argument
 192 sub comp {
 193         my $r = shift;
 194         my $c = shift;
 195         $r =~ tr/aAcCgGtT/tTgGcCaA/ unless (defined($c) && $c);
 196         return $r;
 197 }
 198
 199 # Insert $src just before the $off'th character of $dst
 200 sub ins {
 201         my ($dst, $src, $off) = @_;
 202         $off < length($dst) || die "ins offset $off is not less than dst length ".length($dst)."\n";
 203         my $s = substr($dst, 0, $off) . $src . substr($dst, $off);
 204         return $s;
 205 }
 206
 207 # Delete the $off'th character of $dst
 208 sub del {
 209         my ($dst, $off) = @_;
 210         my $s = substr($dst, 0, $off).substr($dst, $off+1);
 211         return $s;
 212 }
 213
 214 # Substitute $src for the $off'th character of $dst
 215 sub subst {
 216         my ($dst, $src, $off) = @_;
 217         my $s = substr($dst, 0, $off).$src.substr($dst, $off+1);
 218         return $s;
 219 }
 220
 221 ##
 222 # Sanity check whether a read and a list of edits corresponds correctly
 223 # to a substring of the reference.
 224 #
 225 sub checkAlignmentRef {
 226         my ($ref, $read, $fw, $off, $edits, $alnuc) = @_;
 227         return unless $alnuc;
 228         my $orig = $read;
 229         $off-- unless $alnuc;
 230         if($edits ne '-' && $edits ne "") {
 231                 my $adjust = 0;
 232                 my @es = split(/[,]/, $edits);
 233                 for my $e (@es) {
 234                         my $colonoff = index($e, ":");
 235                         my $caratoff = index($e, ">");
 236                         my $pos = substr($e, 0, $colonoff);
 237                         my $ref = substr($e, $colonoff+1, $caratoff-$colonoff-1);
 238                         my $qry = substr($e, $caratoff+1);
 239                         length($qry) == 1 || die "Query char in edit $e, \"$qry\", isn't 1 char\n";
 240                         $ref = comp($ref, $alnuc) unless $fw;
 241                         if($qry eq "-") {
 242                                 # insertion
 243                                 #my $inspos = ($fw ? $pos + $adjust : length($read) - $pos - $adjust);
 244                                 my $inspos = $pos + $adjust;
 245                                 $read = ins($read, $ref, $inspos);
 246                                 $adjust += length($ref);
 247                         } elsif($ref eq "-") {
 248                                 # deletion
 249                                 #my $delpos = ($fw ? $pos + $adjust : length($read) - $pos - $adjust);
 250                                 my $delpos = $pos + $adjust;
 251                                 $read = del($read, $delpos);
 252                                 $adjust--;
 253                         } else {
 254                                 # mismatch
 255                                 #my $mmpos = ($fw ? $pos + $adjust : length($read) - $pos - $adjust);
 256                                 my $mmpos = $pos + $adjust;
 257                                 $read = subst($read, $ref, $mmpos);
 258                         }
 259                 }
 260         }
 261         $read = revcomp($read, $alnuc) unless $fw;
 262         my $rstr = substr($ref, $off, length($read));
 263         $fw = $fw ? '+' : '-';
 264         $read eq $rstr || die "\n\nAlignment failed to pass sanity check:\n".
 265                               "FW: $fw, Off: $off, Edits: $edits\n".
 266                               "Orig: $orig\n".
 267                               "Qry:  $read\n".
 268                               "Ref:  $rstr\n";
 269 }
 270
 271 ##
 272 # Given a string in nucleotide space, convert to colorspace.
 273 #
 274 sub colorize($$) {
 275         my ($s, $nucs) = @_;
 276         defined($s) || die;
 277         my %cmap = (
 278                 "AA" => "0", "CC" => "0", "GG" => "0", "TT" => "0",
 279                 "AC" => "1", "CA" => "1", "GT" => "1", "TG" => "1",
 280                 "AG" => "2", "GA" => "2", "CT" => "2", "TC" => "2",
 281                 "AT" => "3", "TA" => "3", "CG" => "3", "GC" => "3",
 282                 "NA" => ".", "NC" => ".", "NG" => ".", "NT" => ".",
 283                 "AN" => ".", "CN" => ".", "GN" => ".", "TN" => ".",
 284                 "NN" => "."
 285         );
 286         my %nmap = ("0" => "A", "1" => "C", "2" => "G", "3" => "T", "." => "N");
 287         my $ret = "";
 288         for(my $i = 0; $i < length($s)-1; $i++) {
 289                 my $di = uc substr($s, $i, 2);
 290                 $di =~ tr/-NnMmRrWwSsYyKkVvHhDdBbXx/N/;
 291                 defined($cmap{$di}) || die "Bad dinuc: $di\n";
 292                 $ret .= ($nucs ? $nmap{$cmap{$di}} : $cmap{$di});
 293         }
 294         return $ret;
 295 }
 296
 297 ##
 298 # Encode colors as proxy nucleotides
 299 #
 300 sub nucencode($) {
 301         my $s = shift;
 302         my %nmap = ("0" => "A", "1" => "C", "2" => "G", "3" => "T", "." => "N",
 303                     "A" => "A", "C" => "C", "G" => "G", "T" => "T", "N" => "N");
 304         my $ret = "";
 305         for(my $i = 0; $i < length($s); $i++) {
 306                 my $c = uc substr($s, $i, 1);
 307                 defined($nmap{$c}) || die;
 308                 $ret .= $nmap{$c};
 309         }
 310         return $ret;
 311 }
 312
 313 # Add some random quality values to encourage excercising the
 314 # backtracking code
 315 sub addQual($) {
 316         my $r = shift;
 317         my $len = length($r);
 318         $r .= ":";
 319         for(my $i = 0; $i < $len; $i++) {
 320                 my $c = "-";
 321                 while(not $c =~ /[0-9A-Z\/=@%]/) {
 322                         $c = chr(33 + int(rand(41)));
 323                 }
 324                 $r .= $c;
 325         }
 326         return $r;
 327 }
 328
 329 # Trim whitespace from a string argument
 330 sub trim($) {
 331         my $string = shift;
 332         $string =~ s/^\s+//;
 333         $string =~ s/\s+$//;
 334         return $string;
 335 }
 336
 337 sub run {
 338         my $cmd = shift;
 339         open CMDTMP, ">.randtmp$seed.cmd" || die;
 340         print CMDTMP "$cmd\n";
 341         close(CMDTMP);
 342         return system($cmd);
 343 }
 344
 345 sub runBacktick {
 346         my $cmd = shift;
 347         open CMDTMP, ">.randtmp$seed.cmd" || die;
 348         print CMDTMP "$cmd\n";
 349         close(CMDTMP);
 350         return `$cmd`;
 351 }
 352
 353 # Build an Ebwt based on given arguments
 354 sub build {
 355         my($t, $color, $offRate, $ftabChars) = @_;
 356         my $ret = 0;
 357
 358         my $file1 = "-c";
 359         my $file2 = "\"$t\"";
 360         if(substr($t, 0, 1) eq '-') {
 361                 # Add backslash to escape first dash
 362                 $file2 = "\"\\$t\"";
 363         }
 364
 365         $color = ($color ? "-C" : "");
 366
 367         # Write reference sequences to a FASTA file
 368         open FA, ">.randtmp$seed.fa" || die "Could not open temporary fasta file";
 369         my @seqs = split(/,/, $t);
 370         for(my $i = 0; $i <= $#seqs; $i++) {
 371                 print FA ">$i\n";
 372                 print FA "$seqs[$i]\n";
 373         }
 374         close(FA);
 375
 376         $cmd = "perl $Bin/test/inspect.pl --debug --seed $seed --bowtie-build2=$bowtie_build_old --ref=.randtmp$seed.fa";
 377         print "$cmd\n";
 378         system($cmd) == 0
 379                 || die "inspect.pl died with exitlevel $?";
 380
 381         # Make a version of the FASTA file where all non-A/C/G/T characters
 382         # are Ns.  This is useful if we'd like to compare to the output of
 383         # bowtie-inspect.
 384         #
 385         # We make one version (FAN) that's colorized when the index is
 386         # colorspace so that we can compare the output of bowtie-inspect -e,
 387         # and another version (FANO) that's always in nucleotide space, so
 388         # we can test default bowtie-inspect.
 389         #
 390         open FAN, ">.randtmp$seed.ns.fa" || die "Could not open temporary fasta file";
 391         open FANO, ">.randtmp$seed.ns.orig.fa" || die "Could not open temporary fasta file";
 392         for(my $i = 0; $i <= $#seqs; $i++) {
 393                 print FAN ">$i\n";
 394                 print FANO ">$i\n";
 395                 my $t = nonACGTtoN($seqs[$i]);
 396                 defined($t) || die;
 397                 print FANO "$t\n";
 398                 $t = colorize($t, 1) if $color;
 399                 print FAN "$t\n";
 400         }
 401         close(FAN);
 402         close(FANO);
 403
 404         my $fasta = int(rand(2)) == 0;
 405         if($fasta) {
 406                 # Use the FASTA file as input
 407                 $file1 = "-f";
 408                 $file2 = ".randtmp$seed.fa";
 409         }
 410
 411         my $bucketArg = "";
 412         my $bucketRand = int(rand(3));
 413         if($bucketRand == 0) {
 414                 $bucketArg = "--bmaxdivn ";
 415                 $bucketArg .= (int(rand(30))+1);
 416         } elsif($bucketRand == 1) {
 417                 $bucketArg = "-a ";
 418         }
 419
 420         $offRate   = "--offrate $offRate"   if $offRate ne "";
 421         $ftabChars = "--ftabchars $ftabChars" if $ftabChars ne "";
 422         my $noauto = "";
 423         $noauto = "-a" if $offRate eq "" && $ftabChars eq "";
 424         my $args = "-q --sanity $color $file1 $noauto $offRate $ftabChars $bucketArg $file2";
 425
 426         # Do unpacked version
 427         my $cmd = "${bowtie_build}-debug $args .tmp$seed";
 428         run("echo \"$cmd\" > .tmp$seed.cmd");
 429         print "$cmd\n";
 430         my $out = trim(runBacktick("$cmd 2>&1"));
 431         $out =~ s/Warning: Encountered reference sequence with only gaps//g;
 432         $out = trim($out);
 433         if($out eq "") {
 434                 $ret++;
 435         } else {
 436                 print "Expected no output, got:\n$out\n";
 437                 if($exitOnFail) {
 438                         exit 1;
 439                 }
 440         }
 441
 442         # Do packed version and assert that it matches unpacked version
 443         # (sometimes, but not all the time because it takes a while)
 444         if(int(rand(4)) == 0) {
 445                 $cmd = "${bowtie_build}-debug -a -p $args .tmp$seed.packed";
 446                 print "$cmd\n";
 447                 $out = trim(runBacktick("$cmd 2>&1"));
 448                 $out =~ s/Warning: Encountered reference sequence with only gaps//g;
 449                 $out = trim($out);
 450                 if($out eq "") {
 451                         if(run("diff .tmp$seed.1.ebwt .tmp$seed.packed.1.ebwt") != 0) {
 452                                 die if $exitOnFail;
 453                         } elsif(run("diff .tmp$seed.2.ebwt .tmp$seed.packed.2.ebwt") != 0) {
 454                                 die if $exitOnFail;
 455                         } else {
 456                                 $ret++;
 457                         }
 458                 } else {
 459                         print "Expected no output, got:\n$out\n";
 460                         if($exitOnFail) {
 461                                 exit 1;
 462                         }
 463                 }
 464         }
 465
 466         return $ret;
 467 }
 468
 469 sub deleteReadParts {
 470         run("rm -f .tmp.un$seed". ".* .tmp.un$seed". "_1.* .tmp.un$seed". "_2.*");
 471         run("rm -f .tmp.max$seed".".* .tmp.max$seed"."_1.* .tmp.max$seed"."_2.*");
 472         run("rm -f .tmp.al$seed". ".* .tmp.al$seed". "_1.* .tmp.al$seed". "_2.*");
 473 }
 474
 475 sub search($$$$$$$$$$) {
 476         my($tstr, $cstr, $pe, $color, $p1, $p2, $policy, $oneHit, $requireResult, $offRate) = @_;
 477         my $ret = doSearch($tstr, $cstr, $pe, $color, $p1, $p2, $policy, $oneHit, $requireResult, $offRate);
 478         deleteReadParts();
 479         return $ret;
 480 }
 481
 482 ##
 483 # Make sure that a verbose alignment is consistent with the content of
 484 # the reference
 485 #
 486 # $alnuc: Alignments are in nucleotides (not colors, as with -C --col-cseq)
 487 #
 488 sub checkRefVerbose($$$$) {
 489         my ($l, $textsr, $ctextsr, $alnuc) = @_;
 490         my @texts = @{$textsr};
 491         my @ctexts = @{$ctextsr};
 492         scalar(@texts) > 0 || die;
 493         my @ls = split(/\t/, $l);
 494         my ($fw, $ref, $off, $seq, $mms) = ($ls[1] eq '+', $ls[2], $ls[3], $ls[4], $ls[7]);
 495         $mms = "-" unless(defined($mms));
 496         defined($ref) || die "Malformed verbose output line: $l\n";
 497         $ref == int($ref) || die;
 498         $ref <= $#texts || die "Ref idx $ref exceeds number of texts ".($#texts+1)."\n";
 499         my $orig = $seq;
 500         $orig = revcomp($orig, $alnuc) unless $fw;
 501         print STDERR "Checking alignemnt: $l\n";
 502         checkAlignmentRef($alnuc ? $texts[$ref] : $ctexts[$ref],
 503                           $orig, $fw, $off, $mms, $alnuc);
 504         return 1;
 505 }
 506
 507 # Search for a pattern in an existing Ebwt
 508 sub doSearch($$$$$$$$$$) {
 509         my($nstr, $cstr, $pe, $color, $p1, $p2, $policy, $oneHit, $requireResult, $offRate) = @_;
 510
 511         my @nts = split(/,/, $nstr); # nucleotide texts
 512         my @cts = split(/,/, $cstr); # color texts
 513
 514         my $patarg = "-c";
 515         my $patstr = "\"$p1\"";
 516         $patstr = "-1 \"$p1\" -2 \"$p2\"" if $pe;
 517         my $outfile = ".tmp$seed.out";
 518
 519         my $alnuc = 1;
 520         if($color) {
 521                 $color = "-C";
 522                 if(int(rand(3))) {
 523                         $color .= " --col-cseq ";
 524                         $alnuc = 0;
 525                 }
 526                 $color .= " --col-cqual" if int(rand(3)) == 0;
 527         }
 528
 529         my $ext = ".fq";
 530         my $format = int(rand(5));
 531         if($format == 0) {
 532                 # FASTA
 533                 open FA, ">.randread$seed.fa" || die "Could not open temporary fasta file";
 534                 $ext = ".fa";
 535                 my @cs = split /[,]/, $p1;
 536                 my $idx = 0;
 537                 foreach my $c (@cs) {
 538                         my @cms = split /[:]/, $c;
 539                         print FA ">$idx\n$cms[0]\n";
 540                         $idx++;
 541                 }
 542                 close FA;
 543                 $patarg = "-f";
 544                 $patstr = ".randread$seed.fa";
 545                 if($pe) {
 546                         run("mv .randread$seed.fa .randread$seed"."_1.fa");
 547                         $patstr = "-1 .randread$seed"."_1.fa";
 548                         open FA, ">.randread$seed"."_2.fa" || die "Could not open temporary fasta file";
 549                         @cs = split /[,]/, $p2;
 550                         my $idx = 0;
 551                         foreach my $c (@cs) {
 552                                 my @cms = split /[:]/, $c;
 553                                 print FA ">$idx\n$cms[0]\n";
 554                                 $idx++;
 555                         }
 556                         close FA;
 557                         $patstr .= " -2 .randread$seed"."_2.fa";
 558                 }
 559         } elsif($format == 1) {
 560                 # FASTQ with ASCII qualities
 561                 open FQ, ">.randread$seed.fq" || die "Could not open temporary fastq file";
 562                 $ext = ".fq";
 563                 my @cs = split /[,]/, $p1;
 564                 my $idx = 0;
 565                 foreach my $c (@cs) {
 566                         my @cms = split /[:]/, $c;
 567                         $#cms == 1 || die "Should be pair with : separating seq from quals: $c\n$p1\n";
 568                         print FQ "\@$idx\n$cms[0]\n+\n$cms[1]\n";
 569                         $idx++;
 570                 }
 571                 close FQ;
 572                 $patarg = "-q";
 573                 $patstr = ".randread$seed.fq";
 574                 if($pe) {
 575                         run("mv .randread$seed.fq .randread$seed"."_1.fq");
 576                         $patstr = "-1 .randread$seed"."_1.fq";
 577                         open FQ, ">.randread$seed"."_2.fq" || die "Could not open temporary fastq file";
 578                         @cs = split /[,]/, $p2;
 579                         my $idx = 0;
 580                         foreach my $c (@cs) {
 581                                 my @cms = split /[:]/, $c;
 582                                 $#cms == 1 || die "Should be pair with : separating seq from quals: $c\n$p2\n";
 583                                 print FQ "\@$idx\n$cms[0]\n+\n$cms[1]\n";
 584                                 $idx++;
 585                         }
 586                         close FQ;
 587                         $patstr .= " -2 .randread$seed"."_2.fq";
 588                 }
 589         } elsif($format == 2) {
 590                 # FASTQ with integer qualities
 591                 open FQ, ">.randread$seed.integer.fq" || die "Could not open temporary solexa fastq file";              my $ext = ".fa";
 592                 $ext = ".fq";
 593                 my @cs = split /[,]/, $p1;
 594                 my $idx = 0;
 595                 foreach my $c (@cs) {
 596                         my @cms = split /[:]/, $c;
 597                         print FQ "\@$idx\n$cms[0]\n+\n";
 598                         for(my $i = 0; $i < length($cms[1]); $i++) {
 599                                 my $q = substr($cms[1], $i, 1);
 600                                 $q = ord($q) - 33;
 601                                 print FQ "$q ";
 602                         }
 603                         print FQ "\n";
 604                         $idx++;
 605                 }
 606                 close FQ;
 607                 $patarg = "-q --integer-quals";
 608                 $patstr = ".randread$seed.integer.fq";
 609                 if($pe) {
 610                         run("mv .randread$seed.integer.fq .randread$seed.integer_1.fq");
 611                         $patstr = "-1 .randread$seed.integer_1.fq";
 612                         open FQ, ">.randread$seed.integer_2.fq" || die "Could not open temporary fastq file";
 613                         @cs = split /[,]/, $p2;
 614                         my $idx = 0;
 615                         foreach my $c (@cs) {
 616                                 my @cms = split /[:]/, $c;
 617                                 print FQ "\@$idx\n$cms[0]\n+\n";
 618                                 for(my $i = 0; $i < length($cms[1]); $i++) {
 619                                         my $q = substr($cms[1], $i, 1);
 620                                         $q = ord($q) - 33;
 621                                         print FQ "$q ";
 622                                 }
 623                                 $idx++;
 624                                 print FQ "\n";
 625                         }
 626                         close FQ;
 627                         $patstr .= " -2 .randread$seed.integer_2.fq";
 628                 }
 629         } elsif($format == 3) {
 630                 # Raw
 631                 open RAW, ">.randread$seed.raw" || die "Could not open temporary raw file";
 632                 $ext = ".raw";
 633                 my @cs = split /[,]/, $p1;
 634                 foreach my $c (@cs) {
 635                         my @cms = split /[:]/, $c;
 636                         print RAW "$cms[0]\n";
 637                 }
 638                 close RAW;
 639                 $patarg = "-r";
 640                 $patstr = ".randread$seed.raw";
 641                 if($pe) {
 642                         run("mv .randread$seed.raw .randread$seed"."_1.raw");
 643                         $patstr = "-1 .randread$seed"."_1.raw";
 644                         open RAW, ">.randread$seed"."_2.raw" || die "Could not open temporary fastq file";
 645                         @cs = split /[,]/, $p2;
 646                         foreach my $c (@cs) {
 647                                 my @cms = split /[:]/, $c;
 648                                 print RAW "$cms[0]\n";
 649                         }
 650                         close RAW;
 651                         $patstr .= " -2 .randread$seed"."_2.raw";
 652                 }
 653         }
 654
 655         # Perhaps dump unaligned reads using --un argument
 656         my $unalignArg = "";
 657         my $unalignReconArg = "";
 658         my $unalign = int(rand(5));
 659         if($unalign == 0 || $unalign == 2) {
 660                 $unalignArg .= "--un .tmp.un$seed$ext ";
 661                 if($pe) {
 662                         $unalignReconArg .= " .tmp.un$seed"."_1$ext .tmp.un$seed"."_2$ext";
 663                 } else {
 664                         $unalignReconArg .= " .tmp.un$seed$ext";
 665                 }
 666         }
 667         if($unalign == 1 || $unalign == 2) {
 668                 $unalignArg .= "--max .tmp.max$seed$ext ";
 669                 if($unalign == 2) {
 670                         if($pe) {
 671                                 $unalignReconArg .= " .tmp.max$seed"."_1$ext .tmp.max$seed"."_2$ext";
 672                         } else {
 673                                 $unalignReconArg .= " .tmp.max$seed$ext";
 674                         }
 675                 }
 676         }
 677         if($unalign == 2 || $unalign == 3) {
 678                 $unalignArg .= "--al .tmp.al$seed$ext ";
 679                 if($unalign == 2) {
 680                         if($pe) {
 681                                 $unalignReconArg .= " .tmp.al$seed"."_1$ext .tmp.al$seed"."_2$ext";
 682                         } else {
 683                                 $unalignReconArg .= " .tmp.al$seed$ext";
 684                         }
 685                 }
 686         }
 687
 688         my $khits = "-k 1";
 689         my $mhits = 0;
 690         if(int(rand(2)) == 0) {
 691                 $khits = "-a";
 692         } else {
 693                 $khits = "-k " . (int(rand(20))+2);
 694         }
 695         if(int(rand(3)) == 0) {
 696                 $khits .= " --strata --best";
 697         }
 698         if(int(rand(2)) == 0) {
 699                 $requireResult = 0;
 700                 $mhits = (int(rand(20))+2);
 701         }
 702
 703         if($mhits > 0) {
 704                 if(int(rand(2)) == 0) {
 705                         $khits .= " -m $mhits";
 706                 } else {
 707                         $khits .= " -M $mhits";
 708                 }
 709         }
 710
 711         my $strand = "";
 712         if(int(rand(4)) == 0) {
 713                 if(int(rand(2)) == 0) {
 714                         # Reverse-complement reference strand only
 715                         $strand = "--nofw";
 716                 } else {
 717                         # Forward reference strand only
 718                         $strand = "--norc";
 719                 }
 720         }
 721
 722         if($oneHit || 1) {
 723                 $oneHit = "";
 724         } else {
 725                 $oneHit = "-a";
 726         }
 727         my $offRateStr = "";
 728         if($offRate ne "" && int(rand(3)) == 0) {
 729                 $offRate == int($offRate) || die "Bad offrate: $offRate\n";
 730                 $offRateStr = "--offrate " . ($offRate + 1 + int(rand(4)));
 731         }
 732         defined($policy) || die;
 733         defined($color) || die;
 734         defined($strand) || die;
 735         defined($unalignArg) || die;
 736         defined($khits) || die;
 737         defined($offRateStr) || die;
 738         defined($nstr) || die;
 739         my $cmd = "${bowtie}-debug $policy $color $strand $unalignArg $khits $offRateStr --cost --orig \"$nstr\" $oneHit --sanity $patarg .tmp$seed $patstr";
 740         print "$cmd\n";
 741         my $out = trim(runBacktick("$cmd 2>.tmp$seed.stderr | tee .tmp$seed.stdout"));
 742
 743         # Bad exitlevel?
 744         if($? != 0) {
 745                 print "Exitlevel: $?\n";
 746                 if($exitOnFail) {
 747                         my $err = `cat .tmp$seed.stderr 2> /dev/null`;
 748                         print "Stdout:\n$out\nStderr:\n$err\n";
 749                         exit 1;
 750                 }
 751                 return 0;
 752         }
 753         my $err = `cat .tmp$seed.stderr 2> /dev/null`;
 754
 755         # No output?
 756         if($out eq "" && $requireResult) {
 757                 print "Expected results but got \"No Results\"\n";
 758                 if($exitOnFail) {
 759                         print "Stdout:\n$out\nStderr:\n$err\n";
 760                         exit 1;
 761                 }
 762                 return 0;
 763         }
 764
 765         # Parse output to see if any of it is bad
 766         print $out;
 767         my @outlines = split(/[\r\n]+/, $out);
 768         my %outhash = ();
 769         my %readcount = ();
 770         my %readStratum = (); # for unpaired reads
 771         my %readStratum1 = (); # for mate 1
 772         my %readStratum2 = (); # for mate 2
 773         my $lastread = "";
 774
 775         for(my $i = 0; $i <= $#outlines; $i++) {
 776                 # Get the alignment for the first mate (or the unpaired read)
 777                 my $l = $outlines[$i];
 778                 chomp($l);
 779                 checkRefVerbose($l, \@nts, \@cts, $alnuc);
 780                 my $key = "$l";
 781                 my $l2 = "";
 782                 if($pe) {
 783                         # Get the alignment for the second mate
 784                         $l2 = $outlines[++$i];
 785                         defined($l2) || die "Odd number of output lines";
 786                         chomp($l2);
 787                         checkRefVerbose($l2, \@nts, \@cts, $alnuc);
 788                         $key .= ", $l2";
 789                 }
 790                 print "$key\n";
 791                 # No two results should be the same
 792                 if(!$pe) {
 793                         !defined($outhash{$key}) || die "Result $key appears in output twice";
 794                 }
 795                 $outhash{$key} = 1;
 796                 # Parse out the read id
 797                 my @ls = split(/[\t]/, $l);
 798                 my ($mate1, $fw1, $stratum1) = ($ls[0], $ls[1], $ls[8]);
 799                 my ($mate2, $fw2, $stratum2) = (undef, undef, undef);
 800                 if($pe) {
 801                         my @l2s = split(/[\t]/, $l2);
 802                         ($mate2, $fw2, $stratum2) = ($l2s[0], $l2s[1], $l2s[8]);
 803                 }
 804                 my $peor = ($color ? "ff" : "fr");
 805                 my $m1fw = ($peor eq 'fr' || $peor eq 'ff') ? '+' : '-';
 806                 my $m2fw = ($peor eq 'ff') ? '+' : '-';
 807                 if($pe && $peor eq 'ff') {
 808                         $fw1 eq $fw2 || die "Saw different orientations for mates\n";
 809                 } elsif($pe) {
 810                         $fw1 ne $fw2 || die "Saw same orientation for mates\n";
 811                 }
 812                 if($strand =~ /nofw/) {
 813                         if($pe) {
 814                                 my $m = substr($mate1, index($mate1, "/")+1);
 815                                 if($peor eq 'ff') {
 816                                         ($fw1 eq '-' && $fw2 eq '-') ||
 817                                                 die "Saw non-rc alignment with --nofw specified\n";
 818                                 } else {
 819                                         "$m$fw1" eq "2+" ||
 820                                                 die "Saw a forward alignment on line ".($i+1)." when --nofw was specified";
 821                                 }
 822                         } else {
 823                                 $fw1 eq "-" ||
 824                                         die "Saw a forward alignment on line ".($i+1)." when --nofw was specified";
 825                         }
 826                 } elsif($strand =~ /norc/) {
 827                         if($pe) {
 828                                 if($peor eq 'ff') {
 829                                         ($fw1 eq '+' && $fw2 eq '+') ||
 830                                                 die "Saw non-fw alignment with --nofw specified\n";
 831                                 } else {
 832                                         my $m = substr($mate1, index($mate1, "/")+1);
 833                                         "$m$fw1" eq "1+" ||
 834                                                 die "Saw a rev-comp alignment on line ".($i+1)." when --norc was specified";
 835                                 }
 836                         } else {
 837                                 $fw1 eq "+" ||
 838                                         die "Saw a rev-comp alignment on line ".($i+1)." when --norc was specified";
 839                         }
 840                 }
 841                 if($mate1 ne $lastread && !$pe) {
 842                         die "Read $mate1 appears multiple times non-consecutively" if defined($readcount{$mate1});
 843                 }
 844                 if(!$pe && $policy =~ /--strata /) {
 845                         !defined($readStratum{$mate1}) ||
 846                                 $readStratum{$mate1} == $stratum1 ||
 847                                 die "Incompatible strata: old: $readStratum{$mate1}, cur: $stratum1";
 848                         $readStratum{$mate1} = $stratum1;
 849                 } elsif($policy =~ /--strata /) {
 850                         # Paired-end case; use minimum of 2 mates
 851                         if($mate2 =~ /\/1$/) {
 852                                 # Swticheroo
 853                                 my $tmp = $mate1; $mate1 = $mate2; $mate2 = $tmp;
 854                                 $tmp = $stratum1; $stratum1 = $stratum2; $stratum2 = $tmp;
 855                         }
 856                         if(defined($readStratum1{$mate1})) {
 857                                 # One or ther other mate has to have the same stratum as we recorded previously
 858                                 defined($readStratum2{$mate2}) || die "$mate1 was seen before, but not $mate2";
 859                                 $readStratum1{$mate1} == $stratum1 ||
 860                                 $readStratum2{$mate2} == $stratum2 ||
 861                                 $readStratum1{$mate1} == $stratum2 ||
 862                                 $readStratum2{$mate2} == $stratum1 ||
 863                                         die "Incompatible strata: old: <$readStratum1{$mate1}, $readStratum2{$mate2}> cur: <$stratum1, $stratum2>";
 864                         }
 865                         $readStratum1{$mate1} = $stratum1;
 866                         $readStratum2{$mate2} = $stratum2;
 867                 }
 868                 $lastread = $mate1;
 869                 $readcount{$mate1}++;
 870                 if($mhits > 0) {
 871                         $readcount{$mate1} <= $mhits || die "Read $mate1 matched more than $mhits times";
 872                 }
 873         }
 874
 875         {
 876                 $cmd = "${bowtie} $policy $color $strand $unalignArg $khits $offRateStr --cost --orig \"$nstr\" $oneHit --sanity $patarg .tmp$seed $patstr";
 877                 print "$cmd\n";
 878                 my $out2 = trim(runBacktick("$cmd 2>.tmp$seed.stderr"));
 879                 $out2 eq $out || die "Normal bowtie output did not match debug bowtie output";
 880
 881                 $cmd = "${bowtie} $policy $color $strand $unalignArg $khits --orig \"$nstr\" $oneHit --sanity $patarg .tmp$seed $patstr $outfile";
 882                 print "$cmd\n";
 883                 my $out3 = trim(runBacktick("$cmd 2>.tmp$seed.stderr"));
 884
 885                 $cmd = "${bowtie} --mm $policy $color $strand $unalignArg $khits --orig \"$nstr\" $oneHit --sanity $patarg .tmp$seed $patstr $outfile";
 886                 print "$cmd\n";
 887                 my $out4 = trim(runBacktick("$cmd 2>.tmp$seed.stderr"));
 888                 $out3 eq $out4 || die "Normal bowtie output did not match memory-mapped bowtie output";
 889         }
 890
 891         # Now do another run with verbose output so that we can check the
 892         # mismatch strings
 893         $cmd = "${bowtie}-debug $policy $color $strand $unalignArg $khits $offRateStr --orig \"$nstr\" $oneHit --sanity $patarg .tmp$seed $patstr";
 894         print "$cmd\n";
 895         $out = trim(runBacktick("$cmd 2>.tmp$seed.stderr"));
 896         # Parse output to see if any of it is bad
 897         @outlines = split('\n', $out);
 898         for my $l (@outlines) {
 899                 checkRefVerbose($l, \@nts, \@cts, $alnuc);
 900         }
 901
 902         if($pe) {
 903                 # If search was for paired-end alignments, then verify the
 904                 # outcome using pe_verify.pl
 905                 my $pol = $policy;
 906                 $pol =~ s/--.*//;
 907                 my $patstr2 = $patstr;
 908                 $patstr2 =~ s/-1//;
 909                 $patstr2 =~ s/-2//;
 910                 my $col = ($color ? "-C" : "");
 911                 open TMP, ">.tmp$seed.pe_verify.cmd" || die;
 912                 $cmd = "perl scripts/pe_verify.pl --args=\"--quiet\" $col -d $pol .tmp$seed $patstr2";
 913                 print "$cmd\n";
 914                 print TMP "$cmd\n";
 915                 $out = trim(runBacktick("$cmd 2>.tmp$seed.pe_verify.stderr"));
 916                 close(TMP);
 917                 # Bad exitlevel?
 918                 if($? != 0) {
 919                         print "scripts/pe_verify.pl exitlevel: $?\n";
 920                         if($exitOnFail) {
 921                                 my $err = `cat .tmp$seed.pe_verify.stderr 2> /dev/null`;
 922                                 print "Stdout:\n$out\nStderr:\n$err\n";
 923                                 exit 1;
 924                         }
 925                         return 0;
 926                 }
 927         }
 928
 929         if(!$pe && $policy =~ /--best/) {
 930                 my $col = ($color ? "-C" : "");
 931                 $cmd = "perl scripts/best_verify.pl -d $policy $col .tmp$seed $patstr";
 932                 print "$cmd\n";
 933                 $out = trim(runBacktick("$cmd 2>.tmp$seed.best_verify.stderr"));
 934
 935                 if($? != 0) {
 936                         print "scripts/best_verify.pl exitlevel: $?\n";
 937                         if($exitOnFail) {
 938                                 my $err = `cat .tmp$seed.best_verify.stderr 2> /dev/null`;
 939                                 print "Stdout:\n$out\nStderr:\n$err\n";
 940                                 exit 1;
 941                         }
 942                         return 0;
 943                 }
 944         }
 945
 946         # Run the reconciler to ensure that --un, --max, and --al had
 947         # sane output
 948         # .tmp$seed.verbose.out
 949         if($format >= 0 && $format <= 2 && $unalignReconArg ne "") {
 950                 deleteReadParts();
 951                 my $cmd = "${bowtie} $policy $color $strand $unalignArg $khits $offRateStr --orig \"$nstr\" $oneHit --sanity $patarg .tmp$seed $patstr .tmp$seed.verbose.out";
 952                 print "$cmd\n";
 953                 run($cmd) == 0 || die "Error performing reconciler run\n";
 954                 $khits =~ s/--strata --best//;
 955                 $patarg =~ s/--integer-quals//;
 956                 $unalignReconArg =~ /\.tmp\.un/ || die;
 957                 if($pe) {
 958                         $patstr =~ s/-1//;
 959                         $patstr =~ s/-2//;
 960                         $cmd = "perl scripts/reconcile_alignments_pe.pl $color $patarg $khits $patstr .tmp$seed.verbose.out $unalignReconArg";
 961                         print "$cmd\n";
 962                         run($cmd) == 0 || die "Failed to reconcile";
 963                 } else {
 964                         $cmd = "perl scripts/reconcile_alignments.pl $color $patarg $khits $patstr .tmp$seed.verbose.out $unalignReconArg";
 965                         print "$cmd\n";
 966                         run($cmd) == 0 || die "Failed to reconcile";
 967                 }
 968         }
 969
 970         # Success
 971         return 1;
 972 }
 973
 974 my $pass = 0;
 975 my $tests = 0;
 976 my $fail = 0;
 977
 978 for(; $outer > 0; $outer--) {
 979
 980         # Generate random parameters
 981         my $offRate = "";
 982         if(int(rand(2)) == 0) {
 983                 $offRate = int(rand(16));
 984         }
 985         my $ftabChars = "";
 986         if(int(rand(2)) == 0) {
 987                 $ftabChars = 1 + int(rand(8));
 988         }
 989
 990         # Generate random text(s)
 991         my $nt = int(rand(10)) + 1;
 992         my $tstr = '', $cstr = '';
 993         my @nts;
 994         my @cts;
 995
 996         # For each reference sequence...
 997         for(my $i = 0; $i < $nt; $i++) {
 998                 my $tlen = $tbase + int(rand($trand));
 999                 # Add some all-gap references to try to confuse bowtie-build
1000                 # and bowtie-inspect
1001                 my $tt;
1002                 if(rand() < 0.15) {
1003                         $tt = 'N' x $tlen;
1004                         push(@nts, $tt);
1005                         push(@cts, colorize($tt, 1));
1006                         my $newt = randGap() . $tt . randGap();
1007                         $tstr .= $newt;
1008                         $cstr .= colorize($newt, 1);
1009                         $tstr .= ",";
1010                         $cstr .= ",";
1011                 }
1012                 $tt = randDna($tlen); # add text meat
1013                 push(@nts, $tt);
1014                 push(@cts, colorize($tt, 1));
1015                 my $newt = randGap() . $tt . randGap();
1016                 $tstr .= $newt;
1017                 $cstr .= colorize($newt, 1);
1018                 if($i < $nt-1) {
1019                         $tstr .= ",";
1020                         $cstr .= ",";
1021                 }
1022         }
1023
1024         my $color = (int(rand(2)) == 0);
1025
1026         # Run the command to build the Ebwt from the random text
1027         $pass += build($tstr, $color, $offRate, $ftabChars);
1028         last if(++$tests > $limit);
1029
1030         my $in = $inner;
1031         for(; $in >= 0; $in--) {
1032                 # Paired-end?
1033                 my $pe = (int(rand(2)) == 0);
1034                 $pe = 1 if $pairedEndOnly;
1035                 # Generate random pattern(s) based on text
1036                 my $pfinal1 = '';
1037                 my $pfinal2 = '';
1038                 # Decide number of patterns to generate
1039                 my $np = int(rand(30)) + 1;
1040                 for(my $i = 0; $i < $np; $i++) {
1041                         # Pick a text
1042                         my $tt = $nts[int(rand($#nts))];
1043                         # Pick a length
1044                         my $pl;
1045                         my $plen;
1046                         my $pr;
1047                         my $p1;
1048                         my $p2;
1049                         if($pe) {
1050                                 # Pick a left-hand offset for the insert
1051                                 my $il = int(rand(length($tt))) - 10;
1052                                 # Pick a length for the insert
1053                                 my $ilen = int(rand($irand)) + $ibase;
1054                                 $il = min($il, length($tt)-$ilen);
1055                                 my $ir = min($il + $ilen, length($tt));
1056                                 my $is = substr $tt, $il, $ir - $il;
1057                                 # Pick a length for the #1 mate
1058                                 my $plen1 = int(rand($prand)) + $pbase;
1059                                 # Pick a length for the #2 mate
1060                                 my $plen2 = int(rand($prand)) + $pbase;
1061                                 $p1 = substr $is, 0, $plen1;
1062                                 $is = reverse $is;
1063                                 $p2 = substr $is, 0, $plen2;
1064                                 $p2 = reverse $p2;
1065                                 $p2 = revcomp($p2);
1066                         } else {
1067                                 # Pick a length for the read
1068                                 $plen = int(rand($prand)) + $pbase;
1069                                 $pl = int(rand(length($tt))) - 10;
1070                                 $pl = max($pl, $color ? 5 : 4);
1071                                 $pl = min($pl, length($tt));
1072                                 $pr = min($pl + $plen, length($tt));
1073                                 $p1 = substr $tt, $pl, $pr - $pl;
1074                         }
1075                         # Check for empty pattern or pattern that spans a comma
1076                         if(length($p1) < ($color ? 5 : 4) || index($p1, ",") != -1) {
1077                                 $i--; next;
1078                         }
1079                         if($pe && (length($p2) < ($color ? 5 : 4) || index($p2, ",") != -1)) {
1080                                 $i--; next;
1081                         }
1082                         # Optionally add nucleotide changes to pattern
1083                         if($i > 0) {
1084                                 my $nummms = int(rand($color ? 3 : 5));
1085                                 for(my $j = 0; $j < $nummms; $j++) {
1086                                         substr($p1, int(rand(length($p1))), 1) = randDna(1);
1087                                 }
1088                                 if($pe) {
1089                                         $nummms = int(rand(4));
1090                                         for(my $j = 0; $j < $nummms; $j++) {
1091                                                 substr($p2, int(rand(length($p2))), 1) = randDna(1);
1092                                         }
1093                                 }
1094                         }
1095                         # Possibly reverse complement it
1096                         if((int(rand(2)) == 0)) {
1097                                 $p1 = revcomp($p1);
1098                                 if($pe) {
1099                                         $p2 = revcomp($p2);
1100                                         my $ptmp = $p1;
1101                                         $p1 = $p2;
1102                                         $p2 = $ptmp;
1103                                 }
1104                         }
1105                         $p1 =~ tr/MRWSYKVHDBX/N/;
1106                         if($color) {
1107                                 defined($p1) || die;
1108                                 !$pe || defined($p2) || die;
1109                                 my ($nize1, $nize2) = (int(rand(2)) == 0, int(rand(2)) == 0);
1110                                 $p1 = colorize($p1, $nize1);
1111                                 $p2 = colorize($p2, $nize2) if $pe;
1112                                 # Optionally add color changes to pattern
1113                                 if($i > 0) {
1114                                         my $nummms = int(rand(3));
1115                                         for(my $j = 0; $j < $nummms; $j++) {
1116                                                 my $r = randColor(1);
1117                                                 $r = nucencode($r) if $nize1;
1118                                                 substr($p1, int(rand(length($p1))), 1) = $r;
1119                                         }
1120                                         if($pe) {
1121                                                 $nummms = int(rand(4));
1122                                                 for(my $j = 0; $j < $nummms; $j++) {
1123                                                         my $r = randColor(1);
1124                                                         $r = nucencode($r) if $nize2;
1125                                                         substr($p2, int(rand(length($p2))), 1) = $r;
1126                                                 }
1127                                         }
1128                                 }
1129                         }
1130                         # Add valid random quality values
1131                         $p1 = addQual($p1);
1132                         $pfinal1 .= $p1;
1133                         $pfinal1 .= "," if($i < $np-1);
1134                         if($pe) {
1135                                 $p2 =~ tr/MRWSYKVHDBX/N/;
1136                                 $p2 = addQual($p2);
1137                                 $pfinal2 .= $p2;
1138                                 $pfinal2 .= "," if($i < $np-1);
1139                         }
1140                 }
1141
1142                 # Run the command to search for the pattern from the Ebwt
1143                 my $oneHit = (int(rand(3)) == 0);
1144                 my $policy = pickPolicy($pe);
1145                 my $expectResult = 1;
1146                 if(!$pe) {
1147                         for(my $i = 0; $i < length($pfinal1); $i++) {
1148                                 my $c = substr($pfinal1, $i, 1);
1149                                 last if($c eq ',');
1150                                 if($c ne 'A' && $c ne 'C' && $c ne 'G' && $c ne 'T') {
1151                                         $expectResult = 0;
1152                                         last;
1153                                 }
1154                         }
1155                 } else {
1156                         $expectResult = 0;
1157                 }
1158                 $pass += search($tstr, $cstr, $pe, $color, $pfinal1, $pfinal2, $policy, $oneHit, $expectResult, $offRate); # require 1 or more results
1159                 last if(++$tests > $limit);
1160         }
1161
1162         $in = $inner;
1163         for(; $in >= 0; $in--) {
1164                 my $pe = (int(rand(2)) == 0);
1165                 # Generate random pattern *not* based on text
1166                 my $pfinal1 = '';
1167                 my $pfinal2 = '';
1168                 my $np = int(rand(10)) + 1;
1169                 for(my $i = 0; $i < $np; $i++) {
1170                         my $plen = int(rand($prand)) + $pbase;
1171                         my $p1 = randDna($plen);
1172                         $plen = int(rand($prand)) + $pbase;
1173                         my $p2 = randDna($plen);
1174                         $p1 =~ tr/MRWSYKVHDBX/N/;
1175                         $p1 = colorize($p1, int(rand(2)) == 0) if $color;
1176                         $p1 = addQual($p1);
1177                         $p2 =~ tr/MRWSYKVHDBX/N/ if $pe;
1178                         $p2 = colorize($p2, int(rand(2)) == 0) if $color && $pe;
1179                         $p2 = addQual($p2) if $pe;
1180                         $pfinal1 .= $p1;
1181                         $pfinal2 .= $p2 if $pe;
1182                         if($i < $np-1) {
1183                                 $pfinal1 .= ",";
1184                                 $pfinal2 .= "," if $pe;
1185                         }
1186                 }
1187
1188                 # Run the command to search for the pattern from the Ebwt
1189                 my $oneHit = (int(rand(3)) == 0);
1190                 my $policy = pickPolicy($pe);
1191                 $pass += search($tstr, $cstr, $pe, $color, $pfinal1, $pfinal2, $policy, $oneHit, 0, $offRate); # do not require any results
1192                 last if(++$tests > $limit);
1193         }
1194 }
1195
1196 print "$pass tests passed, $fail failed\n";
1197 exit 1 if $fail > 0;