0
|
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 }
|