misc/novo2sam.pl

   1 #!/usr/bin/perl -w
   2
   3 # Contact: lh3
   4 # Version: 0.1.3
   5
   6 #Modified by Zayed Albertyn(zayed.albertyn@gmail.com) & Colin Hercus(colin@novocraft.com)
   7
   8 #use strict;
   9 #use warnings;
  10 use Data::Dumper;
  11 use Getopt::Std;
  12
  13 &novo2sam;
  14 exit;
  15
  16 sub mating {
  17   my ($s1, $s2) = @_;
  18   my $isize = 0;
  19   if ($s1->[2] ne '*' && $s1->[2] eq $s2->[2]) { # then calculate $isize
  20         my $x1 = ($s1->[1] & 0x10)? $s1->[3] + length($s1->[9]) : $s1->[3];
  21         my $x2 = ($s2->[1] & 0x10)? $s2->[3] + length($s2->[9]) : $s2->[3];
  22         $isize = $x2 - $x1;
  23   }
  24   # update mate coordinate
  25   if ($s2->[2] ne '*') {
  26         @$s1[6..8] = (($s2->[2] eq $s1->[2])? "=" : $s2->[2], $s2->[3], $isize);
  27         $s1->[1] |= 0x20 if ($s2->[1] & 0x10);
  28   } else {
  29         $s1->[1] |= 0x8;
  30   }
  31   if ($s1->[2] ne '*') {
  32         @$s2[6..8] = (($s1->[2] eq $s2->[2])? "=" : $s1->[2], $s1->[3], -$isize);
  33         $s2->[1] |= 0x20 if ($s1->[1] & 0x10);
  34   } else {
  35         $s2->[1] |= 0x8;
  36   }
  37 }
  38
  39 sub novo2sam {
  40   my %opts = ();
  41   getopts("p", \%opts);
  42   die("Usage: novo2sam.pl [-p] <aln.novo>\n") if (@ARGV == 0);
  43   my $is_paired = defined($opts{p});
  44   # core loop
  45   my @s1 = ();
  46   my @s2 = ();
  47   my ($s_last, $s_curr) = (\@s1, \@s2);
  48   while (<>) {
  49         next if (/^#/);
  50         next if (/(QC|NM)\s*$/ || /(R\s+\d+)\s*$/);
  51         &novo2sam_aux($_, $s_curr, $is_paired);
  52         if (@$s_last != 0 && $s_last->[0] eq $s_curr->[0]) {
  53           &mating($s_last, $s_curr);
  54           print join("\t", @$s_last), "\n";
  55           print join("\t", @$s_curr), "\n";
  56           @$s_last = (); @$s_curr = ();
  57         } else {
  58           print join("\t", @$s_last), "\n" if (@$s_last != 0);
  59           my $s = $s_last; $s_last = $s_curr; $s_curr = $s;
  60         }
  61   }
  62   print join("\t", @$s_last), "\n" if (@$s_last != 0);
  63 }
  64
  65 sub novo2sam_aux {
  66   my ($line, $s, $is_paired) = @_;
  67
  68   chomp($line);
  69   my @t = split(/\s+/, $line);
  70   my @variations =  @t[13 .. $#t];
  71   @$s = ();
  72   return if ($t[4] ne 'U');
  73   my $len = length($t[2]);
  74   # read name
  75   $s->[0] = substr($t[0], 1);
  76   $s->[0] =~ s/\/[12]$//g;
  77   # initial flag (will be updated later)
  78   $s->[1] = 0;
  79   $s->[1] |= 1 | 1<<($t[1] eq 'L'? 6 : 7);
  80   $s->[1] |= 2 if ($t[10] eq '.');
  81   # read & quality
  82   if ($t[9] eq 'R') {
  83         $s->[9] = reverse($t[2]);
  84         $s->[10] = reverse($t[3]);
  85         $s->[9] =~ tr/ACGTRYMKWSNacgtrymkwsn/TGCAYRKMWSNtgcayrkmwsn/;
  86   } else {
  87         $s->[9] = $t[2]; $s->[10] = $t[3];
  88   }
  89   # cigar
  90    my $cigarstring ="";
  91   if (scalar @variations ==0 ) {
  92          $s->[5] = $len . "M"; # IMPORTANT: this cigar is not correct for gapped alignment
  93   } else {
  94         #convert to correct CIGAR
  95         my $tmpstr =  join" ",@variations ;
  96         if ( $tmpstr=~ /\+|\-/ ) {
  97                 $cigarstring  = cigar_method($line,\@variations,$len);
  98                 $s->[5]=$cigarstring;
  99         } else {
 100                 $s->[5]=$len. "M";
 101         }
 102 }
 103
 104 # coor
 105   $s->[2] = substr($t[7], 1); $s->[3] = $t[8];
 106   $s->[1] |= 0x10 if ($t[9] eq 'R');
 107   # mapQ
 108   $s->[4] = $t[5] > $t[6]? $t[5] : $t[6];
 109   # mate coordinate
 110   $s->[6] = '*'; $s->[7] = $s->[8] = 0;
 111   # aux
 112   push(@$s, "NM:i:".(@t-13));
 113   my $md = '';
 114   $md = mdtag($md,$line,\@variations,$len);
 115   push(@$s, "MD:Z:$md");
 116
 117 }
 118
 119 sub mdtag {
 120         my $oldmd = shift;
 121         my $line = shift;
 122         my $ref =shift;
 123         my $rdlen  = shift;
 124         my @variations = @$ref;
 125         my $string="";
 126         my $mdtag="";
 127         my $t=1;
 128         my $q=1;
 129         my $deleteflag=0;
 130         my $len =0;
 131         foreach $string (@variations) {
 132                 my ($indeltype,$insert) = indeltype($string);
 133                 if ($indeltype eq "+") {
 134                         $len = length ($insert);
 135                         $q+=$len;
 136                   next;
 137                 }
 138                 my $pos = $1 if $string =~ /^(\d+)/;
 139                 $len = $pos - $t;
 140                 if ($len !=0 || ($deleteflag eq 1 && $indeltype eq ">")) {
 141                         $mdtag.=$len;
 142                 }
 143                 $t+=$len;
 144                 $q+=$len;
 145                 if ($indeltype eq ">") {
 146                         $mdtag.=$insert;
 147                         $deleteflag=0;
 148                       $t+=1;
 149                       $q+=1;
 150                 }
 151                 if ($indeltype eq "-") {
 152                         my $deletedbase = $2 if $string =~ /(\d+)\-([A-Za-z]+)/;
 153                         if ($deleteflag == 0 ) {
 154                                 $mdtag.="^";
 155                         }
 156                         $mdtag.=$deletedbase;
 157                         $deleteflag=1;
 158                         $t+=1;
 159                 }
 160         }
 161         $len = $rdlen - $q + 1;
 162         if ($len > 0) {
 163                 $mdtag.="$len";
 164         }
 165 #       print "In:$line\n";
 166 #       print "MD: OLD => NEW\nMD: $oldmd => $mdtag\n\n";
 167
 168         return $mdtag;
 169 }
 170
 171 sub indeltype {
 172         my $string =  shift;
 173         my $insert="";
 174         my $indeltype;
 175                if ($string =~ /([A-Za-z]+)\>/) {
 176                         $indeltype=">";
 177                         $insert=$1;
 178                 } elsif ($string =~ /\-/) {
 179                         $indeltype="-";
 180                 } elsif ($string =~ /\+([A-Za-z]+)/) {
 181                         $indeltype="+";
 182                         $insert=$1;
 183                 }
 184          return ($indeltype,$insert);
 185
 186 }
 187
 188
 189 sub cigar_method {
 190         my $line = shift;
 191         my $ref =shift;
 192         my $rdlen  = shift;
 193         my @variations = @$ref;
 194         my $string="";
 195         my $type="";
 196         my $t =1;
 197         my $q=1;
 198         my $indeltype="";
 199         my $cigar=  "";
 200         my $insert = "";
 201         my $len=0;
 202         my @cig=();
 203         foreach $string (@variations) {
 204                 next if $string =~  />/;
 205                 my $pos = $1 if $string =~ /^(\d+)/;
 206
 207                 if ($string =~ /\+([A-Za-z]+)/) {
 208                         $indeltype="+";
 209                         $insert = $1;
 210                 }elsif ($string =~ /\-([A-Za-z]+)/) {
 211                         $indeltype="-";
 212                         $insert = $1;
 213                 }
 214 #print "$pos $indeltype $insert $t $q\n";
 215                 $len = $pos - $t;
 216                 if ( $len > 0) {
 217                         $cigar.=$len."M";
 218                         push(@cig,$len."M");
 219                 }
 220                 $t+=$len;
 221                 $q+=$len;
 222
 223                 if ($indeltype eq "-") {
 224                         $cigar.="D";
 225                         push(@cig,"D");
 226                         $t++;
 227                 }
 228                 if ($indeltype eq "+") {
 229                         $len = length ($insert);
 230                         if ($len == 1) {
 231                                 $cigar.="I";
 232                                 push(@cig,"I");
 233                         }
 234                         if ($len > 1) {
 235                                 $cigar.=$len."I";
 236                                 push(@cig,$len."I")
 237                         }
 238                         $q+=$len;
 239                 }
 240                 $insert="";
 241         }
 242         $len= $rdlen - $q + 1;
 243         if ($len > 0) {
 244                 $cigar.=$len."M";
 245                 push(@cig,$len."M");
 246         }
 247
 248       $cigar = newcigar($cigar,'D');
 249       $cigar = newcigar($cigar,'I');
 250
 251         #print "$line\n";
 252         #print "c CIGAR:\t$cigar\n\n";
 253         return $cigar;
 254
 255 }
 256
 257
 258
 259 sub newcigar {
 260         my $cigar = shift;
 261         my $char = shift;
 262         my $new = "";
 263         my $copy = $cigar;
 264 #print "$cigar\n";
 265         $copy =~ s/^($char+)/$1;/g;
 266 #print "$copy\n";
 267         $copy =~ s/([^0-9$char])($char+)/$1;$2;/g;
 268 #print "$copy\n";
 269         my @parts = split(/;/,$copy);
 270         my $el="";
 271         foreach $el (@parts) {
 272 #print "$el\n";
 273                 if ($el =~ /^$char+$/) {
 274                         $new.=length($el).$char;
 275                 }else {
 276                         $new.=$el;
 277                 }
 278
 279         }
 280      return  $new;
 281 }