]> git.donarmstrong.com Git - samtools.git/blob - misc/bowtie2sam.pl
* samtools-0.1.5-24 (r435)
[samtools.git] / 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 }