#!/bin/perl
use strict;
use warnings;

my $dat_ct="/u/home/praktikum/rna_benni/RNA_STRAND_data/datenbank_ct/"; #Ordner der Rohdaten aus Datenbank
my $cofold_ct="/u/home/praktikum/rna_benni/RNA_STRAND_data/cofold_ct/"; #Ordner der cofold-Daten


opendir(my $datliste, $dat_ct)||die "can not opendir,schrei um Hilfe";
my @datenliste = grep {/\.ct$/} readdir($datliste);
closedir $datliste;

opendir(my $coliste, $cofold_ct)||die "can not opendir, schrei um Hilfe"; #Erstelle Liste mit Seqenzdaten
my @cofoldliste = grep {/\.ct$/} readdir($coliste);
closedir $coliste;


#Lese Cofoldliste in Tabelle/array ein

my %FullNameDatenliste;
my %FullNameCofoldliste;

for (my $i=0; $i<=$#datenliste; $i++) {
    my $input=$datenliste[$i];
    my $name=$input;
    $name =~ s/\..*$//;
    $FullNameDatenliste{$name}=$input;
}

my @jobliste=();

for (my $i=0; $i<=$#cofoldliste; $i++) {
    my $input=$cofoldliste[$i];
    my $name=$input;
    $name =~ s/\..*$//; # input ist der Basisname einer datei in Cofoldliste
    if (exists $FullNameDatenliste{$name}) {
	#print "$input komt in beiden Listen vor\n";
	$FullNameCofoldliste{$name}=$input;
	push @jobliste,$name;
    }
}

#print "------------------------------------------------------------\n";
#print "Eintraege in jobliste: ".(@jobliste)."\n";

#Schreibe Datein mit gleichen Namen, gepaart in eine Liste

foreach my $name (@jobliste) {
    my $inputDaten=$FullNameDatenliste{$name};
    my $inputCofold=$FullNameCofoldliste{$name};
    
    my ($TOTAL,$lenge) = read_ct_seq("$cofold_ct$inputCofold");
    my $TP=  TP_ct("$dat_ct$inputDaten","$cofold_ct$inputCofold", $TOTAL);
    my $FP= FP_ct("$dat_ct$inputDaten","$cofold_ct$inputCofold", $TOTAL);
    my $FN=  FN_ct("$dat_ct$inputDaten","$cofold_ct$inputCofold", $TOTAL);
    
    #   print "vergleiche(\"$dat_ct$inputDaten\",\"$cofold_ct$inputCofold\")\n";
       print "  TP = $TP\n";
       print "  FP = $FP\n";
       print "  FN = $FN\n";
    
    my $TN=($TOTAL-$TP-$FP-$FN);
    print "  TN = $TN\n";
    
#Berechnug von Vergleichsfaktoren!
    
 #   print "lenge=$lenge\n";
    
    
    my $Sensitivity=sprintf("%.5f", ($TP/($TP+$FN)));
    print "Sensitivity = $Sensitivity\n";
    
    my $FPR=($FP/($FP+$TN));
    print "false positive rate = $FPR\n";
    
    my $Specificity=sprintf("%.5f",($TN/($FP+$TN)));
    print "Specificity = $Specificity\n";
    
    my $PPV=($TP/($TP+$FP));
   print "positive predictive value = $PPV\n";
    
    my $NPV=($TN/($TP+$FN));
    print "negative predictive value = $NPV\n";
    
    my $FDR=($FP/($FP+$TP));
    print "false discovery rate = $FDR\n";
    
    my $MCC=sprintf("%.5f",(($TP*$TN)-($FP*$FN))/((($TP+$FP)*($TP+$FN)*($TN+$FP)*($TN+$FN))**0.5));
    print "Matthews correlation coefficient = $MCC\n";
    
    my $N=($TP+$TN+$FN+$FP);
  
    my $P=(($TP+$FP)/$N);
    my $ACC=(($TP+$TN)/($P+$N));
 #   print  "accuracy = $ACC\n";
    
    my $F1=(($TP*$TP)/(($TP*$TP)+$FP+$FN));
    print "F1 score = $F1\n";

 #   print "$lenge\t$NPV\n";
    

}


