#!/usr/bin/perl -w
use strict;
use warnings;
use Cwd;
use List::Util qw(first max maxstr min minstr reduce shuffle sum);
use Pod::Usage;
use Getopt::Long;

### GLOBAL VARIABLES ###
my $clustersWithJust1Structure = 0;
my ($clusterDir , $secondaryStructureDir, $methodFile , $alignmentStructureAgreement , @orderedLabels , $iRMSD , $outputFolder);
my $help                       = 0;

#Cleaning the Cache
#print "removing tcoffee cache..\n";
#my $cacheCleaning = 'rm -rf ~/.t_coffee/cache/*';
#(system "$cacheCleaning") == 0 or die "cannot execute the command $cacheCleaning $!\n";
my $varCheck = `env|grep PDB`;
unless (($varCheck=~/PDB_DIR/) and ($varCheck=~/NO_REMOTE_PDB_DIR=1/)){
   die "ERROR! Before running the script you must export the variables PDB_DIR and NO_REMOTE_PDB_DIR\n";
}


#Parsing parameters
my $result = GetOptions ("c=s"                           => \$clusterDir,
			 "s=s"                           => \$secondaryStructureDir,
			 "m=s"                           => \$methodFile,
			 "o=s"                           => \$outputFolder,
			 "alignmentStructureAgreement=s" => \$alignmentStructureAgreement,
                         "help"                          => \$help);
$help = 1 unless ($result);

if ((! $help) && (! defined $clusterDir)) {
  warn "Error! The cluster directory is missing. Please specify the -c parameter\n";
  $help = 1;
}
if ((! $help) && (! defined $secondaryStructureDir)) {
  warn "Error! The secondary Structure Directory is missing. Please specify the -s parameter\n";
  $help = 1;
}
if ((! $help) && (! defined $methodFile)) {
  warn "Error! The method is missing. Please specify the -m parameter\n";
  $help = 1;
}
if ((! $help) && (! defined $alignmentStructureAgreement)) {
  warn "Error! alignmentStructureAgreement script is missing. Please specify the -alignmentStructureAgreement parameter\n";
  $help = 1;
}
pod2usage(-exitval => 0, -verbose => 2, -noperldoc => 1) if ($help);


if (-d "${outputFolder}/OutDir/"){die "Before going on please remove the OutDir directory\n";}                   #TRICK!!!
my $tempOutDir                 = "${outputFolder}/OutDir/";
my @allClusters = preparing();


