]> git.donarmstrong.com Git - samtools.git/blob - misc/novo2sam.pl
* samtools-0.1.5-24 (r435)
[samtools.git] / 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-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-Z]+)\>/) {
176                         $indeltype=">";
177                         $insert=$1;
178                 } elsif ($string =~ /\-/) {
179                         $indeltype="-";
180                 } elsif ($string =~ /\+([A-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-Z]+)/) {  
208                         $indeltype="+";
209                         $insert = $1;
210                 }elsif ($string =~ /\-([A-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 }