scripts/bs_mapability.pl

   1 #!/usr/bin/perl -w
   2
   3 ##
   4 # mapability.pl
   5 #
   6 # Calculate mapability of each reference position assuming reference is
   7 # bisulfite treated.
   8 #
   9
  10 use strict;
  11 use warnings;
  12 use Getopt::Long;
  13
  14 my $fa = "";
  15 my $win = 50;
  16 my $freq = 1;
  17 my $bowtie = "";
  18 my $bowtie_arg = "";
  19 my $pol = "-v 3";
  20 my $fwidx = "";
  21 my $rcidx = "";
  22 my $btargs = "-t --norc -S --sam-nohead -M 1 --mm";
  23 my $debug = 0;
  24
  25 if(defined($ENV{BOWTIE_HOME})) {
  26         $bowtie = "$ENV{BOWTIE_HOME}/bowtie";
  27         unless(-x $bowtie) { $bowtie = "" };
  28 }
  29 if($bowtie eq "") {
  30         $bowtie = `which bowtie 2>/dev/null`;
  31         chomp($bowtie);
  32         unless(-x $bowtie) { $bowtie = "" };
  33 }
  34 $bowtie = "./bowtie" if ($bowtie eq "" && -x "./bowtie");
  35
  36 GetOptions(
  37         "fasta=s"     => \$fa,
  38         "window=i"    => \$win,
  39         "frequency=i" => \$freq,
  40         "bowtie=s"    => \$bowtie_arg,
  41         "policy=s"    => \$pol,
  42         "fwidx=s"     => \$fwidx,
  43         "rcidx=s"     => \$rcidx,
  44         "debug"       => \$debug
  45 ) || die;
  46
  47 print STDERR "Bowtie: found: $bowtie; given: $bowtie_arg\n";
  48 print STDERR "Input fasta: $fa\n";
  49 print STDERR "FW index: $fwidx\n";
  50 print STDERR "RC index: $rcidx\n";
  51 print STDERR "Alignment policy: $pol\n";
  52 print STDERR "Window size: $win\n";
  53 print STDERR "Frequency: $freq\n";
  54
  55 $fa ne "" || die "Must specify -fasta\n";
  56 $fwidx ne "" || die "Must specify -fwidx\n";
  57 $rcidx ne "" || die "Must specify -rcidx\n";
  58 -f "$fwidx.1.ebwt" || die "Could not find -fwidx index file $fwidx.1.ebwt\n";
  59 -f "$rcidx.1.ebwt" || die "Could not find -rcidx index file $rcidx.1.ebwt\n";
  60
  61 $bowtie = $bowtie_arg if $bowtie_arg ne "";
  62 unless(-x $bowtie) {
  63         # No bowtie? die
  64         if($bowtie_arg ne "") {
  65                 die "Specified -bowtie, \"$bowtie\" doesn't exist or isn't executable\n";
  66         } else {
  67                 die "bowtie couldn't be found in BOWTIE_HOME, PATH, or current directory; please specify -bowtie\n";
  68         }
  69 }
  70
  71 my $running = 0;
  72 my $name = ""; # name of sequence currently being processed
  73 my %lens = ();
  74 my @names = ();
  75 my $totlen = 0;
  76
  77 ##
  78 # Read lengths of all the entries in all the input fasta files.
  79 #
  80 sub readLens($) {
  81         my $ins = shift;
  82         my @is = split(/[,]/, $ins);
  83         for my $i (@is) {
  84                 open IN, "$i" || die "Could not open $i\n";
  85                 my $name = "";
  86                 while(<IN>) {
  87                         chomp;
  88                         if(substr($_, 0, 1) eq '>') {
  89                                 next if /\?[0-9]*$/; # Skip >?50000 lines
  90                                 $name = $_;
  91                                 $name = substr($name, 1); # Chop off >
  92                                 if($name =~ /^FW:/ || $name =~ /^RC:/) {
  93                                         $name = substr($name, 3); # Chop off FW:/RC:
  94                                 }
  95                                 my @ns = split(/\s+/, $name);
  96                                 $name = $ns[0]; # Get short name
  97                                 push @names, $name;
  98                                 print STDERR "Saw name $name\n";
  99                         } else {
 100                                 $name ne "" || die;
 101                                 $lens{$name} += length($_); # Update length
 102                                 $totlen += length($_);
 103                         }
 104                 }
 105                 close(IN);
 106         }
 107 }
 108 print STDERR "Reading fasta lengths\n";
 109 readLens($fa);
 110 print STDERR "  read ".scalar(keys %lens)." sequences with total length $totlen\n";
 111
 112 my @last;
 113 for(my $i = 0; $i < $win; $i++) { push @last, 0 };
 114 sub clearLast {
 115         for(my $i = 0; $i < $win; $i++) { $last[$i] = 0 };
 116 }
 117
 118 print STDERR "Opening bowtie pipes\n";
 119 my $fwcmd = "$bowtie -F $win,$freq $btargs $pol $fwidx $fa";
 120 my $rccmd = "$bowtie -F $win,$freq $btargs $pol $rcidx $fa";
 121 print STDERR "Forward command: $fwcmd\n";
 122 print STDERR "Reverse comp command: $rccmd\n";
 123 open BTFW, "$fwcmd |" || die "Couldn't open pipe '$fwcmd |'\n";
 124 open BTRC, "$rccmd |" || die "Couldn't open pipe '$rccmd |'\n";
 125
 126 # 10_554  4       *       0       0       *       *       0       0       NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN      IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII      XM:i:0
 127 # 10_555  4       *       0       0       *       *       0       0       NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN      IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII      XM:i:0
 128 # 10_556  4       *       0       0       *       *       0       0       NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN      IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII      XM:i:0
 129 # 10_557  4       *       0       0       *       *       0       0       NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN      IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII      XM:i:0
 130 # 10_558  4       *       0       0       *       *       0       0       NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN      IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII      XM:i:0
 131 # 10_559  4       *       0       0       *       *       0       0       NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN      IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII      XM:i:0
 132
 133 print STDERR "Reading...\n";
 134 my $ln = 0;
 135 my $cur = 0;
 136 my $lastc = "\n";
 137 while(1) {
 138         my $fwl = <BTFW>;
 139         my $rcl = <BTRC>;
 140
 141         last unless defined($fwl) && defined($rcl);
 142         $ln++;
 143
 144         my @fws = split(/\t/, $fwl);
 145         my @fws1 = split(/_/, $fws[0]);
 146         my @rcs = split(/\t/, $rcl);
 147         my @rcs1 = split(/_/, $rcs[0]);
 148
 149         my $cname = join("_", @fws1[0..$#fws1-1]);
 150         $cname eq join("_", @fws1[0..$#fws1-1]) || die "Name mismatch on line $ln:\n$fwl\n$rcl\n";
 151         my $off = $fws1[-1];
 152         $off eq $rcs1[-1] || die "Offset mismatch on line $ln:\n$fwl\n$rcl\n";
 153         # If a read aligns, XM:i is not printed
 154         # If a read fails to align, XM:i:0 is printed
 155         # If a read aligns multiple places, XM:i:N is printed where N>0
 156         my $fwxm0 = $fws[-1] =~ /XM:i:0/;
 157         my $rcxm0 = $rcs[-1] =~ /XM:i:0/;
 158         my $fwxm  = $fws[-1] =~ /XM:i:/;
 159         my $rcxm  = $rcs[-1] =~ /XM:i:/;
 160         # For mapable to be true, we need the read to have failed to align
 161         # in one index and aligned uniquely in the other
 162         my $mapable = (($fwxm0 && !$rcxm) || ($rcxm0 && !$fwxm)) ? 1 : 0;
 163
 164         $cname =~ s/\s.*//; # trim everything after first whitespace to get short name
 165         if($cname =~ /^FW:/ || $cname =~ /^RC:/) {
 166                 $cname = substr($cname, 3);
 167         }
 168         if($name ne $cname) {
 169                 if($name ne "") {
 170                         # Flush remaining characters from previous name
 171                         defined($lens{$name}) || die;
 172                         $cur == $lens{$name} - $win + 1 || die "name is $name, cur is $cur, len is $lens{$name}, win is $win\n";
 173                         for(; $cur < $lens{$name}; $cur++) {
 174                                 $running -= $last[$cur % $win];
 175                                 $last[$cur % $win] = 0;
 176                                 $lastc = chr($running + 64);
 177                                 print $lastc;
 178                                 if((($cur+1) % 60) == 0) {
 179                                         $lastc = "\n";
 180                                         print $lastc;
 181                                 }
 182                         }
 183                 }
 184                 $name = $cname;
 185                 defined($lens{$name}) || die "No such name as \"$name\"\n";
 186                 print "\n" unless $lastc eq "\n";
 187                 print ">$name\n";
 188                 $lastc = "\n";
 189                 $cur = 0;
 190                 $running = 0;
 191                 clearLast();
 192         }
 193
 194         $running -= $last[$cur % $win];
 195         $last[$cur % $win] = $mapable;
 196         $running += $last[$cur % $win];
 197
 198         $running <= $win || die "running counter $running got higher than window size $win\n";
 199         $lastc = chr($running + 64);
 200         print $lastc;
 201         if((($cur+1) % 60) == 0) {
 202                 $lastc = "\n";
 203                 print $lastc;
 204         }
 205
 206         $cur++;
 207 }
 208 defined($lens{$name}) || die;
 209 for(; $cur < $lens{$name}; $cur++) {
 210         $running -= $last[$cur % $win];
 211         $last[$cur % $win] = 0;
 212         $lastc = chr($running + 64);
 213         print $lastc;
 214         if((($cur+1) % 60) == 0) {
 215                 $lastc = "\n";
 216                 print $lastc;
 217         }
 218 }
 219 close(BTFW);
 220 $? == 0 || die "Bad exitlevel from forward bowtie: $?\n";
 221 close(BTRC);
 222 $? == 0 || die "Bad exitlevel from reverse-comp bowtie: $?\n";