scripts/best_verify.pl

   1 #!/usr/bin/perl -w
   2
   3 #
   4 # Verifies that Bowtie's --better and --best alignment modes really
   5 # produce alignments that are the best possible in terms of score.
   6 #
   7 # Run this from the Bowtie directory.
   8 #
   9 # E.g.: perl ../scripts/best_verify.pl -v 0
  10 #       perl ../scripts/best_verify.pl -v 0 e_coli genomes/NC_008253.fna reads/e_coli_1000.fq
  11 #
  12
  13 use warnings;
  14 use strict;
  15 use Getopt::Long;
  16
  17 my %options=();
  18 my $varg;
  19 my $narg;
  20 my $g = undef;
  21 my $l = undef;
  22 my $e = undef;
  23 my $C = undef;
  24 my $debug = 0;
  25 my $round = "";
  26 my $best = undef;
  27 my $nomaqround = 0;
  28 my $result = GetOptions(
  29         "v=i" => \$varg,
  30         "n=i" => \$narg,
  31         "l=i" => \$l,
  32         "e=i" => \$e,
  33         "d"   => \$debug,
  34         "g"   => \$g,
  35         "best"=> \$best,
  36         "C"   => \$C,
  37         "nomaqround" => \$nomaqround) || die "Bad options";
  38 my $match_mode = "-n 2";
  39 $match_mode  = " -v " . $varg if defined($varg);
  40 $match_mode  = " -n " . $narg if defined($narg);
  41 $match_mode .= " -l $l " if defined($l);
  42 $match_mode .= " -e $e " if defined($e);
  43 $match_mode .= " -g " if defined($g);
  44 $match_mode .= " -C " if defined($C);
  45 $match_mode .= " --cost ";
  46 $round = "--nomaqround" if $nomaqround;
  47 my $seedLen = 28;
  48 $seedLen = $l if defined($l);
  49
  50 print "Maq-like rounding is: ".($nomaqround ? "off" : "on") . "\n";
  51 print "Using match mode: $match_mode\n";
  52
  53 if(defined($C)) {
  54         print "Aborting with exitlevel 0 because colorspace alignment was selected\n";
  55         exit 0;
  56 }
  57
  58 my $bowtie_dir = ".";
  59 my $bowtie_exe = "bowtie";
  60 $bowtie_exe .= "-debug" if $debug;
  61
  62 my $index  = "e_coli";
  63 $index = $ARGV[0] if defined($ARGV[0]);
  64 my $reads = "reads/e_coli_1000.fq";
  65 $reads = $ARGV[1] if defined($ARGV[1]);
  66
  67 if($reads =~ /\.fa$/) {
  68         $reads = " -f $reads ";
  69 } elsif($reads =~ /\.raw$/) {
  70         $reads = " -r $reads ";
  71 } elsif($reads =~ /\.integer.fq$/) {
  72         $reads = " -q --integer-quals $reads ";
  73 } elsif($reads =~ /\.fq$/) {
  74         $reads = " -q $reads ";
  75 } else {
  76         $reads = " -c $reads ";
  77 }
  78
  79 my $vmode = ($match_mode =~ /[-]v/);
  80
  81 system("make -C $bowtie_dir $bowtie_exe") == 0 || die;
  82
  83 # Run Bowtie to get best alignments
  84 my $bowtie_best_cmd = "$bowtie_dir/$bowtie_exe -y $round $match_mode --best --refidx $index $reads";
  85
  86 # Run Bowtie to get all alignments
  87 my $bowtie_all_cmd = "$bowtie_dir/$bowtie_exe -y $round $match_mode -a --refidx $index $reads";
  88
  89 print "$bowtie_best_cmd\n";
  90 open BOWTIE_BEST, "$bowtie_best_cmd |";
  91 my %nameToBestScore = ();
  92 my %nameToBestAlignment = ();
  93 my $bestAls = 0;
  94 while(<BOWTIE_BEST>) {
  95         next if /^Reported/;
  96         chomp;
  97         my $line = $_;
  98         my @ls = split(/[\t]/, $line);
  99         my $btcost = $ls[-1];
 100         $#ls >= 5 || die "Alignment not formatted correctly: $line";
 101         my $name = $ls[0];
 102         defined($nameToBestAlignment{$name}) && die "Read with name $name appeared more than once in best-hit output";
 103         defined($nameToBestScore{$name}) && die "Read with name $name appeared more than once in best-hit output";
 104         my $len = length($ls[4]);
 105         my $quals = $ls[5];
 106         my $mmstr = "";
 107         $mmstr = $ls[7] if defined($ls[7]);
 108         my $cost = 0;
 109         my $fw = $ls[1] eq "+";
 110         $seedLen = $len if $vmode;
 111         if($mmstr ne "-") {
 112                 my @mms = split(/,/, $mmstr);
 113                 for my $mm (@mms) {
 114                         my @mmss = split(/:/, $mm);
 115                         my $mmoff = int($mmss[0]);
 116                         if($mmoff < $seedLen) {
 117                                 $cost += (1 << 14);
 118                         }
 119                         if(!$fw) {
 120                                 $mmoff = $len - $mmoff - 1;
 121                         }
 122                         my $q = substr($quals, $mmoff, 1);
 123                         $q = int(ord($q) - 33);
 124                         if(!$nomaqround) {
 125                                 $q = int(int($q + 5) / 10);
 126                                 $q = 3 if $q > 3;
 127                                 $q *= 10;
 128                         }
 129                         my $qcost += $q;
 130                         $cost += $qcost;
 131                 }
 132         }
 133         print "$line: $cost\n";
 134         $cost == $btcost || die "script-calculated cost $cost doesn't match bowtie-calculated cost $btcost\n";
 135         $nameToBestAlignment{$name} = $line;
 136         $nameToBestScore{$name} = $cost;
 137         $bestAls++;
 138 }
 139 close(BOWTIE_BEST);
 140 if($? != 0) {
 141         die "bowtie --best process failed with exitlevel $?\n";
 142 }
 143
 144 print "$bowtie_all_cmd\n";
 145 open BOWTIE_ALL, "$bowtie_all_cmd |";
 146 my $allAls = 0;
 147 while(<BOWTIE_ALL>) {
 148         next if /^Reported/;
 149         chomp;
 150         my $line = $_;
 151         my @ls = split(/[\t]/, $line);
 152         my $btcost = $ls[-1];
 153         $#ls >= 5 || die "Alignment not formatted correctly: $line";
 154         my $name = $ls[0];
 155         my $len = length($ls[4]);
 156         my $quals = $ls[5];
 157         my $mmstr = "";
 158         $mmstr = $ls[7] if defined($ls[7]);
 159         my $cost = 0;
 160         my $fw = $ls[1] eq "+";
 161         $seedLen = $len if $vmode;
 162         if($mmstr ne "-") {
 163                 my @mms = split(/,/, $mmstr);
 164                 for my $mm (@mms) {
 165                         my @mmss = split(/:/, $mm);
 166                         my $mmoff = int($mmss[0]);
 167                         if($mmoff < $seedLen) {
 168                                 $cost += (1 << 14);
 169                         }
 170                         if(!$fw) {
 171                                 $mmoff = $len - $mmoff - 1;
 172                         }
 173                         my $q = substr($quals, $mmoff, 1);
 174                         $q = int(ord($q) - 33);
 175                         if(!$nomaqround) {
 176                                 $q = int(int($q + 5) / 10);
 177                                 $q = 3 if $q > 3;
 178                                 $q *= 10;
 179                         }
 180                         my $qcost += $q;
 181                         $cost += $qcost;
 182                 }
 183         }
 184         print "$line: $cost\n";
 185         $cost == $btcost || die "script-calculated cost $cost doesn't match bowtie-calculated cost $btcost\n";
 186         defined($nameToBestAlignment{$name}) ||
 187                 die "Read with alignment:\n$line\nhas no corresponding alignment in best-hit mode\n";
 188         int($cost) >= int($nameToBestScore{$name}) ||
 189                 die "Alignment:\n$line\n$cost\nis better than:\n$nameToBestAlignment{$name}\n$nameToBestScore{$name}\n";
 190         $allAls++;
 191 }
 192 close(BOWTIE_ALL);
 193 if($? != 0) {
 194         die "bowtie -a process failed with exitlevel $?\n";
 195 }
 196
 197 print "Checked $bestAls best-alignments against $allAls all-alignments\n";
 198 print "PASSED\n";