#!/usr/bin/perl
# ID = pid_timestamp
# Something wrong when ID = fasta_file.pid, if fasta_file use absolute path.
use warnings;
use strict;
use Time::HiRes qw(gettimeofday);

########################################
if(@ARGV != 7){
  die "Need seven arguments: 
           (1) protein sequence fasta file as input ; 
           (2) profile file as the first output; 
           (3) protein contact map (DCA) file as the second output;
           (4) path to DCA_exe;
           (5) path to hhsuite;
           (6) path to hhsuite protein sequence database;
           (7) path to temporary work folder.";
}
my $fasta_file = $ARGV[0];
my $profile_file = $ARGV[1];
my $output_file = $ARGV[2];
my $dir1 = $ARGV[3];
my $dir2 = $ARGV[4];
my $protein_db_hhblits = $ARGV[5];
my $dir3 = $ARGV[6];

my $pl_fas2msa = "$dir1/fas2msa.pl";
my $exe_DCA = "$dir1/DCA_exe";
my $exe_hhblits = "$dir2/bin/hhblits";
my $pl_reformat = "$dir2/scripts/reformat.pl";

my $pid = $$;
my $timestamp = int (gettimeofday * 1000);
my $id = "${pid}_$timestamp";
my $num_CPUs = 3; # you can change it

####################################
# run HHblits to generate a MSA file
my $a3m_file = "$dir3/$id.a3m";
my $fas_file = "$dir3/$id.fas"; 
my $hhblits_output_file = "$dir3/$id.hhblits.out";
my $msa_file = "$dir3/$id.msa"; 
`$exe_hhblits -i $fasta_file -d $protein_db_hhblits -oa3m $a3m_file -cpu $num_CPUs -n 3 -o $hhblits_output_file`;

`perl $pl_reformat $a3m_file $fas_file`;
`perl $pl_fas2msa $fas_file > $msa_file`;

###############################
# check if MSA is less than 10
open IN, "$msa_file" or die $!;
my $num_lines = 0;
while(<IN>){ $num_lines++; }
close IN;
if($num_lines < 8){
  open IN, "$hhblits_output_file" or die $!;
  my $n = 0;
  my %hash;
  while(my $line = <IN>){
    chomp $line;
    if(substr($line,0,7) eq "Probab="){
      $n++;
      my @items = split(/\s+/, $line);
      $items[0] =~ s/Probab\=//;
      $items[1] =~ s/E-value\=//;
      push(@{$hash{$n}}, $items[0]);
      push(@{$hash{$n}}, $items[1]);
    }
  }
  close IN;
  my $ep;
  for my $key (sort {$a<=>$b} keys %hash){
    if($hash{$key}->[0] > 50){
      $ep = $hash{$key}->[1];
    }
  }
  `$exe_hhblits -i $fasta_file -d $protein_db_hhblits -oa3m $a3m_file -cpu $num_CPUs -e $ep -n 3 -o $hhblits_output_file`;
  `perl $pl_reformat $a3m_file $fas_file`;
  `perl $pl_fas2msa $fas_file > $msa_file`;
}


#####################################################
# run DCA to generate the contact map prediction file
`$exe_DCA $msa_file $profile_file $output_file`;

`rm $hhblits_output_file`;
`rm $a3m_file`;
`rm $fas_file`;
`rm $msa_file`;

# End