sub TP_ct {
    my ($refname,$predname, $total)=@_;
   # print "vergleiche(\"$refname\",\"$predname\")\n";
    
    my @reflist = read_ct_bps($refname);
    my @predlist = read_ct_bps($predname);

    #my $total = read_ct_seq($predname);

    # bestimme true positives
    my $tp=0;
    for (my $i=0; $i<@reflist; $i++) {
	if ($reflist[$i] > $i) {
	    if ($predlist[$i]==$reflist[$i]) {
		$tp++;
	    }
	}
    }
  # print "  TP = $tp\n";
    return $tp;}
    
sub FP_ct {
    # bestimme false positives
    my ($refname,$predname, $total)=@_;
   # print "vergleiche(\"$refname\",\"$predname\")\n";
    
    my @reflist = read_ct_bps($refname);
    my @predlist = read_ct_bps($predname);

    # my $total = read_ct_seq($predname);
    my $fp=0;
    for (my $i=0; $i<@predlist; $i++) {
	if ($predlist[$i] > $i) {
	    if ($predlist[$i]!=$reflist[$i]) {
		$fp++;
	    }
	}
    }#print "  FP = $fp\n";
    return $fp;
     }

sub FN_ct {
    # bestimme false positives
    my ($refname,$predname, $total)=@_;
   # print "vergleiche(\"$refname\",\"$predname\")\n";
    
    my @reflist = read_ct_bps($refname);
    my @predlist = read_ct_bps($predname);

    #my $total = read_ct_seq($predname);
    my $fn=0;
    for (my $i=0; $i<@predlist; $i++) {
	if ($predlist[$i] > $i) {
	    if ($predlist[$i]!=$reflist[$i]) {
		$fn++;
	    }
	}
    }#print "  FN = $fn\n";
    return $fn;
     }
sub read_ct_bps {
    my ($filename) = @_;
    
    my @bps=();
    
    open(my $in, $filename);
    my @lines = (<$in>);
    close $in;

    foreach my $line (@lines) {
	next if($line=~/^$/);	
if ($line!~/^#/ && $line!~/ENERGY/i ) {
	    my @entries = split /\s+/,$line;
	    push @bps, $entries[5];
	}
    }
        
    # print "read_ct_bps($filename): @bps\n";
    
    return @bps;
}
sub read_ct_seq {
    my ($filename) = @_;
    my %count=("A"=>0,"UT"=>0,"G"=>0,"C"=>0)
	       ;

    my @seq=();
    
    open(my $in, $filename);
    my @lines = (<$in>);
    close $in;
    
    foreach my $line (@lines) {
	next if($line=~/^$/);
	if ($line!~/^#/ && $line!~/ENERGY/i ) {
	    my @entries = split /\s+/,$line;
	    push @seq, $entries[2];
	}
    }
#count A
    for my $seq (@seq){
	if ($seq eq "A"|| $seq eq "a"){
	    $count{"A"}++;
	}
    }
     print "A_base = $count{'A'}\n";
    
#count T
    
    for my $seq (@seq){
	if ($seq eq "T" || $seq eq "U"|| $seq eq "u"||$seq eq "t"){
	    $count{"UT"}++;
	}
    }
print "U_base = $count{'UT'}\n";
    
#count G
    
    for my $seq (@seq){
	if ($seq eq "G"|| $seq eq "g"){
	    $count{"G"}++;
	}
    }
    
print "G_base = $count{'G'}\n";
    
#count C
    
    for my $seq (@seq){
	if ($seq eq "C"|| $seq eq "c"){
	    $count{"C"}++;
	}
    }
     print "C_base = $count{'C'}\n";
    
#Berechne TN,als TN= AxT=GxC-TP-FP-FN ,TO DO
    
    my $TN=($count{"A"}*$count{"UT"})+($count{"G"}*$count{"C"});
    
    #print " TN = $TN\n";
    my $lenge=($count{"A"}+$count{"UT"}+$count{"G"}+$count{"C"});

  #  print "$lenge\t$count{'A'}\t$count{'UT'}\t$count{'G'}\t$count{'C'}\n";

    # print "read_ct_bps($filename): @bps\n";
    return ($TN,$lenge);
}