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];
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);
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);
39 die("Usage: soap2sam.pl [-p] <aln.soap>\n") if (@ARGV == 0);
40 my $is_paired = defined($opts{p});
44 my ($s_last, $s_curr) = (\@s1, \@s2);
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 = ();
53 print join("\t", @$s_last), "\n" if (@$s_last != 0);
54 my $s = $s_last; $s_last = $s_curr; $s_curr = $s;
57 print join("\t", @$s_last), "\n" if (@$s_last != 0);
61 my ($line, $s, $is_paired) = @_;
63 my @t = split("\t", $line);
64 return -1 if (@t < 9);
68 $s->[0] =~ s/\/[12]$//g;
69 # initial flag (will be updated later)
71 $s->[1] |= 1 | 1<<($t[4] eq 'a'? 6 : 7);
72 $s->[1] |= 2 if ($is_paired);
74 $s->[9] = $t[1]; $s->[10] = $t[2];
76 $s->[5] = length($s->[9]) . "M";
78 $s->[2] = $t[7]; $s->[3] = $t[8];
79 $s->[1] |= 0x10 if ($t[6] eq '-');
81 $s->[4] = $t[3] == 1? 30 : 0;
83 $s->[6] = '*'; $s->[7] = $s->[8] = 0;
85 push(@$s, "NM:i:$t[9]");
90 push(@x, sprintf("%.3d,$1", $2)) if ($t[$_] =~ /^([ACGT])->(\d+)/i);
95 my ($y, $z) = split(",");
96 $md .= (int($y)-$a) . $z;
99 $md .= length($t[1]) - $a;
103 push(@$s, "MD:Z:$md");