misc/export2sam.pl

   1 #!/usr/bin/perl -w
   2
   3 # Contact: lh3
   4 # Version: 0.1.2 (03JAN2009)
   5
   6 use strict;
   7 use warnings;
   8 use Getopt::Std;
   9
  10 &export2sam;
  11 exit;
  12
  13 sub export2sam {
  14   my ($fh1, $fh2, $is_paired);
  15   $is_paired = (@ARGV >= 2);
  16   die("export2sam.pl <read1.export> [<read2.export>]\n") if (@ARGV == 0);
  17   open($fh1, $ARGV[0]) || die;
  18   if ($is_paired) {
  19         open($fh2, $ARGV[1]) || die;
  20   }
  21   # conversion table
  22   my @conv_table;
  23   for (-64..64) {
  24         $conv_table[$_+64] = chr(int(33 + 10*log(1+10**($_/10.0))/log(10)+.499));
  25   }
  26   # core loop
  27   while (<$fh1>) {
  28         my (@s1, @s2);
  29         &export2sam_aux($_, \@s1, \@conv_table, $is_paired);
  30         if ($is_paired) {
  31           $_ = <$fh2>;
  32           &export2sam_aux($_, \@s2, \@conv_table, $is_paired);
  33           if (@s1 && @s2) { # then set mate coordinate
  34                 my $isize = 0;
  35                 if ($s1[2] ne '*' && $s1[2] eq $s2[2]) { # then calculate $isize
  36                   my $x1 = ($s1[1] & 0x10)? $s1[3] + length($s1[9]) : $s1[3];
  37                   my $x2 = ($s2[1] & 0x10)? $s2[3] + length($s2[9]) : $s2[3];
  38                   $isize = $x2 - $x1;
  39                 }
  40                 # update mate coordinate
  41                 if ($s2[2] ne '*') {
  42                   @s1[6..8] = (($s2[2] eq $s1[2])? "=" : $s2[2], $s2[3], $isize);
  43                   $s1[1] |= 0x20 if ($s2[1] & 0x10);
  44                 } else {
  45                   $s1[1] |= 0x8;
  46                 }
  47                 if ($s1[2] ne '*') {
  48                   @s2[6..8] = (($s1[2] eq $s2[2])? "=" : $s1[2], $s1[3], -$isize);
  49                   $s2[1] |= 0x20 if ($s1[1] & 0x10);
  50                 } else {
  51                   $s2[1] |= 0x8;
  52                 }
  53           }
  54         }
  55         print join("\t", @s1), "\n" if (@s1);
  56         print join("\t", @s2), "\n" if (@s2 && $is_paired);
  57   }
  58   close($fh1);
  59   close($fh2) if ($is_paired);
  60 }
  61
  62 sub export2sam_aux {
  63   my ($line, $s, $ct, $is_paired) = @_;
  64   chomp($line);
  65   my @t = split("\t", $line);
  66   @$s = ();
  67   return if ($t[21] ne 'Y');
  68   # read name
  69   $s->[0] = $t[1]? "$t[0]_$t[1]:$t[2]:$t[3]:$t[4]:$t[5]" : "$t[0]:$t[2]:$t[3]:$t[4]:$t[5]";
  70   # initial flag (will be updated later)
  71   $s->[1] = 0;
  72   $s->[1] |= 1 | 1<<(5 + $t[7]) if ($is_paired);
  73   # read & quality
  74   $s->[9] = $t[8]; $s->[10] = $t[9];
  75   if ($t[13] eq 'R') { # then reverse the sequence and quality
  76         $s->[9] = reverse($t[8]);
  77         $s->[9] =~ tr/ACGTacgt/TGCAtgca/;
  78         $s->[10] = reverse($t[9]);
  79   }
  80   $s->[10] =~ s/(.)/$ct->[ord($1)]/eg; # change coding
  81   # cigar
  82   $s->[5] = length($s->[9]) . "M";
  83   # coor
  84   my $has_coor = 0;
  85   $s->[2] = "*";
  86   if ($t[10] eq 'NM' || $t[10] eq 'QC') {
  87         $s->[1] |= 0x4; # unmapped
  88   } elsif ($t[10] =~ /(\d+):(\d+):(\d+)/) {
  89         $s->[1] |= 0x4; # TODO: should I set BAM_FUNMAP in this case?
  90         push(@$s, "H0:i:$1", "H1:i:$2", "H2:i:$3")
  91   } else {
  92         $s->[2] = $t[10];
  93         $has_coor = 1;
  94   }
  95   $s->[3] = $has_coor? $t[12] : 0;
  96   $s->[1] |= 0x10 if ($has_coor && $t[13] eq 'R');
  97   # mapQ (TODO: should I choose the larger between $t[15] and $t[16]?)
  98   $s->[4] = 0;
  99   $s->[4] = $t[15] if ($t[15] ne '');
 100   $s->[4] = $t[16] if ($t[16] ne '' && $s->[4] < $t[16]);
 101   # mate coordinate
 102   $s->[6] = '*'; $s->[7] = $s->[8] = 0;
 103   # aux
 104   push(@$s, "BC:Z:$t[6]") if ($t[6]);
 105   push(@$s, "MD:Z:$t[14]") if ($has_coor);
 106   push(@$s, "SM:i:$t[15]") if ($is_paired && $has_coor);
 107 }