#!/usr/bin/perl -w

use strict;
use File::Temp qw(tempfile);
use Getopt::Long;
use warnings;
use feature qw/say switch/;

sub usage{
    print "USAGE:\n./ICI_sno.pl  -p|--pred <target-prediction-file> -f|--fasta <path to mfastas> -o|--out <output-file> [-a|-aln <file with start&end-pos of seq in alignment>]\n";
    exit;
}

    
my ($pred,$out,$path,$aln_start) = ('','','','');
my $help;

GetOptions('h|help|?' => \$help,
	   'p|pred=s' => \$pred,
	   'f|fasta=s'  => \$path,
	   'a|aln=s' => \$aln_start,
	   'o|out=s'  => \$out);

if(defined $help || $pred eq ''){&usage;}

if($path eq ''){$path = "/scr/vodka/sebastian/snoRNA/mfasta/"}
# if($aln_start eq ''){$aln_start = "/scr/rioja/steffi/PROJECTS/TARGETS/RNAsalsa_new/aln_seq_starts";}
#alle spezies mit homologen der Familie
my @tmp = split(/\//,$pred);
@tmp = split(/_/,$tmp[-1]);
my $homs = &findAllHoms($tmp[0]."_".$tmp[1]);
my $nrHoms = scalar(split/,/,$homs);

open(PRED , "<".$pred) or die $!; 

my $hash;

while(<PRED>){
    my @line = split(/\s+/,$_);
    my $spez = $line[0];
    my $pos = $line[3];
    my $mfe = $line[4];
    if($mfe eq 'NA'){next;}
    #neuer Organismus
    if(!defined $hash->{$spez}){
	#individuelle Interaction speichern
	$hash->{$spez}->{$pos}->{mfe} = $mfe;

	#mfes pro viech aufsummieren und zaehlen
	$hash->{$spez}->{sumMFE} = $mfe;
	$hash->{$spez}->{cnt} = 1;

	#um average energy pro viech zu kennen
	#$hash->{$spez}->{avgMFE} = $mfe;
    }
    else{
	#individuelle Interaction speichern
	$hash->{$spez}->{$pos}->{mfe} = $mfe;
	
	#mfes summieren und anzahl der interactions zaehlen
	$hash->{$spez}->{sumMFE} += $mfe;
	$hash->{$spez}->{cnt}++;

	#$hash->{$spez}->{avgMFE} =  $hash->{$spez}->{sumMFE} / $hash->{$spez}->{cnt};
    }
    #print "$spez \n";	
}

close(PRED);

#calculate the ICI values for each target that can bind to the snoRNA in at least one species
my $ICI; 

#my $anz = 0;

foreach my $org (keys %{$hash}){
    #in wievielen Spezies hab ich predictions
     #$anz++;
    #if( $hash->{$org}->{cnt} > 1){
	#calculate the average interaction energy in one species
	my $avgMFE = $hash->{$org}->{sumMFE} / $hash->{$org}->{cnt};
	#print "$org \t $avgMFE\n";

	#for all interactions in the given target-file
	foreach my $mod (keys %{$hash->{$org}}){
	    #mod if the alignment position, this must be defined
	    if($mod =~ /\-/ && $mod !~ 'NA'){
		#print "$mod \n";
		#sum the individual mfe normalized by the average energy in the species
		$ICI->{$mod}->{val} += $hash->{$org}->{$mod}->{mfe} / $avgMFE ;
		#print "$ICI->{$mod} \n"
	    }
	#}
    }
}

#anz = #spezies - kein homologes - keine rRNA sequenz
#my $sort;
#my $interact;
#sort ICIs
foreach my $key ( keys(%{$ICI}) ) {
    #$key .= "\t";
    #print $key."\n";
    #check how many species have this target
    my $interact = `awk '{print \$1,\$4}' $pred |grep $key\$ | awk '{print \$1}'|sort|uniq |wc -l`;
    #$key =~ s/\t//;
    #print $key;
    chomp($interact);
    #print "$interact \n";
    #if there is at least some conservation
   if($interact > 1){
       my $anz;
       if($aln_start ne ''){
	   #check if there is an rRNA at the position
	   my $lack = &check_rRNA($homs, $key);
	   my $nrLack = scalar(split(/,/,$lack));
     
	   #normalization factor homologe - keine rRNA
	   $anz = $nrHoms - $nrLack;	
	   #print "$homs \n $lack \n $nrHoms: $nrLack : $anz : $interact\n";
       }
       else{
	   $anz = $nrHoms;
       }
	
	$ICI->{$key}->{val} = $ICI->{$key}->{val} / $anz;
	$ICI->{$key}->{cnt} = $interact;

    }
    #keine conservierung
    else{
	$ICI->{$key}->{val} = 0;
	$ICI->{$key}->{cnt} = 1;
    }
    #print "\t $key \t $ICI->{$key}\n";
}

open(OUT,">$out");
#print "$pred : \n";
foreach my $key (sort hashValueDescendingNum (keys(%{$ICI}))) {
    if(	$ICI->{$key} > 0){
	print OUT"\t $tmp[0]_$tmp[1] \t $key \t $ICI->{$key}->{val} \t $ICI->{$key}->{cnt} \n";
    }
}
close(OUT);

sub hashValueDescendingNum {
    $ICI->{$b}->{val} <=> $ICI->{$a}->{val};
}

sub findAllHoms{
    my $fam = shift;
    #$fam = $type eq 'm' ? s/_D(prime)*$// : s/_HP\d$// ;
    #$fam =~ s/_D(prime)*$// if($type eq 'm');
    #$fam =~ s/_HP\d$// if($type eq 'p');
    #such alle spezien in denen die snoRNA vorkaommt im multifasta file, pfad entsprechend snoBoard
    my $grep = `grep $fam $path/$fam.fa | awk -F "_" '{print \$3}' | sort| uniq `;
    $grep =~ s/\n/,/g;
    $grep =~ s/\./_/g;
    #print "HALLO: $fam $grep \n";
    return $grep;
}

sub check_rRNA{
    my ($homs,$mod) = @_;
    #my $homs = shift;
  
    my $lack = '';

    my @spez = split(/,/,$homs);
    my $pos = (split(/\-/,$mod))[1];

    #case-variante mit snRNAs
    my $column = 0;

    given($mod){
	when( /^18S/ ) { $column = 2 }
	when( /^28S/ ) { $column = 4 }
	when( /^5.8S-/ ) { $column = 18 }
	when( /^U1-/ ) { $column = 6 }
	when( /^U2-/ ) { $column = 8 }
	when( /^U4-/ ) { $column = 10 }
	when( /^U4atac-/ ) { $column = 22 }
	when( /^U5-/ ) { $column = 12 }
	when( /^U6-/ ) { $column = 14 }
	when( /^U6atac-/ ) { $column = 24 }
	when( /^U11-/ ) { $column = 20 }
	when( /^U12-/ ) { $column = 16 }
	default { print $mod; exit(0) }
    }


    foreach my $org (@spez){
	$org =~ s/\./_/;
	
	my $awk = "awk '{print \$".$column.",\$".($column+1)."}'";
	
	my $ends = `grep $org $aln_start | $awk `;
	chomp($ends);
	#print $column,"\t",$ends,"\n";

	my @ends = split(/\s/,$ends);
	if(!defined($ends[0]) || $ends[0] eq '-' || $pos < $ends[0] || $pos > $ends[1]){
	    $lack .= $org.",";
		#print "$org has no rRNA sequence at this position\n" 
	}
    }
    return $lack;
}