#READ THE METHODS
my (%methodsFinalCount , %methodsONclusters , %structureAgreementFinalCount , %structureAgreementONclusters , %averageCoreIndexes , %time);
open (M,"<$methodFile") or die "cannot open the method file $!\n";
foreach my $commandLine (<M>){
  next if (($commandLine=~/^\s*$/) or ($commandLine=~/^#/));
  chomp $commandLine;
  my ($label , $command, $commandRecall);
  if ($commandLine=~/^([^\;]+)\s*\;([^\;]+)$/){
    $command = $commandRecall = $1;
    $label   = $2;
    $label   =~ s/\s//g;
    mkdir $tempOutDir . "/$label";
  }
  else {die "ERROR! label not found in $commandLine !!$!\n";}
  print OUT "$label\n";
  $methodsONclusters{$label} = 1;
  push (@orderedLabels , $label);
  #doShellScript
  open (S,">${outputFolder}/runMethod.sh") or die "cannot create runMethod.sh for\n$command\n$!\n";
  print S '#!/bin/sh'."\n";

  my $clusterNumber = 0;
  foreach my $cluster (@allClusters){
    next if (($cluster eq '.') or ($cluster eq '..') or ($cluster=~/~$/));
    my $command          = $commandRecall;
    my $clusterFile      = $clusterDir . "/$cluster";
    my $fileOut          = $tempOutDir . "/$label/$cluster".".aln";
    my $fileStructureOut = "$tempOutDir/$label/$cluster".".str";
    my $dndOut           = "${outputFolder}" . '/dndThresh' . "/$cluster";
    $command =~ s/<TARGET>/$clusterFile/g;
    $command =~ s/<OUT>/$fileOut/g;
    $command =~ s/<DND>/$dndOut/g;

    #ALIGNING
    print S "$command" . "\n";
  }
  close S;
  (system "/usr/bin/time -o ${outputFolder}/time_$label sh ${outputFolder}/runMethod.sh") == 0 or die "cannot execute /usr/bin/time -o time_$label ${outputFolder}/runMethod.sh\n$!\n";            ##    TRICK!!
  system "rm ${outputFolder}/runMethod.sh"        ;
  system ("rm -rf *rfold 2> /dev/null")           ;
  system ("rm -rf *template_list  2> /dev/null")  ;

 foreach my $cluster (@allClusters){
    next if (($cluster eq '.') or ($cluster eq '..') or ($cluster=~/~$/));
    my $command          = $commandRecall;
    my $clusterFile      = $clusterDir . "/$cluster";
    my $fileOut          = $tempOutDir . "/$label/$cluster".".aln";
    my $fileStructureOut = "$tempOutDir/$label/$cluster".".str";
    my $dndOut           = 'dndThresh' . "/$cluster";

    system "rm $clusterDir/${cluster}.dnd" if (-f $clusterDir . "/${cluster}.dnd");   #clean clustalw tree
    open (STRUCTURE,">$fileStructureOut") or die "cannot create the $fileStructureOut file $!\n";

    #taking the cluster members SEQUENCES
    my (%clusters , %pairs , %gaps);
    open (F,"<$tempOutDir/$label/$cluster".".aln") or die "cannot open $tempOutDir/$label/$cluster for taking the sequences names $!\n";
    foreach my $line (<F>){
      chomp $line;
      next if (($line=~/CLUSTAL/) or ($line=~/PROBCONS/) or ($line=~/Probalign/) or ($line =~/^\s*$/) or ($line=~/^\s+/));
      if ($line=~/^(\S+)\s+(\S+)/){
	$clusters{$1} .= $2;
      }
      else {die "cannot parse the $tempOutDir/$label/$cluster line $line \n";}
    }
    close F;

    $clusterNumber++;


    #mapping the GAPS
    foreach my $seq (keys %clusters){
      foreach my $index (1..length ($clusters{$seq})){
	my $symbol = (substr $clusters{$seq},$index-1,1);
	push (@{$gaps{$seq}}, $index)if ($symbol eq '-');
      }
    }
    #taking members secondary structure PAIRS
    foreach my $seq (keys %clusters){
      if ($seq=~/^([^\_]+)\_(\S+)/){
	my $seqId = $1;
	my $chain = $2;
	next unless (-f "$secondaryStructureDir"."/$seqId".'.pdb.'."$chain".'.pairs');
	open (S,"$secondaryStructureDir"."/$seqId".'.pdb.'."$chain".'.pairs') or die "cannot open the $seq secondary structure file $!\n";
	foreach my $line (<S>){

	  if ($line =~ /\*\s+(\d+)\s+(\d+)/){
	    my %hash = ('L' => $1, 'R' => $2);
	    push (@{$pairs{$seq}}, \%hash);
	  }
	}
	close S;
      }
      else {die "cannot recognize the FASTA name $seq\n";}
    }
    #sanity check: skipp when the cluster has not sequence with a proper secondary structure
    my $check = 0;
    foreach my $seq (keys %clusters){
     $check = 1 if (defined @{$pairs{$seq}});
    }
    if ($check == 0){
      print OUT "There are not secondary structures kaksi fixed for the cluster $cluster , skipping\n";
      next;
    }
    #update the paris according to the gaps
    my $realStructureInTheCluster = 0;
    foreach my $seq (keys %clusters){
      print STRUCTURE "#$seq\n";
      $realStructureInTheCluster++ if (scalar @{$pairs{$seq}} > 0);
      foreach my $oldPair (@{$pairs{$seq}}){
	foreach my $gap (@{$gaps{$seq}}){
	  if ($gap <= "$oldPair->{L}"){
	    $oldPair->{'L'}++;
	  }

	  if ($gap <= "$oldPair->{R}"){
	    $oldPair->{'R'}++;
	  }
	}
	print STRUCTURE "$oldPair->{L} $oldPair->{R}\n";
      }
    }

    #STRUCTURE AGREEMENT BENCHMARCK
    my ($averageStructureAgreement , @score);
    if ($realStructureInTheCluster > 1){
      @score = `$alignmentStructureAgreement  -a $fileOut -m $fileStructureOut`;
      $averageStructureAgreement = $score[0];
      chomp $averageStructureAgreement;
      $structureAgreementFinalCount{$label}          +=  $averageStructureAgreement;
      $structureAgreementONclusters{$label}{$cluster} =  $averageStructureAgreement;
    }
    else {$clustersWithJust1Structure++;}


    #Printing output for each cluster
    print OUT "$cluster";
    if (defined $averageStructureAgreement){
      print OUT ";$averageStructureAgreement";
    }
    else {print OUT ";just1Structure";}
    print OUT "\n";
    close STRUCTURE;
  }
  $averageCoreIndexes{$label} /= $clusterNumber if (%averageCoreIndexes);
}
















#OUTPUTS
print OUT "\n\n#STRUCTURAL AGREEMENT\n";
foreach my $label (keys %structureAgreementFinalCount){
  print OUT "3SP score for $label is ";
  my $s = $structureAgreementFinalCount{$label} / 34;
  $s = sprintf("%.2f", $s);
  print OUT "$s\n";
}
#initializing the hashes
my (%winningMethod_structure , %evenRecall_structure);
foreach my $label (keys %structureAgreementONclusters){
  $winningMethod_structure{$label} = 0;
  $evenRecall_structure{$label}    = 0;
}
#Computing the average performance for the methods and printing them
my $noAgreement = 0;
foreach my $label (keys %structureAgreementONclusters){
  foreach my $cluster (keys %{$structureAgreementONclusters{$label}}){
    my $score   = $structureAgreementONclusters{$label}{$cluster};
    my $evenSpy = 0;
    my %evenTest;

    foreach my $otherLabel (keys %structureAgreementFinalCount){
      next if ($otherLabel eq $label);
      $score = max ($score , $structureAgreementONclusters{$otherLabel}{$cluster});
    }
    if ($score == 0){$noAgreement++;next;}
    #looks for even winners
    foreach my $otherLabel (keys %structureAgreementFinalCount){
      if ($score == $structureAgreementONclusters{$otherLabel}{$cluster}){
	$winningMethod_structure{$otherLabel}++;
	$evenSpy++;
	$evenTest{$otherLabel} = 1;
      }
    }
    if ($evenSpy > 1){
      foreach my $evenLabel (keys %evenTest){
	$evenRecall_structure{$evenLabel}++;
	$winningMethod_structure{$evenLabel}--;
      }
    }
  }
  last;
}
print OUT "\n";
foreach my $label (keys %winningMethod_structure){
  print OUT "$label is better in producing structure agreement in $winningMethod_structure{$label} clusters\n";
}
print OUT "\n";
foreach my $label (keys %evenRecall_structure){
  print OUT "$label is even in producing structure agreement in $evenRecall_structure{$label} clusters\n";
}
my $methodNumber = keys %methodsONclusters;
$clustersWithJust1Structure = $clustersWithJust1Structure / $methodNumber;
print OUT "$clustersWithJust1Structure clusters have just 1 real structure. No possible to compute the structure agreement here\n";
if ($noAgreement){
  print OUT "\nscore equal 0 for $noAgreement clusters. Probable it has been impossible predict any structure\n";
}


#PAIRWISE COMPARISON MATRIX STRUCTURAL CONSISTENCY
print OUT "\n##MATRIX STRUCTURAL CONSISTENCY\n";
foreach my $label (@orderedLabels){print OUT "$label ";}
print OUT "\n";
foreach my $label (@orderedLabels){
  print OUT "$label";
  foreach my $otherLabel (@orderedLabels){
    if ($otherLabel eq $label){print OUT "  -"; next;}
    my $paiwiseComparison = 0;
    foreach my $cluster (keys %{$structureAgreementONclusters{$label}}){
      my $score   = $structureAgreementONclusters{$label}{$cluster};
      my $otherScore = $structureAgreementONclusters{$otherLabel}{$cluster};
      $paiwiseComparison++ if ($otherScore > $score);       #just > and not >= so that we exlude the evens
    }
    if ($paiwiseComparison > 9){
      print OUT " $paiwiseComparison";
    }
    else{
      print OUT "  $paiwiseComparison";
    }
  }
  print OUT "\n";
}
close OUT;





system ("rm *.rfold 2> /dev/null");
system ("rm -rf ${outputFolder}/dndThresh 2> /dev/null");
system ("rm *.template_list 2> /dev/null");
convert2StandardClustal();


#FUNCTIONS
sub preparing {
  mkdir $tempOutDir;
  mkdir "${outputFolder}/dndThresh" unless (-d "${outputFolder}/dndThresh");
  open (OUT,">${outputFolder}/3SP.out") or die "cannot create the 3SP output file\n";
  #READ THE CLUSTER DIRECTORY
  opendir (C,$clusterDir) or die "cannot open the cluster dir $!\n";
  my @allClusters;
  while (my $c = readdir C){
	next if (($c eq '.') or ($c eq '..') or ($c=~/~$/));
  	push (@allClusters,$c);
  }
  closedir C;
  return @allClusters;
}

sub convert2StandardClustal {
  foreach my $dir (@orderedLabels){
    opendir (DIR , "${tempOutDir}/$dir") or die "cannot open ${tempOutDir}/$dir \n$!\n";
    while (my $file = readdir(DIR)){
      if ($file=~/.aln$/){
	open (N ,">${tempOutDir}/${dir}/${file}.new") or die "cannot open ${tempOutDir}/${dir}/${file}.new $!\n";
	open (O ,"<${tempOutDir}/${dir}/${file}") or die "cannot open ${tempOutDir}/${dir}/${file} $!\n";
	foreach my $l (<O>){
	  $l=~s/PROBCONS/CLUSTAL/g;
	  $l=~s/Probalign/CLUSTAL/g;
	  print N $l;
	}
	close N;
	close O;
	(system "mv ${tempOutDir}/${dir}/${file}.new ${tempOutDir}/${dir}/${file}") == 0 or die "cannot mv ${tempOutDir}/${dir}/${file}.new ${tempOutDir}/${dir}/${file} \n$!\n";;
      }
    }
    closedir DIR;
  }
}







=head1 NAME

runSaraBenchmark_V10.pl - This script run the benchmarck of SaraTcoffee. The V9 was working perfectly. This update is just meant to clean up the code to run just the needed analysis. Moreover it calls the iRMSD scripts

=head1 SYNOPSIS

perl runSaraBenchmark_V10.pl  -c clusterDir -s secondaryStructureDir -m methodFile -alignmentStructureAgreement fileName -help

=head1 DESCRIPTION

=over

=item * This script run different MSA methods on the input clusters, compute the consensus secondary structure,

=item * compare this prediction with the real secondary structures of the sequencences composing the alignment which has been computed with rnass.pl (from emidio Capriotti)and then refined with kaksiSecondaryStructurePDBfix.pl. This is a measure of the "modelling capacity" of a method.

=item * Each identical secondary structure pair match score a point. The total score for each cluster is normalized dividing by the number of real structure pairs

=item * If the method file you specify contains "-output clustalw,html" the script will also produce the .html files. These are parsed in order to obtain the core index score which will be reported in the output.

=item * if "alignmentStructureAgreement" you must specify the alignmentStructureAgreement.py script path. This script creates for each cluster a color cache file in the output folder. That kind of file can be used to produce colored alignment (in this way: t_coffee -other_pg seq_reformat -in=./15.rcoffee.aln -in2 ./rcoffee_secoStruc -input2 alifold -action +add_alifold -output stockholm_aln -struc_in=./15.rcoffee.aln.color_cache -output=color_html -out=test.html)

=item * Using this option the script will also return the structure agreement benchmark scores, indicating which method is the best in producing alignments that support the real secondary structure alignment. this is the first score returned by alignmentStructureAgreement.py. It is the sum of the square of the ratio between number of occurrencies of a certain pair in the cluster and the number of sequences, all divided by the number of different pairs.

=item * This option also return a file called "pairwise_strucConsistency_matrix" and "pairwise_modellingCapacity_matrix". These matrices contains a pairwise method comparison of the methods here listed and sorted as follows: "Tcoffee_vanila","Rcoffee","BestPair","RcoffeeReal","ConsanCoffee","SaraCoffee".

=item * It print the number of times in which the structural consistency score of the METHOD ON THE COLUMNS is equal to the one on the ROWS

=item * By default the modelling capacity score for each cluster is divided by the number of real structures pairs in that cluster. Use the flag "notWeigthed" if you don't want it.

=back

=head1 OPTIONS

=over

=item B<-c> I<Directory>

Directory containing the cluster files

=item B<-s> I<Directory>

Directory containing the secondary structure files produced with kaksiSecondaryStructurePDBfix.pl

=item B<-m> I<fileName>

file Containing the compared methods. An example is shown as a comment in this script.

=item B<-alignmentStructureAgreement> I<fileName>

specify the alignmentStructureAgreement.py fileName tu run the structure agreement benchmark, produce the color_cache files and the pairwise method comparison matrix

Help

=cut



#METHOD FILE

#~/Desktop/Apps/tCoffee-subversion/trunk/t_coffee/src/t_coffee -in <TARGET> -method sara_pair -template_file=TEMPLATE_FILE,RNA -extend_mode rna2 -outfile <OUT> -output clustalw -newtree <DND>  ; SaraCoffee
#t_coffee <TARGET> -outfile <OUT> -output clustalw -newtree <DND> > /dev/null 2>/dev/null;Tcoffee_vanila
#t_coffee -in <TARGET> -mode rcoffee -outfile <OUT> -output clustalw -newtree <DND>  > /dev/null 2>/dev/null ;Rcoffee
