#!/usr/bin/perl
use strict;
use Getopt::Std;

####################################################################
#
#  Evaluator.pl
#
#  Enrique Blanco for ABS v 1.0
#
#  INPUT: Length, annotations.gff, predictions.gff, PID
#
#  STDOUT: accuracy assesment at nucleotide and site levels
#  OUTPUT_FILE: ANNOT+PRED for gff2psimage
####################################################################

### CONSTANTS
my $TRUE=1;
my $FALSE=0;
my $OVERLAP=4;
my $GAME="EVALUATOR";
my $sANNOTATION="ANNOTATION";
my $sPREDICTION="PREDICTION";

### STEP 0. READING VALUES AND OPTIONS
my %opts;
my ($LENGTH,$annot_file,$pred_file,$PID);
my $OUTPUT_FILE;                                                               

getopts('vh',\%opts);
print_help();

if (scalar(@ARGV) != 4)
{
    print_error("Wrong Number of Parameters (4 is OK)");
}

$LENGTH = shift(@ARGV);
$annot_file = shift(@ARGV);
$pred_file = shift(@ARGV);
$PID = shift(@ARGV);
$OUTPUT_FILE ="/tmp/EVALUATOR_$PID.gff";


### STEP 1. LOADING ANNOTATIONS 
print_mess("Step 1. Loading Annotations");

my ($line,@record);
my %ACOVERAGE;
my %ANNOT;
my $i;
my ($locus,$pos1,$pos2);

open(SITES,$annot_file) or die (print_error("Can't open $annot_file"));
while($line=<SITES>)
{
    @record = split(/\s+/,$line);

    $locus = $record[0];
    $pos1 = $record[3];
    $pos2 = $record[4];

    # Initialize the ACOVERAGE
    if (!defined($ACOVERAGE{$locus}))
    {
	for($i=0; $i<$LENGTH; $i++)
	{
	    ${$ACOVERAGE{$locus}}[$i]=$FALSE;
        }
    }

    # Mark-up the positions of this motif
    # position corrections: -1
    for($i=$pos1-1; $i<$pos2-1+1; $i++)
    {
	${$ACOVERAGE{$locus}}[$i]=$TRUE;
    }

    # Save the original gff lines for posterior image
    push(@{$ANNOT{$locus}},$line);

}
close(SITES);


### STEP 2. LOADING PREDICTIONS
print_mess("Step 2. Loading Predictions");

my %COVERAGE;
my %PRED;

open(SITES,$pred_file) or die (print_error("Can't open $pred_file"));
while($line=<SITES>)
{
    @record = split(/\s+/,$line);

    $locus = $record[0];
    $pos1 = $record[3];
    $pos2 = $record[4];

    # Initialize the ACOVERAGE
    if (!defined($COVERAGE{$locus}))
    {
	for($i=0; $i<$LENGTH; $i++)
	{
	    ${$COVERAGE{$locus}}[$i]=$FALSE;
        }
    }

    # Mark-up the positions of this motif
    # position corrections: -1
    for($i=$pos1-1; $i<$pos2-1+1; $i++)
    {
	    ${$COVERAGE{$locus}}[$i]=$TRUE;
    }

    # Save the original gff lines for posterior image
    push(@{$PRED{$locus}},$line);

}
close(SITES);


### STEP 3. COMPUTING ACC at nucleotide level
print_mess("Step 3. Computing Accuracy at nucleotide level");

my ($TP,$FN,$FP,$TN,$CC);
my ($SN,$SP,$CC);
my ($allSN,$allSP,$allCC);
my $key;

print "<BR><BR><B><FONT COLOR=blue>Nucleotide Level Accuracy:</FONT></B><BR>";
print "<TABLE BORDER=1 BORDERCOLOR=red><TR><TD>";
print "<TABLE BORDER=0 CELLPADDING=10 CELLSPACING=10><TR><TH BGCOLOR=blue><FONT COLOR=white>LOCUS</TH><TH BGCOLOR=blue><FONT COLOR=white>TP</TH><TH BGCOLOR=blue><FONT COLOR=white>FP</TH><TH BGCOLOR=blue><FONT COLOR=white>FN</TH><TH BGCOLOR=blue><FONT COLOR=white>TN</TH><TH></TH><TH BGCOLOR=blue><FONT COLOR=white>SN</TH><TH BGCOLOR=blue><FONT COLOR=white>SP</TH><TH BGCOLOR=blue><FONT COLOR=white>CC</TH>";

