#!/usr/local/bin/perl -w

# takes 2 standard files as input: one real, one prediction 
# and a data list file
# compute the TN value for each non coding region of a sequence:
# number of nucleic acid in a non coding region that are predicted 
# as non coding

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

sub max
  {
    return $_[0] > $_[1] ? $_[0] : $_[1];
  }

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

sub min
  {
    return $_[0] < $_[1] ? $_[0] : $_[1];
  }

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

if ($#ARGV != 2)
  {
    die "usage: param1: standard real data param2: standard prediction, param3: data list file";
  }

open(RDATA, $ARGV[0]) || die "can't open $ARGV[0]";
open(PDATA, $ARGV[1]) || die "can't open $ARGV[1]";
open(LDATA, $ARGV[2]) || die "can't open $ARGV[2]";

# read sequence length
while (<LDATA>)
  {
    if (($seq_id, $Length) = (/^seq(\d+)\s+(\d+)/))
      {
	$seq_length{int($seq_id)} = $Length;
      }
  }
close(LDATA);

# read real data and build non coding regions
$prevSeq = 0;
while (<RDATA>)
  {
    if (($seq_id, $Lend, $Rend) =
	(/^seq(\d+)\s+.{0,4}\s+[+-]\s+(\d+)\s+(\d+)\s+\d+\s+[+-]?\d?\s+[+-]\d/))
      {
	$seq_id = int($seq_id);
	
	if ($prevSeq != $seq_id)
	  {
	    if ($prevSeq != 0 && $prevLend != $seq_length{$prevSeq})
	      {
		# handle end of previous sequence
		$rintron{$prevSeq}{$exonNb}{'Lend'} = $prevLend;
		$rintron{$prevSeq}{$exonNb}{'Rend'} = $seq_length{$prevSeq};
	      }
	    $exonNb = 0;
	    $prevLend = 1;
	    $prevSeq = $seq_id;
	  }
	
	# compute borders of non coding region
	if ($Lend != 1)
	  {
	    # exon not starting at the begining of the sequence
	    $rintron{$seq_id}{$exonNb}{'Lend'} = $prevLend;
	    $rintron{$seq_id}{$exonNb}{'Rend'} = $Lend - 1;
	  }
	$prevLend = $Rend + 1;
	
	$exonNb++;
      }
  }
close(RDATA);

# handle end of last sequence
$rintron{$prevSeq}{$exonNb}{'Lend'} = $prevLend;
$rintron{$prevSeq}{$exonNb}{'Rend'} = $seq_length{$prevSeq};

# read prediction
$prevSeq = 0;
while (<PDATA>)
  {
    if (($seq_id, $Lend, $Rend) =
	(/^seq(\d+)\s+.{0,4}\s+[+-]\s+(\d+)\s+(\d+)\s+\d+\s+[+-]?\d?\s+[+-]\d/))
      {
	$seq_id = int($seq_id);
	if ($prevSeq != $seq_id)
	  {
	    $exonNb = 0;
	    $prevSeq = $seq_id;
	  }
	$pdata{$seq_id}{$exonNb}{'Lend'} = $Lend;
	$pdata{$seq_id}{$exonNb}{'Rend'} = $Rend;
	$exonNb++;
      }
  }
close(PDATA);

# concat overlapping prediction
@sequences = sort {$a <=> $b} keys %pdata;
$prevSeq = 0;
foreach $s (@sequences)
  {
    @Pexons = sort {$pdata{$s}{$a}{'Lend'} <=> $pdata{$s}{$b}{'Lend'}} keys %{$pdata{$s}};
    $e1 = shift(@Pexons);
    @Lexons = ($e1);
    foreach $e2 (@Pexons)
      {
	if ($pdata{$s}{$e1}{'Lend'} <= $pdata{$s}{$e2}{'Lend'} &&
	    $pdata{$s}{$e2}{'Lend'} <= $pdata{$s}{$e1}{'Rend'})
	  {
	    # overlapping regions
	    $pdata{$s}{$e1}{'Rend'} = max($pdata{$s}{$e1}{'Rend'},$pdata{$s}{$e2}{'Rend'});
	  }
	else
	  {
	    # really new predicted coding region
	    push(@Lexons, $e2); 
	    $e1 = $e2;
	  }
      }
    
    # build non coding regions
    @Lexons = sort {$pdata{$s}{$a}{'Lend'} <=> $pdata{$s}{$b}{'Lend'}} @Lexons;
    foreach $e (@Lexons)
      {
	if ($prevSeq != $s)
	  {
	    if ($prevSeq != 0 && $prevLend != $seq_length{$prevSeq})
	      {
		# handle end of previous sequence
		$pintron{$prevSeq}{$exonNb}{'Lend'} = $prevLend;
		$pintron{$prevSeq}{$exonNb}{'Rend'} = $seq_length{$prevSeq};
	      }
	    $exonNb = 0;
	    $prevLend = 1;
	    $prevSeq = $s;
	  }
	
	# compute borders of non coding region
	if ($pdata{$s}{$e}{'Lend'} != 1)
	  {
	    # exon not starting at the begining of the sequence
	    $pintron{$s}{$exonNb}{'Lend'} = $prevLend;
	    $pintron{$s}{$exonNb}{'Rend'} = $pdata{$s}{$e}{'Lend'} - 1;
	  }
	$prevLend = $pdata{$s}{$e}{'Rend'} + 1;
	
	$exonNb++;
      }

    # handle end of last sequence
    $pintron{$prevSeq}{$exonNb}{'Lend'} = $prevLend;
    $pintron{$prevSeq}{$exonNb}{'Rend'} = $seq_length{$prevSeq};
  }

# compute ovelapping introns
@sequences = sort {$a <=> $b} keys %rintron;
foreach $s (@sequences)
  {
    @IR = sort {$rintron{$s}{$a}{'Lend'} <=> $rintron{$s}{$b}{'Lend'}} keys %{$rintron{$s}};
    foreach $ir (@IR)
      {
	if (!exists($pintron{$s}))
	  {
	    next;
	  }
	@IP = sort {$pintron{$s}{$a}{'Lend'} <=> $pintron{$s}{$b}{'Lend'}} keys %{$pintron{$s}};
	foreach $ip (@IP)
	  {
	    if (($rintron{$s}{$ir}{'Lend'} <= $pintron{$s}{$ip}{'Lend'} &&
		 $pintron{$s}{$ip}{'Lend'} <= $rintron{$s}{$ir}{'Rend'}) ||
		($pintron{$s}{$ip}{'Lend'} <= $rintron{$s}{$ir}{'Lend'} &&
		 $rintron{$s}{$ir}{'Lend'} <= $pintron{$s}{$ip}{'Rend'}))
	      {
		push (@{$overlap{$s}{$ir}}, $ip);
	      }
	  }
      }
  }

# generates output
print "Contig\t";
print "Non coding (real)\t";
print "Length\t";
print "Non coding (pred)\t";
print "TN\t";
print "%TN\t";
print "FN\t";
print "%FN\n";
@sequences = sort {$a <=> $b} keys %rintron;
foreach $s (@sequences)
  {
    @introns = sort {$rintron{$s}{$a}{'Lend'} <=> $rintron{$s}{$b}{'Lend'}} keys %{$rintron{$s}};
    foreach $ir (@introns)
      {
	$L = $rintron{$s}{$ir}{'Lend'};
	$R = $rintron{$s}{$ir}{'Rend'};
	
	if (!exists($overlap{$s}))
	  {
	    # sequenced with no prediction
	    printf "seq%03d\t[%5d,%5d]\t%d\t0\t0.00\n", $s, $L, $R, $R - $L + 1;
	  }
	else 
	  {
	    foreach $pi (@{$overlap{$s}{$ir}})
	      {
		# Contig
		printf "seq%03d\t", $s;
		# Non coding (real)
		printf "[%5d,%5d]\t", $L, $R;
		# Length
		$Length =  $R - $L + 1;
		printf "%d\t", $Length;
		# Non coding (pred)
		$Lp = $pintron{$s}{$pi}{'Lend'};
		$Rp = $pintron{$s}{$pi}{'Rend'};
		printf  "[%5d,%5d]\t", $Lp, $Rp;
		# TN
		$TN = &min($R,  $Rp) - &max($L,  $Lp) + 1;
		printf "%d\t", $TN;
		# %TN
		printf "%.2f\t", $TN / $Length * 100;
		$FN = $Rp - $Lp + 1 - $TN;
		# FN
		printf "%d\t", $FN;
		# %FN
		printf "%.2f\n", $FN / $Length * 100;
	      }
	  }
      }
  }
