4 # Verifies that Bowtie's --better and --best alignment modes really
5 # produce alignments that are the best possible in terms of score.
7 # Run this from the Bowtie directory.
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
28 my $result = GetOptions(
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;
48 $seedLen = $l if defined($l);
50 print "Maq-like rounding is: ".($nomaqround ? "off" : "on") . "\n";
51 print "Using match mode: $match_mode\n";
54 print "Aborting with exitlevel 0 because colorspace alignment was selected\n";
59 my $bowtie_exe = "bowtie";
60 $bowtie_exe .= "-debug" if $debug;
63 $index = $ARGV[0] if defined($ARGV[0]);
64 my $reads = "reads/e_coli_1000.fq";
65 $reads = $ARGV[1] if defined($ARGV[1]);
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 ";
76 $reads = " -c $reads ";
79 my $vmode = ($match_mode =~ /[-]v/);
81 system("make -C $bowtie_dir $bowtie_exe") == 0 || die;
83 # Run Bowtie to get best alignments
84 my $bowtie_best_cmd = "$bowtie_dir/$bowtie_exe -y $round $match_mode --best --refidx $index $reads";
86 # Run Bowtie to get all alignments
87 my $bowtie_all_cmd = "$bowtie_dir/$bowtie_exe -y $round $match_mode -a --refidx $index $reads";
89 print "$bowtie_best_cmd\n";
90 open BOWTIE_BEST, "$bowtie_best_cmd |";
91 my %nameToBestScore = ();
92 my %nameToBestAlignment = ();
94 while(<BOWTIE_BEST>) {
98 my @ls = split(/[\t]/, $line);
100 $#ls >= 5 || die "Alignment not formatted correctly: $line";
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]);
107 $mmstr = $ls[7] if defined($ls[7]);
109 my $fw = $ls[1] eq "+";
110 $seedLen = $len if $vmode;
112 my @mms = split(/,/, $mmstr);
114 my @mmss = split(/:/, $mm);
115 my $mmoff = int($mmss[0]);
116 if($mmoff < $seedLen) {
120 $mmoff = $len - $mmoff - 1;
122 my $q = substr($quals, $mmoff, 1);
123 $q = int(ord($q) - 33);
125 $q = int(int($q + 5) / 10);
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;
141 die "bowtie --best process failed with exitlevel $?\n";
144 print "$bowtie_all_cmd\n";
145 open BOWTIE_ALL, "$bowtie_all_cmd |";
147 while(<BOWTIE_ALL>) {
151 my @ls = split(/[\t]/, $line);
152 my $btcost = $ls[-1];
153 $#ls >= 5 || die "Alignment not formatted correctly: $line";
155 my $len = length($ls[4]);
158 $mmstr = $ls[7] if defined($ls[7]);
160 my $fw = $ls[1] eq "+";
161 $seedLen = $len if $vmode;
163 my @mms = split(/,/, $mmstr);
165 my @mmss = split(/:/, $mm);
166 my $mmoff = int($mmss[0]);
167 if($mmoff < $seedLen) {
171 $mmoff = $len - $mmoff - 1;
173 my $q = substr($quals, $mmoff, 1);
174 $q = int(ord($q) - 33);
176 $q = int(int($q + 5) / 10);
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";
194 die "bowtie -a process failed with exitlevel $?\n";
197 print "Checked $bestAls best-alignments against $allAls all-alignments\n";