$allSN = $allSP = $allCC = 0;
foreach $key (sort keys %ACOVERAGE)
{
    print_mess("Processing annotations from sequence $key");
    $TP = $FN = $FP = $TN = $CC = 0;

    # NO PREDICTIONS FOR THIS SEQUENCE
    if (!defined($COVERAGE{$key}))
    {
	$TP = $FP = 0;
	# all of the annotations are missed
	for($i=0; $i<$LENGTH; $i++)
	{
	    if (${$ACOVERAGE{$key}}[$i] == $TRUE)
	    {
		$FN++;;
	    }
        }
        $TN = $LENGTH - $FN;
    }
    else
    {
	# Scan both arrays to measure TP,FN,TN,FP
	for($i=0; $i<$LENGTH; $i++)
	{
	    if (${$ACOVERAGE{$key}}[$i] == $TRUE && ${$COVERAGE{$key}}[$i] == $TRUE)
            {
		$TP++;
	    }
            if (${$ACOVERAGE{$key}}[$i] == $FALSE && ${$COVERAGE{$key}}[$i] == $TRUE)
            {
		$FP++;
	    }
            if (${$ACOVERAGE{$key}}[$i] == $TRUE && ${$COVERAGE{$key}}[$i] == $FALSE)
            {
		$FN++;
	    }
            if (${$ACOVERAGE{$key}}[$i] == $FALSE && ${$COVERAGE{$key}}[$i] == $FALSE)
            {
		$TN++;
            }
        }
     }

     # Computing SN,SP and CC for current sequence
     if (($TP+$FP) == 0)
     {
	 $SN = $SP = $CC = 0;
     }
     else
     {
	$SN =  $TP/($TP+$FN); 
	$SP =  $TP/($TP+$FP); 
	$CC = (($TP*$TN) - ($FN*$FP)) / sqrt((($TP+$FN)*($TN+$FP)*($TP+$FP)*($TN+$FN)));
     }

     print_mess(join ("\t","TP","FP","FN","TN","|","SN","SP","CC"));     
     print_mess(join ("\t",$TP,$FP,$FN,$TN,"|",$SN,$SP,$CC,"\n"));    

     # Output values for current sequence
     printf("<TR><TD align=right>$key</TD><TD align=right>%d</TD><TD align=right>%d</TD><TD align=right>%d</TD><TD align=right>%d</TD><TD>&nbsp;</TD><TD align=right>%.2f</TD><TD align=right>%.2f</TD><TD align=right>%.2f</TD>" ,$TP,$FP,$FN,$TN,$SN,$SP,$CC);    

     $allSN = $allSN + $SN;
     $allSP = $allSP + $SP;
     $allCC = $allCC + $CC;
}

printf("<TR><TH>AVERAGE</TH><TD></TD><TD></TD><TD></TD><TD></TD><TD>&nbsp;</TD><TD align=right>%.2f</TD><TD align=right>%.2f</TD><TD align=right>%.2f</TD>" , $allSN/scalar(keys %ACOVERAGE), $allSP/scalar(keys %ACOVERAGE), $allCC/scalar(keys %ACOVERAGE));


print "</TABLE></TD></TABLE><BR><BR>";

### STEP 4. COMPUTING ACC at site level
print_mess("Step 4. Computing Accuracy at site level");

my $matches;
my $j;

# Check all of the sequences
print "<B><FONT COLOR=blue>Site Level Accuracy:</FONT></B><BR>";
print "<TABLE BORDER=1 BORDERCOLOR=red><TR><TD>";
print "<TABLE BORDER=0 CELLPADDING=10 CELLSPACING=10><TR><TH BGCOLOR=blue><FONT COLOR=white>LOCUS</TH><TH BGCOLOR=blue><FONT COLOR=white>sTP</TH><TH BGCOLOR=blue><FONT COLOR=white>sFP</TH><TH BGCOLOR=blue><FONT COLOR=white>sFN</TH><TH></TH><TH BGCOLOR=blue><FONT COLOR=white>sSN</TH><TH BGCOLOR=blue><FONT COLOR=white>sSP</TH><TH BGCOLOR=blue><FONT COLOR=white>sSN+sSP/2</TH>";

