misc/bowtie2sam.pl

   1 #!/usr/bin/perl -w
   2
   3 # Contact: lh3
   4 # Version: 0.1.1
   5
   6 use strict;
   7 use warnings;
   8 use Getopt::Std;
   9
  10 &bowtie2sam;
  11 exit;
  12
  13 sub bowtie2sam {
  14   my %opts = ();
  15   die("Usage: bowtie2sam.pl <aln.bowtie>\n") if (@ARGV == 0 && -t STDIN);
  16   # core loop
  17   my (@s, $last, @staging, $k, $best_s, $subbest_s, $best_k);
  18   $last = '';
  19   while (<>) {
  20         my ($name, $nm) = &bowtie2sam_aux($_, \@s); # read_name, number of mismatches
  21         if ($name eq $last) {
  22           # I do not know whether the multiple hits are ordered on the
  23           # number of mismatches. I assume they are not and so I have to
  24           # keep all these multiple hits in memory.
  25           @{$staging[$k]} = @s;
  26           if ($best_s > $nm) {
  27                 $subbest_s = $best_s;
  28                 $best_s = $nm;
  29                 $best_k = $k;
  30           } elsif ($subbest_s > $nm) {
  31                 $subbest_s = $nm;
  32           }
  33           ++$k;
  34         } else {
  35           if ($last) {
  36                 if ($best_s == $subbest_s) {
  37                   $staging[$best_k][4] = 0;
  38                 } elsif ($subbest_s - $best_s == 1) {
  39                   $staging[$best_k][4] = 15 if ($staging[$best_k][4] > 15);
  40                 }
  41                 print join("\t", @{$staging[$best_k]}), "\n";
  42           }
  43           $k = 1; $best_s = $nm; $subbest_s = 1000; $best_k = 0;
  44           @{$staging[0]} = @s;
  45           $last = $name;
  46         }
  47   }
  48   print join("\t", @{$staging[$best_k]}), "\n" if ($best_k >= 0);
  49 }
  50
  51 sub bowtie2sam_aux {
  52   my ($line, $s) = @_;
  53   chomp($line);
  54   my @t = split("\t", $line);
  55   my $ret;
  56   @$s = ();
  57   # read name
  58   $s->[0] = $ret = $t[0];
  59   $s->[0] =~ s/\/[12]$//g;
  60   # initial flag (will be updated later)
  61   $s->[1] = 0;
  62   # read & quality
  63   $s->[9] = $t[4]; $s->[10] = $t[5];
  64   # cigar
  65   $s->[5] = length($s->[9]) . "M";
  66   # coor
  67   $s->[2] = $t[2]; $s->[3] = $t[3] + 1;
  68   $s->[1] |= 0x10 if ($t[1] eq '-');
  69   # mapQ
  70   $s->[4] = $t[6] == 0? 25 : 0;
  71   # mate coordinate
  72   $s->[6] = '*'; $s->[7] = $s->[8] = 0;
  73   # aux
  74   my $nm = @t - 7;
  75   push(@$s, "NM:i:" . (@t-7));
  76   push(@$s, "X$nm:i:" . ($t[6]+1));
  77   my $md = '';
  78   if ($t[7]) {
  79         $_ = $t[7];
  80         my $a = 0;
  81         while (/(\d+):[ACGTN]>([ACGTN])/gi) {
  82           my ($y, $z) = ($1, $2);
  83           $md .= (int($y)-$a) . $z;
  84           $a += $y - $a + 1;
  85         }
  86         $md .= length($s->[9]) - $a;
  87   } else {
  88         $md = length($s->[9]);
  89   }
  90   push(@$s, "MD:Z:$md");
  91   return ($ret, $nm);
  92 }