scripts/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 $idx = "";
  21 my $btargs = "-t -S --sam-nohead -M 1 --mm";
  22 my $debug = 0;
  23
  24 if(defined($ENV{BOWTIE_HOME})) {
  25         $bowtie = "$ENV{BOWTIE_HOME}/bowtie";
  26         unless(-x $bowtie) { $bowtie = "" };
  27 }
  28 if($bowtie eq "") {
  29         $bowtie = `which bowtie 2>/dev/null`;
  30         chomp($bowtie);
  31         unless(-x $bowtie) { $bowtie = "" };
  32 }
  33 $bowtie = "./bowtie" if ($bowtie eq "" && -x "./bowtie");
  34
  35 GetOptions(
  36         "fasta=s"     => \$fa,
  37         "window=i"    => \$win,
  38         "frequency=i" => \$freq,
  39         "bowtie=s"    => \$bowtie_arg,
  40         "policy=s"    => \$pol,
  41         "idx=s"       => \$idx,
  42         "debug"       => \$debug
  43 ) || die;
  44
  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";
  51
  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";
  55
  56 $bowtie = $bowtie_arg if $bowtie_arg ne "";
  57 unless(-x $bowtie) {
  58         # No bowtie? die
  59         if($bowtie_arg ne "") {
  60                 die "Specified -bowtie, \"$bowtie\" doesn't exist or isn't executable\n";
  61         } else {
  62                 die "bowtie couldn't be found in BOWTIE_HOME, PATH, or current directory; please specify -bowtie\n";
  63         }
  64 }
  65
  66 my $running = 0;
  67 my $name = ""; # name of sequence currently being processed
  68 my %lens = ();
  69 my @names = ();
  70 my $totlen = 0;
  71
  72 ##
  73 # Read lengths of all the entries in all the input fasta files.
  74 #
  75 sub readLens($) {
  76         my $ins = shift;
  77         my @is = split(/[,]/, $ins);
  78         for my $i (@is) {
  79                 open IN, "$i" || die "Could not open $i\n";
  80                 my $name = "";
  81                 while(<IN>) {
  82                         chomp;
  83                         if(substr($_, 0, 1) eq '>') {
  84                                 next if /\?[0-9]*$/; # Skip >?50000 lines
  85                                 $name = $_;
  86                                 $name = substr($name, 1); # Chop off >
  87                                 if($name =~ /^FW:/ || $name =~ /^RC:/) {
  88                                         $name = substr($name, 3); # Chop off FW:/RC:
  89                                 }
  90                                 my @ns = split(/\s+/, $name);
  91                                 $name = $ns[0]; # Get short name
  92                                 push @names, $name;
  93                                 print STDERR "Saw name $name\n";
  94                         } else {
  95                                 $name ne "" || die;
  96                                 $lens{$name} += length($_); # Update length
  97                                 $totlen += length($_);
  98                         }
  99                 }
 100                 close(IN);
 101         }
 102 }
 103 print STDERR "Reading fasta lengths\n";
 104 readLens($fa);
 105 print STDERR "  read ".scalar(keys %lens)." fasta sequences with total length $totlen\n";
 106
 107 my @last;
 108 for(my $i = 0; $i < $win; $i++) { push @last, 0 };
 109 sub clearLast {
 110         for(my $i = 0; $i < $win; $i++) { $last[$i] = 0 };
 111 }
 112
 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";
 117
 118 print STDERR "Reading...\n";
 119 my $ln = 0;
 120 my $cur = 0;
 121 my $lastc = "\n";
 122 while(<BT>) {
 123         $ln++;
 124
 125         my @s = split(/\t/, $_);
 126         my @s1 = split(/_/, $s[0]);
 127
 128         my $cname = join("_", @s1[0..$#s1-1]);
 129         my $off = $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);
 134
 135         $cname =~ s/\s.*//; # trim everything after first whitespace to get short name
 136         if($cname =~ /^FW:/ || $cname =~ /^RC:/) {
 137                 $cname = substr($cname, 3);
 138         }
 139         if($name ne $cname) {
 140                 if($name ne "") {
 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);
 148                                 print $lastc;
 149                                 if((($cur+1) % 60) == 0) {
 150                                         $lastc = "\n";
 151                                         print $lastc;
 152                                 }
 153                         }
 154                 }
 155                 $name = $cname;
 156                 defined($lens{$name}) || die "No such name as \"$name\"\n";
 157                 print "\n" unless $lastc eq "\n";
 158                 print ">$name\n";
 159                 $lastc = "\n";
 160                 $cur = 0;
 161                 $running = 0;
 162                 clearLast();
 163         }
 164
 165         $running -= $last[$cur % $win];
 166         $last[$cur % $win] = $mapable;
 167         $running += $last[$cur % $win];
 168
 169         $running <= $win || die "running counter $running got higher than window size $win\n";
 170         $lastc = chr($running + 64);
 171         print $lastc;
 172         if((($cur+1) % 60) == 0) {
 173                 $lastc = "\n";
 174                 print $lastc;
 175         }
 176
 177         $cur++;
 178 }
 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);
 184         print $lastc;
 185         if((($cur+1) % 60) == 0) {
 186                 $lastc = "\n";
 187                 print $lastc;
 188         }
 189 }
 190 close(BT);
 191 $? == 0 || die "Bad exitlevel from forward bowtie: $?\n";