$allSN = $allSP = $allCC = 0;
foreach $key (sort keys %ANNOT)
{
    $TP=$FN=$FP=0;

    # A. Check all of the annotations of current sequence (TP,FN)
    for($i=0; $i<scalar(@{$ANNOT{$key}}); $i++)
    {
	$line = ${$ANNOT{$key}}[$i];
        
        @record = split(/\s+/,$line);

        $pos1 = $record[3];
        $pos2 = $record[4];

        $matches = 0;

        if (defined($COVERAGE{$key}))
        {
	    for($j=$pos1-1; $j<$pos2-1+1; $j++)
	    {
		if (${$COVERAGE{$key}}[$j] == $TRUE)
                {
		    $matches++;
		}
	    }
        }

        if ($matches >= (($pos2-$pos1+1)/$OVERLAP))
        {
	    $TP++;
	}
        else
        {
	    $FN++;
	}
    }


    # B. Check all of the predictions of current sequence, if any (FP)
    if (defined($PRED{$key}))
    {
        for($i=0; $i<scalar(@{$PRED{$key}}); $i++)
        {
	    $line = ${$PRED{$key}}[$i];
        
            @record = split(/\s+/,$line);

            $pos1 = $record[3];
            $pos2 = $record[4];

            $matches = 0;

	    for($j=$pos1-1; $j<$pos2-1+1; $j++)
	    {
		if (${$ACOVERAGE{$key}}[$j] == $TRUE)
                {
		    $matches++;
		}
	    }

            if ($matches < (($pos2-$pos1+1)/$OVERLAP))
            {
               $FP++;
            }
        }
    }

    # Computing SN,SP and CC for current sequence
    if (($TP+$FP) == 0)
    {
       $SN = $SP = $CC = 0;
    }
    else
    {
       $SN =  $TP/($TP+$FN); 
       $SP =  $TP/($TP+$FP); 
       $CC =  ($SN+$SP)/2;
    }

    print_mess(join ("\t","TP","FP","FN","|","SN","SP","SN+SP/2"));     
    print_mess(join ("\t",$TP,$FP,$FN,"|",$SN,$SP,$CC,"\n"));    

    # Output values for current sequence
    printf("<TR><TD>$key</TD><TD>%d</TD><TD>%d</TD><TD>%d</TD><TD>&nbsp;</TD><TD align=right>%.2f</TD><TD align=right>%.2f</TD><TD align=right>%.2f</TD>" ,$TP,$FP,$FN,$SN,$SP,$CC);    

    $allSN = $allSN + $SN;
    $allSP = $allSP + $SP;
    $allCC = $allCC + $CC;
}

printf("<TR><TH>AVERAGE</TH><TD></TD><TD></TD><TD></TD><TD>&nbsp;</TD><TD align=right>%.2f</TD><TD align=right>%.2f</TD><TD align=right>%.2f</TD>" , $allSN/scalar(keys %ACOVERAGE), $allSP/scalar(keys %ACOVERAGE), $allCC/scalar(keys %ACOVERAGE));


print "</TABLE></TD></TABLE><BR><BR>";
print "<B><FONT COLOR=green>Information about accuracy measures:</FONT></B><BR><BR>";

print "<B>Accuracy at Nucleotide Level</B><BR><br>";
print "<TABLE>";
print "<TR><TH ALIGN=LEFT>- TP (True positives)</TH><TD>Number of nucleotide positions in both known sites and predicted sites</TD>";
print "<TR><TH ALIGN=LEFT>- FN (False negatives)</TH><TD>Number of nucleotide positions in known sites but not in predicted sites</TD>";
print "<TR><TH ALIGN=LEFT>- FP (False positives)</TH><TD>Number of nucleotide positions not in known sites but in predicted sites</TD>";
print "<TR><TH ALIGN=LEFT>- TN (True negatives)</TH><TD>Number of nucleotide positions in neither known sites nor predicted sites</TD>";
print "<TR><TH ALIGN=LEFT></TH><TD></TD>";
print "<TR><TH ALIGN=LEFT>- SN (Sensitivity)</TH><TD>Amount of nucleotides from known sites that are supported by nucleotides from predicted sites</TD>";
print "<TR><TH ALIGN=LEFT></TH><TD><IMG SRC=/datasets/abs2005/data/sn.jpg></TD>";
print "<TR><TH ALIGN=LEFT>- SP (Specificity)</TH><TD>Amount of nucleotides from predicted sites that are supported by nucleotides from known sites</TD>";
print "<TR><TH ALIGN=LEFT></TH><TD><IMG SRC=/datasets/abs2005/data/sp.jpg></TD>";
print "<TR><TH ALIGN=LEFT>- CC (Correlation Coefficient)</TH><TD>Correlation value between SN and SP</TD>";
print "<TR><TH ALIGN=LEFT></TH><TD><IMG SRC=/datasets/abs2005/data/cc.jpg></TD>";
print "</TABLE><BR>";

