#!/usr/bin/perl -w
use strict;
my $chain='';
my $opm_coord='';
my $membrane=0;
my $length=0;
my $protein='';
print "Enter file name with OPM annotation\n";
my $fname=<STDIN>;
chomp $fname;
if (-e $fname)
{open FILE,$fname;}
else {print "There is no such file!\n";
		die;}
my @opm=<FILE>;
close FILE;
foreach my$i(@opm)
	{
	if ($i=~m/Topology.+subunit\s(\w)/)
		{$chain=$1;
		last;}
	}
foreach my$i(@opm)
	{
	if ($i=~m/$chain.+Segments: (.+)/)
		{
		$opm_coord=$1;
		last;
		}
	}
$opm_coord=~s/\d+\(//g;
$opm_coord=~s/\D/ /g;
my @opm_coord=split ' ',$opm_coord;
print "Enter file name with TMHMM prediction\n";
$fname=<STDIN>;
chomp $fname;
if (-e $fname)
{open FILE,$fname;}
else {print "There is no such file!\n";
		die;}
my @tmhmm=<FILE>;
close FILE;
my @tmhmm_coord;
foreach my$i(@tmhmm)
	{
	if ($i=~m/#\s(.+)\sLength:\s(\d+)/)
		{
		$protein=$1;
		$length=$2;
		}
	if ($i=~m/TMhelix\s+(\d+)\s+(\d+)/)
		{
		push @tmhmm_coord,$1,$2;
		$membrane+=$2-$1+1;
		}
	}
print "Enter file name with needle alignment\n";
$fname=<STDIN>;
chomp $fname;
if (-e $fname)
{open FILE,$fname;}
else {print "There is no such file!\n";
		die;}
my @lines=<FILE>;
close FILE;
my $seq_tmhmm='';
my $seq_opm='';
foreach my$i(@lines)
	{
	if ($i=~m/$protein\s+\d+\s([\w-]+)\s+\d+/)
		{$seq_tmhmm.=$1;}
	elsif ($i=~m/\w+_\w+\s+\d+\s([\w-]+)\s+\d+/)
		{$seq_opm.=$1;}
	}
my $gap=0;
my $i=0;
while ()
	{
	$gap=index($seq_tmhmm,'-');
	if ($gap != -1)
		{
		for (my $i=0;$i<@tmhmm_coord;$i++)
			{
			if ($tmhmm_coord[$i]>$gap+1)
				{$tmhmm_coord[$i]++;}
			}
		substr($seq_tmhmm,$gap,1)='';
		}
	else {last;}
	}
while ()
	{
	$gap=index($seq_opm,'-');
	if ($gap != -1)
		{
		for (my $i=0;$i<@opm_coord;$i++)
			{
			if ($opm_coord[$i]>$gap+1)
				{$opm_coord[$i]++;}
			}
		substr($seq_opm,$gap,1)='';
		}
	else {last;}
	}
my $tn=0;
my $tp=0;
my $fn=0;
my $fp=0;
my $left=0;
my $right=0;
my $j=0;
for ($i=0;$i+1<@tmhmm_coord;$i+=2)
	{
	for ($j=0;$j+1<@opm_coord;$j+=2)
		{
		if ($tmhmm_coord[$i+1]<$opm_coord[$j] || $tmhmm_coord[$i]>$opm_coord[$j+1])
			{}
		else {
			$tp+=$tmhmm_coord[$i+1]-$tmhmm_coord[$i]+1;
			$left=$tmhmm_coord[$i]-$opm_coord[$j];
			$right=$tmhmm_coord[$i+1]-$opm_coord[$j+1];
			if ($left>0) {$fn+=$left;}
			else {$fp-=$left;
				$tp+=$left;}
			if ($right<0) {$fn-=$right;}
			else {$fp+=$right;
				$tp-=$right;}
			}
		}
	}
$tn=$length-$tp-$fp-$fn;
open FILE,">results.txt";
print FILE "Length: $length\nMembrane aa: $membrane\nTP: $tp\nFP: $fp\nTN: $tn\nFN: $fn\nSensivity: ".$tp/($tp+$fn)."\nSpecificity: ".$tn/($tn+$fp)."\nPrecision: ".$tp/($tp+$fp)."\nOverprediction: ".$fp/($fp+$tp)."\nUnderprediction: ".$fn/($tn+$fn)."\n";
close FILE;
