#!/usr/bin/perl -w

($samFN, $alnFN) = @ARGV;

open $samFH, "$samFN" or die $!;
%sam = ();
while(<$samFH>)
{
	next if /^@/;
	@t = split;
	if(not defined $sam{$t[0]})
	{
		$sam{$t[0]} = "$t[2]-$t[3]";
	}
}
close $samFH;

open $alnFH, "$alnFN" or die $!;
%roc = ();
$maxScore = 0;
$unaligned = 0;
$aligned = 0;
$correct = 0;
$incorrect = 0;
while(<$alnFH>)
{
	next if /^@/;
	$l = $_;
	@t = split;
	$t[0] =~ s/(\S+?)-\S+/$1/ if $t[0] =~ /-/;
	$t[0] =~ s/^(\S+)\/[12]$/$1/ if $t[0] =~ /\/[12]$/;
	if(($t[1] & 4) != 0)
	{
		$roc{$t[4]}{0} = 0 if(not defined $roc{$t[4]}{0});
		$roc{$t[4]}{0} += 1;
		++$unaligned;
		next;
	}
	++$aligned;
	if(not defined $sam{$t[0]})
	{
		die "Read not found $t[0]\n";
	}
	$sam{$t[0]} =~ /(\S+)-(\S+)/;
	$c1 = $1; $d1 = $2;
	$ok = 0;
	$xa = "";
	$xa = $1 if $l =~ /XA:Z:(\S+)/;
	$l =~ /XG:i:(\d+)/;
	$xg = $1;
	if($c1 eq $t[2]) 
	{
		if(($d1-$t[3]) <= 50 and ($d1-$t[3]) >= -50)
		{
			++$ok;
			#print STDERR "$t[0]\t$t[2],$t[3],$t[5],$xg\t$t[2]\t$t[3]\t$t[5]\t$xa\n" if ($xa ne "");
		}
	}
	while($xa =~ /(chr\d+),[+-](\d+),(\S+?),(\d+);/g)
	{
		if($c1 eq $1)
		{
			if(($d1-$2) <= 50 and ($d1-$2) >= -50)
			{
				++$ok;
				#print STDERR "$t[0]\t$1,$2,$3,$4\t$t[2]\t$t[3]\t$t[5]\t$xa\n";
			}
		}
	}
	if($ok == 0)
	{
		$roc{$t[4]}{0} = 0 if(not defined $roc{$t[4]}{0});
		$roc{$t[4]}{0} += 1;
		++$incorrect;
	}
	else
	{
		$roc{$t[4]}{1} = 0 if(not defined $roc{$t[4]}{1});
		$roc{$t[4]}{1} += 1;
		++$correct;
	}
	$maxScore = $t[4] if($t[4] > $maxScore);
}
close $alnFH;

$incorr = $corr = 0;

for(;$maxScore >= 0; --$maxScore)
{
	if(not defined $roc{$maxScore})
	{
		print "$maxScore\t$incorr\t$corr\n";
	}
	else
	{
		print "$maxScore\t";
		if(not defined $roc{$maxScore}{0})
		{
			print "$incorr\t";
		}
		else
		{
			$incorr += $roc{$maxScore}{0};
			print "$incorr\t";
		}
		if(not defined $roc{$maxScore}{1})
		{
			print "$corr\n";
		}
		else
		{
			$corr += $roc{$maxScore}{1};
			print "$corr\n";
		}
	}
}

print "Unaligned\t$unaligned\tAligned\t$aligned\tIncorrect\t$incorrect\tCorrect\t$correct\n";

0;
