4 # Version: 0.1.1 (03JAN2009)
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;
19 open($fh2, $ARGV[1]) || die;
24 $conv_table[$_+64] = chr(int(33 + 10*log(1+10**($_/10.0))/log(10)+.499));
29 &export2sam_aux($_, \@s1, \@conv_table, $is_paired);
32 &export2sam_aux($_, \@s2, \@conv_table, $is_paired);
33 if (@s1 && @s2) { # then set mate coordinate
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];
40 # update mate coordinate
42 @s1[6..8] = (($s2[2] eq $s1[2])? "=" : $s2[2], $s2[3], $isize);
43 $s1[1] |= 0x20 if ($s2[1] & 0x10);
48 @s2[6..8] = (($s1[2] eq $s2[2])? "=" : $s1[2], $s1[3], -$isize);
49 $s2[1] |= 0x20 if ($s1[1] & 0x10);
55 print join("\t", @s1), "\n" if (@s1);
56 print join("\t", @s2), "\n" if (@s2 && $is_paired);
59 close($fh2) if ($is_paired);
63 my ($line, $s, $ct, $is_paired) = @_;
65 my @t = split("\t", $line);
67 return if ($t[21] ne 'Y');
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)
72 $s->[1] |= 1 | 1<<(5 + $t[7]) if ($is_paired);
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]);
80 $s->[10] =~ s/(.)/$ct->[ord($1)]/eg; # change coding
82 $s->[5] = length($s->[9]) . "M";
86 if ($t[10] eq 'NM' || $t[10] eq 'QC') {
87 $s->[1] |= 0x8; # unmapped
88 } elsif ($t[10] =~ /(\d+):(\d+):(\d+)/) {
89 $s->[1] |= 0x8; # TODO: should I set BAM_FUNMAP in this case?
90 push(@$s, "H0:i:$1", "H1:i:$2", "H2:i:$3")
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]?)
99 $s->[4] = $t[15] if ($t[15] ne '');
100 $s->[4] = $t[16] if ($t[16] ne '' && $s->[4] < $t[16]);
102 $s->[6] = '*'; $s->[7] = $s->[8] = 0;
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);