#!/usr/bin/perl -w

#use Round;
use lib "/homes/bierdepot/steffi/lib/perl/";
use strict;
use File::Temp qw(tempfile);
use Getopt::Long;
use warnings;
use feature qw/say switch/;

sub usage{
    print "USAGE:\n./ICI_mod.pl  -m|--mod <alignment position of modification> -t|--type [m (methylated)| p (pseudouridylated)] -d|--dir <directory with target predictions> -f|--fasta <path to mfastas> [-s|--sort <sorted list of investigated species>] [-a|-aln <file with start&end-pos of seq in alignment>] [-o|--out <output-file>] [-p|--ps <ps-file>] [-n|--num <number of delivered ICI values default=3>]\n";
    exit;
}

my ($mod,$dir,$out,$type,$path,$sortfile,$ps,$num,$aln_start) = ('','','','','','','','','');
my $help;
my $ICI;

GetOptions('h|help|?' => \$help,
	   'm|mod=s' => \$mod,
	   'd|dir=s' => \$dir,
	   'a|aln=s' => \$aln_start,
	   'o|out=s'  => \$out,
	   'f|fasta=s'  => \$path,
	   's|sort=s'  => \$sortfile,
	   'p|ps=s' => \$ps,
	   'n|num=s' => \$num,
	   't|type=s'  => \$type);

if(defined $help || $mod eq '' || $dir eq '' || ($type ne 'm' && $type ne 'p')){&usage;}

