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: novo2sam.pl [-p] <aln.novo>\nWarning: gapped alignments are NOT converted properly!\n") if (@ARGV == 0);
40 my $is_paired = defined($opts{p});
44 my ($s_last, $s_curr) = (\@s1, \@s2);
47 next if (/(QC|NM)\s*$/ || /(R\s+\d+)\s*$/);
48 &novo2sam_aux($_, $s_curr, $is_paired);
49 if (@$s_last != 0 && $s_last->[0] eq $s_curr->[0]) {
50 &mating($s_last, $s_curr);
51 print join("\t", @$s_last), "\n";
52 print join("\t", @$s_curr), "\n";
53 @$s_last = (); @$s_curr = ();
55 print join("\t", @$s_last), "\n" if (@$s_last != 0);
56 my $s = $s_last; $s_last = $s_curr; $s_curr = $s;
59 print join("\t", @$s_last), "\n" if (@$s_last != 0);
63 my ($line, $s, $is_paired) = @_;
65 my @t = split(/\s+/, $line);
67 return if ($t[4] ne 'U');
68 my $len = length($t[2]);
70 $s->[0] = substr($t[0], 1);
71 $s->[0] =~ s/\/[12]$//g;
72 # initial flag (will be updated later)
74 $s->[1] |= 1 | 1<<($t[1] eq 'L'? 6 : 7);
75 $s->[1] |= 2 if ($t[10] eq '.');
78 $s->[9] = reverse($t[2]);
79 $s->[10] = reverse($t[3]);
80 $s->[9] =~ tr/ACGTRYMKWSNacgtrymkwsn/TGCAYRKMWSNtgcayrkmwsn/;
82 $s->[9] = $t[2]; $s->[10] = $t[3];
85 $s->[5] = $len . "M"; # IMPORTANT: this cigar is not correct for gapped alignment
87 $s->[2] = substr($t[7], 1); $s->[3] = $t[8];
88 $s->[1] |= 0x10 if ($t[9] eq 'R');
90 $s->[4] = $t[5] > $t[6]? $t[5] : $t[6];
92 $s->[6] = '*'; $s->[7] = $s->[8] = 0;
94 push(@$s, "NM:i:".(@t-13));
99 push(@x, sprintf("%.4d,$2", $1-1)) if ($t[$_] =~ /^(\d+)([ACGT])>([ACGT])/i);
104 my ($y, $z) = split(",");
105 $md .= (int($y)-$a? (int($y)-$a) : '') . $z;
108 $md .= $len - $a if ($len - $a);
112 push(@$s, "MD:Z:$md");