misc/soap2sam.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 &soap2sam;
  11 exit;
  12
  13 sub mating {
  14   my ($s1, $s2) = @_;
  15   my $isize = 0;
  16   if ($s1->[2] ne '*' && $s1->[2] eq $s2->[2]) { # then calculate $isize
  17         my $x1 = ($s1->[1] & 0x10)? $s1->[3] + length($s1->[9]) : $s1->[3];
  18         my $x2 = ($s2->[1] & 0x10)? $s2->[3] + length($s2->[9]) : $s2->[3];
  19         $isize = $x2 - $x1;
  20   }
  21   # update mate coordinate
  22   if ($s2->[2] ne '*') {
  23         @$s1[6..8] = (($s2->[2] eq $s1->[2])? "=" : $s2->[2], $s2->[3], $isize);
  24         $s1->[1] |= 0x20 if ($s2->[1] & 0x10);
  25   } else {
  26         $s1->[1] |= 0x8;
  27   }
  28   if ($s1->[2] ne '*') {
  29         @$s2[6..8] = (($s1->[2] eq $s2->[2])? "=" : $s1->[2], $s1->[3], -$isize);
  30         $s2->[1] |= 0x20 if ($s1->[1] & 0x10);
  31   } else {
  32         $s2->[1] |= 0x8;
  33   }
  34 }
  35
  36 sub soap2sam {
  37   my %opts = ();
  38   getopts("p", \%opts);
  39   die("Usage: soap2sam.pl [-p] <aln.soap>\n") if (@ARGV == 0);
  40   my $is_paired = defined($opts{p});
  41   # core loop
  42   my @s1 = ();
  43   my @s2 = ();
  44   my ($s_last, $s_curr) = (\@s1, \@s2);
  45   while (<>) {
  46         next if (&soap2sam_aux($_, $s_curr, $is_paired) < 0);
  47         if (@$s_last != 0 && $s_last->[0] eq $s_curr->[0]) {
  48           &mating($s_last, $s_curr);
  49           print join("\t", @$s_last), "\n";
  50           print join("\t", @$s_curr), "\n";
  51           @$s_last = (); @$s_curr = ();
  52         } else {
  53           print join("\t", @$s_last), "\n" if (@$s_last != 0);
  54           my $s = $s_last; $s_last = $s_curr; $s_curr = $s;
  55         }
  56   }
  57   print join("\t", @$s_last), "\n" if (@$s_last != 0);
  58 }
  59
  60 sub soap2sam_aux {
  61   my ($line, $s, $is_paired) = @_;
  62   chomp($line);
  63   my @t = split("\t", $line);
  64   return -1 if (@t < 9);
  65   @$s = ();
  66   # read name
  67   $s->[0] = $t[0];
  68   $s->[0] =~ s/\/[12]$//g;
  69   # initial flag (will be updated later)
  70   $s->[1] = 0;
  71   $s->[1] |= 1 | 1<<($t[4] eq 'a'? 6 : 7);
  72   $s->[1] |= 2 if ($is_paired);
  73   # read & quality
  74   $s->[9] = $t[1]; $s->[10] = $t[2];
  75   # cigar
  76   $s->[5] = length($s->[9]) . "M";
  77   # coor
  78   $s->[2] = $t[7]; $s->[3] = $t[8];
  79   $s->[1] |= 0x10 if ($t[6] eq '-');
  80   # mapQ
  81   $s->[4] = $t[3] == 1? 30 : 0;
  82   # mate coordinate
  83   $s->[6] = '*'; $s->[7] = $s->[8] = 0;
  84   # aux
  85   push(@$s, "NM:i:$t[9]");
  86   my $md = '';
  87   if ($t[9]) {
  88         my @x;
  89         for (10 .. $#t) {
  90           push(@x, sprintf("%.3d,$1", $2)) if ($t[$_] =~ /^([ACGT])->(\d+)/i);
  91         }
  92         @x = sort(@x);
  93         my $a = 0;
  94         for (@x) {
  95           my ($y, $z) = split(",");
  96           $md .= (int($y)-$a) . $z;
  97           $a += $y - $a + 1;
  98         }
  99         $md .= length($t[1]) - $a;
 100   } else {
 101         $md = length($t[1]);
 102   }
 103   push(@$s, "MD:Z:$md");
 104   return 0;
 105 }