6 # Calculate mapability of each reference position assuming reference is
21 my $btargs = "-t -S --sam-nohead -M 1 --mm";
24 if(defined($ENV{BOWTIE_HOME})) {
25 $bowtie = "$ENV{BOWTIE_HOME}/bowtie";
26 unless(-x $bowtie) { $bowtie = "" };
29 $bowtie = `which bowtie 2>/dev/null`;
31 unless(-x $bowtie) { $bowtie = "" };
33 $bowtie = "./bowtie" if ($bowtie eq "" && -x "./bowtie");
38 "frequency=i" => \$freq,
39 "bowtie=s" => \$bowtie_arg,
45 print STDERR "Bowtie: found: $bowtie; given: $bowtie_arg\n";
46 print STDERR "Input fasta: $fa\n";
47 print STDERR "Index: $idx\n";
48 print STDERR "Alignment policy: $pol\n";
49 print STDERR "Window size: $win\n";
50 print STDERR "Frequency: $freq\n";
52 $fa ne "" || die "Must specify -fasta\n";
53 $idx ne "" || die "Must specify -idx\n";
54 -f "$idx.1.ebwt" || die "Could not find -idx index file $idx.1.ebwt\n";
56 $bowtie = $bowtie_arg if $bowtie_arg ne "";
59 if($bowtie_arg ne "") {
60 die "Specified -bowtie, \"$bowtie\" doesn't exist or isn't executable\n";
62 die "bowtie couldn't be found in BOWTIE_HOME, PATH, or current directory; please specify -bowtie\n";
67 my $name = ""; # name of sequence currently being processed
73 # Read lengths of all the entries in all the input fasta files.
77 my @is = split(/[,]/, $ins);
79 open IN, "$i" || die "Could not open $i\n";
83 if(substr($_, 0, 1) eq '>') {
84 next if /\?[0-9]*$/; # Skip >?50000 lines
86 $name = substr($name, 1); # Chop off >
87 if($name =~ /^FW:/ || $name =~ /^RC:/) {
88 $name = substr($name, 3); # Chop off FW:/RC:
90 my @ns = split(/\s+/, $name);
91 $name = $ns[0]; # Get short name
93 print STDERR "Saw name $name\n";
96 $lens{$name} += length($_); # Update length
97 $totlen += length($_);
103 print STDERR "Reading fasta lengths\n";
105 print STDERR " read ".scalar(keys %lens)." fasta sequences with total length $totlen\n";
108 for(my $i = 0; $i < $win; $i++) { push @last, 0 };
110 for(my $i = 0; $i < $win; $i++) { $last[$i] = 0 };
113 print STDERR "Opening bowtie pipe\n";
114 my $cmd = "$bowtie -F $win,$freq $btargs $pol $idx $fa";
115 print STDERR "Forward command: $cmd\n";
116 open BT, "$cmd |" || die "Couldn't open pipe '$cmd |'\n";
118 print STDERR "Reading...\n";
125 my @s = split(/\t/, $_);
126 my @s1 = split(/_/, $s[0]);
128 my $cname = join("_", @s1[0..$#s1-1]);
130 # If a read aligns, XM:i is not printed
131 # If a read fails to align, XM:i:0 is printed
132 # If a read aligns multiple places, XM:i:N is printed where N>0
133 my $mapable = ($s[-1] =~ /XM:i:/ ? 0 : 1);
135 $cname =~ s/\s.*//; # trim everything after first whitespace to get short name
136 if($cname =~ /^FW:/ || $cname =~ /^RC:/) {
137 $cname = substr($cname, 3);
139 if($name ne $cname) {
141 # Flush remaining characters from previous name
142 defined($lens{$name}) || die;
143 $cur == $lens{$name} - $win + 1 || die "name is $name, cur is $cur, len is $lens{$name}, win is $win\n";
144 for(; $cur < $lens{$name}; $cur++) {
145 $running -= $last[$cur % $win];
146 $last[$cur % $win] = 0;
147 $lastc = chr($running + 64);
149 if((($cur+1) % 60) == 0) {
156 defined($lens{$name}) || die "No such name as \"$name\"\n";
157 print "\n" unless $lastc eq "\n";
165 $running -= $last[$cur % $win];
166 $last[$cur % $win] = $mapable;
167 $running += $last[$cur % $win];
169 $running <= $win || die "running counter $running got higher than window size $win\n";
170 $lastc = chr($running + 64);
172 if((($cur+1) % 60) == 0) {
179 defined($lens{$name}) || die;
180 for(; $cur < $lens{$name}; $cur++) {
181 $running -= $last[$cur % $win];
182 $last[$cur % $win] = 0;
183 $lastc = chr($running + 64);
185 if((($cur+1) % 60) == 0) {
191 $? == 0 || die "Bad exitlevel from forward bowtie: $?\n";