← Back to branch summary

~ubuntu-branches/ubuntu/quantal/genometools/quantal-backports

« back to all changes in this revision

Viewing changes to src/external/samtools-0.1.18/misc/bowtie2sam.pl

  • Committer: Package Import Robot
  • Author(s): Sascha Steinbiss
  • Date: 2012-07-09 14:10:23 UTC
  • Revision ID: package-import@ubuntu.com-20120709141023-juuu4spm6chqsf9o
Tags: upstream-1.4.1
Import upstream version 1.4.1

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/external/samtools-0.1.18/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
}