+ my $score = -1;
+ print $_ if (/^\@/);
+ $score = $1 if (/AS:i:(\d+)/);
+ my @t = split("\t");
+ next if (@t < 11);
+ if ($score < 0) { # AS tag is unavailable
+ my $cigar = $t[5];
+ my ($mm, $go, $ge) = (0, 0, 0);
+ $cigar =~ s/(\d+)[ID]/++$go,$ge+=$1/eg;
+ $cigar = $t[5];
+ $cigar =~ s/(\d+)M/$mm+=$1/eg;
+ $score = $mm * $opts{a} - $go * $opts{q} - $ge * $opts{r}; # no mismatches...
+ }
+ $score = 1 if ($score < 1);
+ if ($t[0] ne $last) {
+ &unique_aux(\@a, $opts{f}, $recal_Q, $multi_only) if (@a);
+ $last = $t[0];
+ }
+ push(@a, [$score, \@t]);
+ }
+ &unique_aux(\@a, $opts{f}, $recal_Q, $multi_only) if (@a);
+}
+
+sub unique_aux {
+ my ($a, $fac, $is_recal, $multi_only) = @_;
+ my ($max, $max2, $max_i) = (0, 0, -1);
+ for (my $i = 0; $i < @$a; ++$i) {
+ if ($a->[$i][0] > $max) {
+ $max2 = $max; $max = $a->[$i][0]; $max_i = $i;
+ } elsif ($a->[$i][0] > $max2) {
+ $max2 = $a->[$i][0];
+ }
+ }
+ if ($is_recal) {
+ if (!$multi_only || @$a > 1) {
+ my $q = int($fac * ($max - $max2) / $max + .499);
+ $q = 250 if ($q > 250);
+ $a->[$max_i][1][4] = $q < 250? $q : 250;
+ }
+ }
+ print join("\t", @{$a->[$max_i][1]});
+ @$a = ();
+}
+
+#
+# uniqcmp: compare two SAM files
+#
+
+sub uniqcmp {
+ my %opts = (q=>10, s=>100);
+ getopts('pq:s:', \%opts);
+ die("Usage: samtools.pl uniqcmp <in1.sam> <in2.sam>\n") if (@ARGV < 2);
+ my ($fh, %a);
+ warn("[uniqcmp] read the first file...\n");
+ &uniqcmp_aux($ARGV[0], \%a, 0);
+ warn("[uniqcmp] read the second file...\n");
+ &uniqcmp_aux($ARGV[1], \%a, 1);
+ warn("[uniqcmp] stats...\n");
+ my @cnt;
+ $cnt[$_] = 0 for (0..9);
+ for my $x (keys %a) {
+ my $p = $a{$x};
+ my $z;
+ if (defined($p->[0]) && defined($p->[1])) {
+ $z = ($p->[0][0] == $p->[1][0] && $p->[0][1] eq $p->[1][1] && abs($p->[0][2] - $p->[1][2]) < $opts{s})? 0 : 1;
+ if ($p->[0][3] >= $opts{q} && $p->[1][3] >= $opts{q}) {
+ ++$cnt[$z*3+0];
+ } elsif ($p->[0][3] >= $opts{q}) {
+ ++$cnt[$z*3+1];
+ } elsif ($p->[1][3] >= $opts{q}) {
+ ++$cnt[$z*3+2];
+ }
+ print STDERR "$x\t$p->[0][1]:$p->[0][2]\t$p->[0][3]\t$p->[0][4]\t$p->[1][1]:$p->[1][2]\t$p->[1][3]\t$p->[1][4]\t",
+ $p->[0][5]-$p->[1][5], "\n" if ($z && defined($opts{p}) && ($p->[0][3] >= $opts{q} || $p->[1][3] >= $opts{q}));
+ } elsif (defined($p->[0])) {
+ ++$cnt[$p->[0][3]>=$opts{q}? 6 : 7];
+ print STDERR "$x\t$p->[0][1]:$p->[0][2]\t$p->[0][3]\t$p->[0][4]\t*\t0\t*\t",
+ $p->[0][5], "\n" if (defined($opts{p}) && $p->[0][3] >= $opts{q});