else{
    if($num eq ''){$num = 3;}
    if($sortfile eq ''){$sortfile = "/scr/rioja/steffi/PROJECTS/TARGETS/".(split(/-/,$mod))[0]."-sort_viecher_LONG.lst";}
    if($aln_start eq ''){$aln_start = "/scr/rioja/steffi/PROJECTS/TARGETS/RNAsalsa_new/aln_seq_starts";}
    #my $fams;
    my $viecher;
    #search the position in alignment column (before -energy) either to be psi (HACA-files) or meth (CD-files) und sortieren nach bester energy
    my $cand = $type eq 'p' ?  `egrep  '$mod [[:space:]]{1,}-' $dir/HACA_*|sort -nk 5` :  `egrep  '$mod [[:space:]]{1,}-' $dir/CD_*|sort -nk 5`;
    my @cands = split(/\n/,$cand);

    #print $cand;

    #fuer jede predictete interaction
    foreach my $line(@cands){

	my @line = split(/\s+/,$line);
	my $sno = (split(/\-/,$line[1]))[0]."_".(split(/\-/,$line[1]))[-1];
	$sno =~ s/\s+//;
	#print "$sno end\n";
	$ICI->{$sno}->{ICI} = 0;
	$ICI->{$sno}->{cnt} = 0;
	my $spez = (split(/:/,$line[0]))[1];
	$spez =~ s/\s+$//g;
	#print $spez." end\n";
	#beste energy pro spezies zu jeder familie speicher, da sortiert also nur beim ersten mal
	if(!$viecher->{$spez}->{$sno}){
	    $viecher->{$spez}->{$sno}->{box} = (split(/\-/,$line[1]))[-1];
	    $viecher->{$spez}->{$sno}->{energy} = $line[4];
	    $viecher->{$spez}->{$sno}->{struc} = $line[5];
	    $viecher->{$spez}->{$sno}->{seq} = $line[6];
	    if(defined($viecher->{$spez}->{cnt})){$viecher->{$spez}->{cnt}++;}else{$viecher->{$spez}->{cnt}=1}
	    $viecher->{$spez}->{sum_mfe} += $line[4];
	    #$sum_mfe += $line[4];
        }
    }
    
    
    foreach my $vieh (keys (%{$viecher})){
	#average energy of snoRNAs that bind target in this vieh
	$viecher->{$vieh}->{avg_mfe} = $viecher->{$vieh}->{sum_mfe} / $viecher->{$vieh}->{cnt};
	print "$vieh : $viecher->{$vieh}->{sum_mfe} : $viecher->{$vieh}->{cnt} : $viecher->{$vieh}->{avg_mfe} \n";
	
	foreach my $sno (keys (%{$viecher->{$vieh}})){
	    
	    if($sno ne "sum_mfe" && $sno ne "cnt" && $sno ne "avg_mfe" &&  $sno ne "viecher"){
		#cnt in how many organisms the snoRNA interacts
		$ICI->{$sno}->{cnt}++;
		#sum over all normalized mfes
		$ICI->{$sno}->{viecher} .= $vieh.",";
		$ICI->{$sno}->{ICI} += $viecher->{$vieh}->{$sno}->{energy} / $viecher->{$vieh}->{avg_mfe};
		print "\t $sno : $viecher->{$vieh}->{$sno}->{box} : $viecher->{$vieh}->{$sno}->{energy}  \n";
	    }
	}
    }
    
    foreach my $fam (keys (%{$ICI})){
	#such alle spezien wo die snoRNA vorkommt
	$ICI->{$fam}->{homs} = &findAllHoms($fam);
	$ICI->{$fam}->{lack} = &check_rRNA($ICI->{$fam}->{homs});
#	print "$fam \n $lack \n  $homs\n";
	$ICI->{$fam}->{nrLack} = scalar(split(/,/,$ICI->{$fam}->{lack}));
	$ICI->{$fam}->{nrHoms} = scalar(split(/,/,$ICI->{$fam}->{homs}));
#	print "$ICI->{$fam}->{lack} \t  $ICI->{$fam}->{homs}\n";
	$ICI->{$fam}->{anz} = $ICI->{$fam}->{nrHoms} - $ICI->{$fam}->{nrLack} ;
	#if snoRNA occurs more than once, ICI sum is normalized by number of organisms where snoRNA occurs
	$ICI->{$fam}->{ICI_norm} = $ICI->{$fam}->{cnt} > 1 ? ($ICI->{$fam}->{ICI} / $ICI->{$fam}->{anz}) : 0 ;
	
   }
    #my $sort;
    my $cnt = 0;
    my @print;
    $print[0][0] = "#mod";
    $print[1][0] = $mod;

    foreach my $fam (sort hashValueDescendingNum (keys (%{$ICI})) ) {
	if($ICI->{$fam}->{ICI_norm} > 0 && $cnt < $num){
	    $cnt++;
	    
	    print "--------------\n$fam : $ICI->{$fam}->{ICI} : $ICI->{$fam}->{anz} ($ICI->{$fam}->{nrHoms} - $ICI->{$fam}->{nrLack}): $ICI->{$fam}->{ICI_norm} \n \t HOMOLOGS: $ICI->{$fam}->{homs} \n \t INTERACT: $ICI->{$fam}->{viecher} \n \t NO rRNA: $ICI->{$fam}->{lack} \n---------------\n" ;
	    
	    $print[2][0] .= $fam.",";
	    #$print[3][0] = $fam."\n";
	    
	    open(SORT,'<'.$sortfile);
	    my $cn = 0;
	    while(<SORT>){
		if($_ !~ /^#/){
		$cn++;
		    my @tmp = split(/\t+/,$_);
		    $tmp[1] =~ s/\./_/;
		    chomp($tmp[1]);
		    #print " $tmp[1]\n";
		    $print[0][$cn] .= $tmp[0] if($cnt == 1) ;
		    #$header .= "\t".$tmp[0] if($cnt == 1) ;
		    
		    #no rRNA
		    if($ICI->{$fam}->{lack} =~ /$tmp[1]/){
			#print "no rRNA in $tmp[1]\n";
			$print[1][$cn] = "NA";
			next;
		    }
		    #interaction
		    elsif($ICI->{$fam}->{viecher} =~ /$tmp[1]/){
			#print "interaction in $tmp[1]\n";
			#$line .= "\t".$cnt.":".$viecher->{$tmp[1]}->{$fam}->{energy}.",";
			$viecher->{$tmp[1]}->{$fam}->{energy} =~ s/\s+//;
			$print[1][$cn] .= $cnt.":".$viecher->{$tmp[1]}->{$fam}->{energy}.",";
			$print[2][$cnt].= $tmp[1].",".$viecher->{$tmp[1]}->{$fam}->{struc}."\n";
			$print[3][$cnt].= $viecher->{$tmp[1]}->{$fam}->{seq}."\n";
			next;
		    }
		    #no homologoue
		    elsif($ICI->{$fam}->{homs} !~ /$tmp[1]/){
			#print "no homologue in $tmp[1]\n";
			#$line .= "\t".$cnt.":N,";
			$print[1][$cn] .= $cnt.":N,";
			next;
		    }
		    #no interaction despite homologue and rRNA sequence
		    elsif($ICI->{$fam}->{viecher} !~ /$tmp[1]/ && $ICI->{$fam}->{homs} =~ /$tmp[1]/ && $ICI->{$fam}->{lack} !~ /$tmp[1]/ ){
			#print "no interaction in $tmp[1]\n";
			#$line .= "\t-";
			$print[1][$cn] .= "-" if(!defined($print[1][$cn]) || $print[1][$cn] eq '');
			next;
		    }
		}
	    }
	    close(SORT);
	    $print[0][$cn+1] .= "snoRNA\tICI" if($cnt == 1);
	    #$line .= "\t".$fam."\t".$ICI->{$fam}->{ICI_norm}."\n";
	    $print[1][$cn+1] .= $fam.":".$ICI->{$fam}->{ICI_norm}."\t";
	    
	    #print $header if($cnt == 1);
	    #print $line;
	    #print CSV $header if($cnt == 1);
	    #print CSV $line;
	}
    }
    $print[1][-1] .= ":" if($cnt == 2);
    $print[1][-1] .= ":\t:" if($cnt == 1);

    if($ps ne ''){
	&picture(\@print);
	#if($type eq 'm'){
	#    &color($print[2],$print[3])
	#}
    }
}

sub picture{
    my $print = shift;
    open(CSV,'>'.$mod.'.csv');
    for(my $i = 0; $i < 2; $i++){
	foreach(@{$print->[$i]}){
#	    print "$_ \t";
	    print CSV "$_ \t";
	}
#	print "\n";
	print CSV"\n";
    }
    close(CSV);

    my $rna = (split(/-/,$mod))[0];
    print $rna."\n";
    `/homes/bierdepot/steffi/Diplom/Modifications/Bildchen/tree2psCopies.newICI.pl -t $mod.csv -o $ps -r $rna`;

    print "created ps-file: $ps\n";

}


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

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 "$grep \n";
    return $grep;
}

sub check_rRNA{
    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( /^U1-/ ) { $column = 6 }
	when( /^U2-/ ) { $column = 8 }
	when( /^U4-/ ) { $column = 10 }
	when( /^U5-/ ) { $column = 12 }
	when( /^U6-/ ) { $column = 14 }
	when( /^U12-/ ) { $column = 16 }
	when( /^5.8S-/ ) { return $lack }
    }

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

	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;
}

