#!/usr/local/bin/perl -w
# takes as first parameter path to nathalie real data file
# exported from Excel as TEXT standardized
# takes as second parameter the standardized output file of a prediction program
# WARNING!!! this file must NOT contain "concat" kind of sequences
# So clean it before with a "grep -v concat" before

# For each sequence in the real data file, compute gene begin and 
# gene end for the reality and the prediction

#============================================================

# check parameters
if ($#ARGV != 1)
  {
    die "need 2 parameters: real data file, and standardized output from gene prediction program";
  }

open (REALDATA, $ARGV[0]) || die "can't open $ARGV[0]";
open (PREDICTION, $ARGV[1]) || die "can't open $ARGV[1]";

# read real data and store genes
$prev_seq = 0;
$geneNb = 0;
print "Real Data\n\nSeq\tGene\tFrame\tbegin\tend\n";
while (<REALDATA>)
  {
    if (($seq_id, $type, $strand, $Lend, $Rend) = 
	($_ =~ /^seq(\d+)\s+(.{4})\s+([+-])\s+(\d+)\s+(\d+)/))
      {
	$seq_id = int($seq_id);
	$nbPred{$seq_id} = 0;
	$nbSplit{$seq_id} = 0;
	$nbConcat{$seq_id} = 0;
	$nbCorrect{$seq_id} = 0;
	$nbMissing{$seq_id} = 0;
	$nbWrong{$seq_id} = 0;
	
	# store exons
	if ($seq_id != $prev_seq)
	  {
	    $prev_seq = $seq_id;
	    $exonNb = 0;
	  }
	$real{$seq_id}{$strand}{$exonNb}{'type'} = $type;
	$real{$seq_id}{$strand}{$exonNb}{'Lend'} = $Lend;
	$real{$seq_id}{$strand}{$exonNb}{'Rend'} = $Rend;
	
	#compute gene
	push(@{$Rseq{$seq_id}{$geneNb}{'exons'}}, $exonNb);
	
	if ($type eq "Sngl")
	  {
	    printf "seq%03d\t$geneNb\t$strand\t$Lend\t$Rend\n", $seq_id;
	    $Rseq{$seq_id}{$geneNb}{'strand'} = $strand;
	    $Rseq{$seq_id}{$geneNb}{'begin'} = $Lend;
	    $Rseq{$seq_id}{$geneNb}{'end'} = $Rend;
	    $geneNb++;
	  }
	else
	  {
	    $beginLabel = $strand eq '+'?'Init':'Term';
	    $endLabel = $strand eq '+'?'Term':'Init';
	    
	    if ($type eq $beginLabel)
	      {
		printf "seq%03d\t$geneNb\t$strand\t$Lend\t", $seq_id;
		$Rseq{$seq_id}{$geneNb}{'strand'} = $strand;
		$Rseq{$seq_id}{$geneNb}{'begin'} = $Lend;
	      }
	    elsif ($type eq $endLabel)
	      {
		print "$Rend\n";
		$Rseq{$seq_id}{$geneNb}{'end'} = $Rend;
		$geneNb++;
	      }
	  }
      }
    $exonNb++;
  }

close(REALDATA);

# read predictions

$prev_seq = 0;
while (<PREDICTION>)
  {
    if (($seq_id, $type, $strand, $Lend, $Rend) = 
	($_ =~ /^seq(\d+)\s+(.{4})\s+([+-])\s+(\d+)\s+(\d+)/))
      {
	$seq_id = int($seq_id);
	if ($prev_seq != $seq_id)
	  {
	    $prev_seq = $seq_id;
	    $exonNb = 0;
	  }
	$pred{$seq_id}{$strand}{$exonNb}{'type'} = $type;
	$pred{$seq_id}{$strand}{$exonNb}{'Lend'} = $Lend;
	$pred{$seq_id}{$strand}{$exonNb}{'Rend'} = $Rend;
	$exonNb++;
      }
  }
close(PREDICTION);

#  compute predicted genes
print "\nPrediction\n\nSeq\tGene\tFrame\tbegin\tend\n";
@sequences = sort {$a <=> $b} keys %pred;
$geneNb = 0;
foreach $s (@sequences)
  {
    foreach $strand ('+','-')
      {
	if (!exists($pred{$s}{$strand})){next;};
	
	@exons = sort {$pred{$s}{$strand}{$a}{'Lend'} <=> $pred{$s}{$strand}{$b}{'Lend'}} keys %{$pred{$s}{$strand}};
	
	$inGene = 0;
	$beginLabel = $strand eq '+' ? 'Init' : 'Term';
	$endLabel = $strand eq '+' ? 'Term' : 'Init';
	
	foreach $e (@exons)
	  {
	    $Lend = $pred{$s}{$strand}{$e}{'Lend'};
	    $Rend = $pred{$s}{$strand}{$e}{'Rend'};
	    $type = $pred{$s}{$strand}{$e}{'type'};
	    
	    if ($type eq 'Sngl')
	      { # new single gene
		if ($inGene == 1)
		  { # ends previous gene
		    $Pseq{$s}{$geneNb}{'end'} = $pred{$s}{$strand}{$prev_ex}{'Rend'};
		    $geneNb++;
		  }
		$Pseq{$s}{$geneNb}{'strand'} = $strand;
		$Pseq{$s}{$geneNb}{'begin'} = $Lend;
		$Pseq{$s}{$geneNb}{'end'} = $Rend;
		push(@{$Pseq{$s}{$geneNb}{'exons'}}, $e);
		$geneNb++;
		$inGene = 0;
	      }
	    elsif ($type eq $beginLabel)
	      { # new gene
		if ($inGene == 1)
		  { # ends previous gene
		    $Pseq{$s}{$geneNb}{'end'} = $pred{$s}{$strand}{$prev_ex}{'Rend'};
		    $geneNb++;
		  }
		$Pseq{$s}{$geneNb}{'strand'} = $strand;
		$Pseq{$s}{$geneNb}{'begin'} = $Lend;
		push(@{$Pseq{$s}{$geneNb}{'exons'}}, $e);
		$inGene = 1;
	      }
	    elsif ($type eq $endLabel)
	      { #end of a gene
		if ($inGene == 0)
		  { # current gene not initialized
		    $Pseq{$s}{$geneNb}{'strand'} = $strand;
		    $Pseq{$s}{$geneNb}{'begin'} = $Lend;
		  }
		$Pseq{$s}{$geneNb}{'end'} = $Rend;
		push(@{$Pseq{$s}{$geneNb}{'exons'}}, $e);
		$geneNb++;
		$inGene = 0;
	      }
	    else # Intr exon
	      {
		if ($inGene == 0)
		  { # current gene not initialized
		    $Pseq{$s}{$geneNb}{'strand'} = $strand;
		    $Pseq{$s}{$geneNb}{'begin'} = $Lend;
		  }
		push(@{$Pseq{$s}{$geneNb}{'exons'}}, $e);
		$inGene = 1;
	      }
	    $prev_ex = $e;
	  }
	if ($inGene == 1 && exists($Pseq{$s}{$geneNb}{'begin'}) && !exists($Pseq{$s}{$geneNb}{'end'}))
	  {
	    # current gene not ended
	    $Pseq{$s}{$geneNb}{'end'} = $Rend;
	    $geneNb++;
	  }
      }
  }

@sequences = sort {$a <=> $b} keys %Pseq;
foreach $s (@sequences)
  {
    @gene = sort {$Pseq{$s}{$a}{'begin'} <=> $Pseq{$s}{$b}{'begin'}} keys %{$Pseq{$s}};
    foreach $g (@gene)
      {
	$strand = $Pseq{$s}{$g}{'strand'};
	@pexons = sort {$pred{$s}{$strand}{$a}{'Lend'} <=> $pred{$s}{$strand}{$b}{'Lend'}} @{$Pseq{$s}{$g}{'exons'}};
	printf "seq%03d\t$g\t%s\t%d\t%d\t%s\t%s\n", 
	$s, $Pseq{$s}{$g}{'strand'}, $Pseq{$s}{$g}{'begin'}, $Pseq{$s}{$g}{'end'},
	$pred{$s}{$strand}{$pexons[0]}{'type'}, $pred{$s}{$strand}{$pexons[$#pexons]}{'type'};
      }
  }


# compute ovelapping
@sequences = sort {$a <=> $b} keys %Rseq;
foreach $s (@sequences)
  {
    @Rgene = sort {$Rseq{$s}{$a}{'begin'} <=> $Rseq{$s}{$b}{'begin'}} keys %{$Rseq{$s}};
    
    if (!exists($Pseq{$s}))
      {
	# sequence not predicted at all, so no overlapping
	next;
      }
    @Pgene = sort {$Pseq{$s}{$a}{'begin'} <=> $Pseq{$s}{$b}{'begin'}} keys %{$Pseq{$s}};
    
    foreach $rg (@Rgene)
      {
	foreach $pg (@Pgene)
	  {
	    if ($Rseq{$s}{$rg}{'strand'} ne $Pseq{$s}{$pg}{'strand'}) {next;} #skip gene on different strands
	    
	    # check if prediction overlap reality
	    if (($Rseq{$s}{$rg}{'begin'} <= $Pseq{$s}{$pg}{'begin'} &&
		 $Pseq{$s}{$pg}{'begin'} <= $Rseq{$s}{$rg}{'end'}) ||
		($Pseq{$s}{$pg}{'begin'} <= $Rseq{$s}{$rg}{'begin'} &&
		 $Rseq{$s}{$rg}{'begin'} <= $Pseq{$s}{$pg}{'end'}))
	      {
		push(@{$spt{$s}{$rg}}, $pg); # real gene overlap a prediction (split)
		push(@{$cnt{$s}{$pg}}, $rg); # prediction overlap a real gene (concat)
	      }
	  }
      }
  }

print "\nSplit and Concat\n\nContig\tReal gene(s)\tSplit in/concatenated with\n";

# compute split and concat
@sequences = sort {$a <=> $b} keys %spt;
foreach $s (@sequences) 
  {
    #split
    @gene = sort {$a <=> $b} keys %{$spt{$s}};
    foreach $g (@gene)
      {
	@l = @{$spt{$s}{$g}};
	if ($#l != 0)
	  {
	    printf "seq%03d\t%d\t(%s)\n", $s, $g, join(',',@l);
	    $nbSplit{$s}++;
	  }
      }
    #concat
    @gene = sort {$a <=> $b} keys %{$cnt{$s}};
    foreach $g (@gene)
      {
	@l = @{$cnt{$s}{$g}};
	if ($#l != 0)
	  {
	    printf "seq%03d\t(%s)\t%d\n", $s, join(',',@l), $g;
	    $nbConcat{$s} += 1 + $#l;
	  }
      }
  }

# compute Nb of wrong predicted genes (not overlapping a real gene)
foreach $s (keys %Pseq)
  {
    @gene = keys %{$Pseq{$s}};
    $nbPred{$s} = $#gene + 1;
    foreach $g (@gene)
      {
	if (!exists($cnt{$s}{$g}))
	  {
	    $nbWrong{$s}++;
	  }
      }
  }

# compute Nb of genes per sequence, Nb correctly predicted, Nb not prediected at all
print "\nCompute correctness of prediction\n\nContig\tGene\tComment\n";
foreach $s (sort {$a <=> $b} keys %Rseq)
  {
    @rgene = sort {$a <=> $b} keys %{$Rseq{$s}};
    $nbGene{$s} = $#rgene + 1;
    foreach $rg (@rgene)
      {
	if (!exists($spt{$s}{$rg}))
	  {
	    printf "$s\tGene[%+6s,%+6s]\tnot predicted\n", $Rseq{$s}{$rg}{'begin'}, $Rseq{$s}{$rg}{'end'};
	    $nbMissing{$s}++;
	  }
	else
	  {
	    foreach $pg (@{$spt{$s}{$rg}})
	      {
		$strand = $Rseq{$s}{$rg}{'strand'};
		@rexons = sort {$real{$s}{$strand}{$a}{'Lend'} <=> $real{$s}{$strand}{$b}{'Lend'}} @{$Rseq{$s}{$rg}{'exons'}};
		@pexons = sort {$pred{$s}{$strand}{$a}{'Lend'} <=> $pred{$s}{$strand}{$b}{'Lend'}} @{$Pseq{$s}{$pg}{'exons'}};
		if ($#rexons != $#pexons)
		  {
		    # not the same number of exons
		    next;
		  }
		for ($i=0; $i<=$#rexons; $i++)
		  {
		    if ($real{$s}{$strand}{$rexons[$i]}{'Lend'} != $pred{$s}{$strand}{$pexons[$i]}{'Lend'} ||
			$real{$s}{$strand}{$rexons[$i]}{'Rend'} != $pred{$s}{$strand}{$pexons[$i]}{'Rend'})
		      {
			# exon borders differ
			last;
		      }
		  }
		
		if ($i > $#rexons)
		  {
		    # all exons are OK, the gene has been found
		    printf "$s\tGene[%+6s,%+6s]\tperfect match\t", 
		    $Rseq{$s}{$rg}{'begin'}, $Rseq{$s}{$rg}{'end'};
		    
		    # check type of exons border
		    if ($#rexons == 0)
		      {
			if ($pred{$s}{$strand}{$pexons[0]}{'type'} cmp 'Sngl')
			  {
			    print "NOT single predicted\n";
			  }
			else
			  {
			    print "Good single predicted\n";
			  }
		      }
		    else # mulptiple exons per gene
		      {
			if ($real{$s}{$strand}{$rexons[0]}{'type'} cmp $pred{$s}{$strand}{$pexons[0]}{'type'} ||
			    $real{$s}{$strand}{$rexons[$#rexons]}{'type'} cmp $pred{$s}{$strand}{$pexons[$#rexons]}{'type'})
			  {
			    printf "at least one border exon has incorrect type: %s\t%s\n",
			    $pred{$s}{$strand}{$pexons[0]}{'type'}, $pred{$s}{$strand}{$pexons[$#rexons]}{'type'};
			  }
			else
			  {
			    print "Perfect model\n";
			  }
		      }			
		    
		    $nbCorrect{$s}++;
		    last;
		  }
	      }
	  }
      }
  }

# summary
print "\nSummary\n\nContig\t#genes\t#pred\t#Correct\t#Missing\t#partial\t#Wrong\t#split\t#concat\n";
@sequences = sort {$a <=> $b} keys %Rseq;
foreach $s (@sequences)
  {
    printf "seq%03d\t%d\t%d\t%d\t%d\t%d\t%d\t%d\t%d\n", $s, $nbGene{$s}, $nbPred{$s}, $nbCorrect{$s}, $nbMissing{$s}, 
    $nbGene{$s} - $nbCorrect{$s} - $nbMissing{$s}, $nbWrong{$s}, $nbSplit{$s}, $nbConcat{$s};
  }