print "<B>Accuracy at Site Level</B><BR><BR>";
print "<TABLE>";
print "<TR><TH ALIGN=LEFT>- sTP (True positives)</TH><TD>Number of known sites that are predicted</TD>";
print "<TR><TH ALIGN=LEFT>- sFN (False negatives)</TH><TD>Number of known sites that are not predicted</TD>";
print "<TR><TH ALIGN=LEFT>- sFP (False positives)</TH><TD>Number of predicted sites that are not known sites</TD>";
print "<TR><TH ALIGN=LEFT></TH><TD></TD>";
print "<TR><TH ALIGN=LEFT>- sSN (Sensitivity)</TH><TD>Amount of known sites that are supported by any predicted site</TD>";
print "<TR><TH ALIGN=LEFT></TH><TD><IMG SRC=/datasets/abs2005/data/ssn.jpg></TD>";
print "<TR><TH ALIGN=LEFT>- sSP (Specificity)</TH><TD>Amount of predicted sites that are supported by any known sites</TD>";
print "<TR><TH ALIGN=LEFT></TH><TD><IMG SRC=/datasets/abs2005/data/ssp.jpg></TD>";
print "<TR><TH ALIGN=LEFT>- (sSN+sSP)/2</TH><TD>Average value between sSN and sSP</TD>";
print "</TABLE><BR>";

print "-- To consider a match between a known site and a predicted site, an overlap of 0.25 or more is required -- <BR><BR>";

print "<B>References</B><BR><BR>";
print "- Burset, M. and Guigo, R. Evaluation of gene structure prediction programs. Genomics 34: 353-367 (1996).<BR><BR>";
print "- Tompa <I>et al.</I> Assessing computational tools for the discovery of transcription factor binding sites. Nature Biotechnology 23: 137-144 (2005).<BR><BR>";


### STEP 5. Prepare gff-records for posterior image in the OUTPUT_FILE
print_mess("Step 5. Printing into file the GFF records");

# Open OUTPUT_FILE to write the planted sites for a picture
open(OUT,">$OUTPUT_FILE") or die (print_error("Can't open $OUTPUT_FILE"));

foreach $key (sort keys %ANNOT)
{
    # A. Print all of the annotations from this locus name
    for($i=0; $i<scalar(@{$ANNOT{$key}}); $i++)
    {
	$line = ${$ANNOT{$key}}[$i];
        
        @record = split(/\s+/,$line);
        print OUT join("\t",$record[0]."_REAL",$record[1],$sANNOTATION,$record[3],$record[4],".",".",".","\n");
    }

    # B. Print all of the predictions from this locus name, if any
    if (defined($PRED{$key}))
    {
        for($i=0; $i<scalar(@{$PRED{$key}}); $i++)
        {
	    $line = ${$PRED{$key}}[$i];
        
            @record = split(/\s+/,$line);
            print OUT join("\t",$record[0]."_PRED",$record[1],$sPREDICTION,$record[3],$record[4],".",".",".","\n");
        }
    }
}



#######################################
### FUNCTIONS
#######################################

sub print_error
{
    my @mess = @_;

    print STDERR "!> CONSTRUCTOR\t",@mess,"\n";
    exit(0);
}

sub print_mess
{
    my @mess = @_;

    if (defined($opts{v}))
    {
	print STDERR ">",@mess,"\n";
    }
}

sub print_help
{
    if (defined($opts{h}))
    {
	print "#### Evaluation of pattern discovery 1.0\n";
	print "##################################\n";
	
	print "USAGE: evaluator <L> <file1> <file2>\n";
	print "<L>              Length of the sequences\n";
	print "<file1>          annotations (gff format)\n";
	print "<file2>          predictions (gff format)\n";

	exit(0);
    }
}
