#!/usr/bin/perl -w
# Pablo Vinuesa, Bruno Contreras-Moreira...
# 2011 UNAM, Mexico

# marfil_prot.pl v2, based on marfil_prot.pl; extended options set, to control prottest2.4 and phyml_3.0.1_linux64; also muscle3.8.31 instead of the older binaries and code
# New subs and $ENV{VARs} were added to phyTools.pm

use strict;
use Getopt::Std;
use lib "lib";
use phyTools;
set_phyTools_env();

my $progname = "marfil_prot.pl";
my $VERSION = 2; # upgraded by PV, Sept 2011, to be used with prottest2.4 and phyml_3.0.1_linux64; also muscle3.8.31
my $REMOVETMPFILES;
my @files2remove;

######################################################

my $do_not_copy_input_file = 0;
my ($AAfile,$myAAfile,$exclude_file,$seq,$skip_muscle,$skip_ProtTest,$label_tree,%opts);
my ($search , $random_trees,  $user_tree_file, $Puzzle_models,$optimizations,$no_of_random_trees,);
 # params=tlr : tree topology (t), branch length (l) and rate parameters (r) are optimised.
 # params=tl  : tree topology and branch length are optimised.
 # params=lr  : branch length and rate parameters are optimised.
 # params=l   : branch length are optimised.
 # params=r   : rate parameters are optimised.
 # params=n   : no parameter is optimised.

getopts('hlci:e:O:S:U:DsPRm:n:', \%opts);

if(defined($opts{'h'})||(scalar(keys(%opts))==0)) 
{ 
	print "usage: $progname v.$VERSION [options]\n";
	print "-h \t this message\n";
	print "-c \t do not remove temporary files             (optional; by def will remove temp files)\n";
	print "-i \t input file with protein sequences .faa    (REQUIRED)\n";
	print "-e \t exclude file                              (optional)\n";
	print "-s \t skip making alignment, input is aligned   (optional)\n";
	print "-S \t Search algorithm                          (optional, [NNI|SPR|BEST]; def = BEST)\n";
	print "-O \t Optimizations                             (optional, [tlr|tl|lr|l|r|n], def = tlr)\n";
	print "-P \t use subset of Puzzle-compatible matrices  (optional, def = 0)\n";
	print "-R \t Random trees                              (optional,  def = '')\n";
	print "-U \t User-tree file name                       (optional, [filename], def = '')\n";
	print "-n \t number of random start trees              (optional, [integer], def = 5)\n";
	print "-m \t skip ProtTest, run model provided by user (optional,[WAGG|JTTG|...])\n";
	print "-D \t do not copy input file                    (optional, ignored with -e)\n";
	print "-l \t label tree branches with FASTA headers    (optional)\n\n";
	exit; 
}

if(defined($opts{'i'})){ $AAfile = $opts{'i'}; }
else{ die "# $0 : need a valid protein FASTA file\n" }
if(defined($opts{'c'})){ $REMOVETMPFILES = 0; }else{ $REMOVETMPFILES = 1 }
if(defined($opts{'e'})){ $exclude_file = $opts{'e'}; }else{ $exclude_file = ''; }
if(defined($opts{'D'} && $exclude_file eq '')){ $do_not_copy_input_file= 1; }
if(defined($opts{'s'})){ $skip_muscle = 1 }else{ $skip_muscle = 0; }
if(defined($opts{'S'})){ $search = $opts{'S'}; }else{ $search = 'BEST'; }
if(defined($opts{'O'})){ $optimizations = $opts{'O'}; }else{ $optimizations = 'tlr';} 
if(defined($opts{'P'})){ $Puzzle_models = 1; }else{ $Puzzle_models = 0;} 
if(defined($opts{'R'})){ $random_trees = 1 }else{ $random_trees = '' }
if(defined($opts{'U'})){ $user_tree_file = $opts{'U'}; }else{ $user_tree_file =''; }
if(defined($opts{'n'})){ $no_of_random_trees = $opts{'n'}; }else{ $no_of_random_trees = 5; }
if(defined($opts{'l'})){ $label_tree = 1 }else{ $label_tree = 0; }
if(defined($opts{'m'})){ $skip_ProtTest = $opts{'m'}; }
else{ $skip_ProtTest = '' }

print "# $progname : -i $AAfile -e $exclude_file -D $do_not_copy_input_file -s $skip_muscle -S $search -O $optimizations -U $user_tree_file -R $random_trees -n $no_of_random_trees -P $Puzzle_models -m $skip_ProtTest -l $label_tree -c $REMOVETMPFILES\n";
#exit;
######################################################

if($do_not_copy_input_file == 0)
{
	# -2) read exclusion file that contains species that should be removed from the analysis
	my @EXCLIST;
	if($exclude_file)
	{
		open(EXCLIST,$exclude_file) || die "#progname : cannot read $exclude_file\n";
		while(<EXCLIST>)
		{
			next if(/^#/);
			push(@EXCLIST,(split)[0]);
		}
		close(EXCLIST);
	}

	# -1) add sequence number to original FASTA sequences
	my %AAfasta = read_FASTA_sequence($AAfile);

	$AAfile = reverse $AAfile;
	$AAfile = (split(/\//,$AAfile))[0];
	$AAfile = reverse $AAfile;

	$myAAfile = "my_" . (split(/\./,$AAfile))[0] . ".faa";

	open(MYAA,">$myAAfile") || die "#progname : cannot write to $myAAfile\n";

	print "\n# FASTA header equivalences (AA):\n";

	foreach $seq (sort {$a<=>$b} (keys(%AAfasta)))
	{
		if($AAfasta{$seq}{'NAME'} =~ /\[(\S+)\]/){ next if(grep(/$1/,@EXCLIST)); }	
		
		$AAfasta{$seq}{'NAME'} =~ s/>//;
		$AAfasta{$seq}{'SEQ'} =~ s/\*//;
		
		print "# $seq => $AAfasta{$seq}{'NAME'}";
		
		print MYAA ">$seq " . $AAfasta{$seq}{'NAME'};
		print MYAA $AAfasta{$seq}{'SEQ'},"\n";
	}

	close(MYAA);
	print "\n";

	push(@files2remove,$myAAfile);
}
else{ $myAAfile = $AAfile }


# 0) run muscle to obtain multiple alignment of protein sequences
my $AA_aln_file;
if($skip_muscle){ $AA_aln_file = $myAAfile }
else
{
	$AA_aln_file = run_muscle3831($myAAfile,$REMOVETMPFILES); # changed from $myAAfile
	if($AA_aln_file eq $ERROR)
	{
		unlink(@files2remove);
       		die "# $progname : sorry, muscle alignment failed ($AAfile)\n";
	}
}
#push(@files2remove,$myAAfile);
 
# 1) convert DNA alignment to phylip and nexus formats
my $AA_phylip_file = convert_FAS2PHY($AA_aln_file,$REMOVETMPFILES);
if($AA_phylip_file eq $ERROR)
{
	unlink(@files2remove);
	die "# $progname : sorry, FAS2PHY failed ($AA_aln_file)\n";
}

push(@files2remove,$AA_phylip_file);

my ($PROTTEST_bestAICmodel,$PROTTEST_outfile,$PHYMLparams);
unless ($skip_ProtTest)
{       
	# 2) run PROTTEST and get best fit model
	($PROTTEST_bestAICmodel,$PROTTEST_outfile) = run_PROTTEST2($AA_phylip_file,$REMOVETMPFILES,$Puzzle_models);
	if($PROTTEST_bestAICmodel eq $ERROR)
	{
		unlink(@files2remove);
		die "# $progname : sorry, run_PROTTEST2 failed ($AAfile)\n";
	}#print "|$PROTTEST_bestAICmodel|\n";#exit; 
}
else
{
	# 2.1) use model provided by user
	$PROTTEST_bestAICmodel = $skip_ProtTest;
	$PROTTEST_outfile = '-';
}

# 3.1) select PHYML3 parameters for this model
$PHYMLparams = getPHYML3_parameters_prot($PROTTEST_bestAICmodel);
if($PHYMLparams eq $ERROR)
{
	unlink(@files2remove);
	die "# $progname : sorry, getPHYML3parameters_prot failed with model $PROTTEST_bestAICmodel ($AAfile)\n";
} #print "$PHYMLparams\n";


## 4) run PHYML3 with this parameters

#my ($PHYML_tree,$PHYML_lk,$PHYML_stat,$mean_aLRT,$median_aLRT) = run_PHYML($AA_phylip_file,$PHYMLparams,$REMOVETMPFILES);
my ($PHYML_tree,$PHYML_stat,$mean_aLRT,$median_aLRT) = run_PHYML3($AA_phylip_file,$PHYMLparams,$REMOVETMPFILES,$search,$optimizations,$user_tree_file,$random_trees,$no_of_random_trees);
#print "# $progname: run_PHYML3 returned: $PHYML_tree,$PHYML_stat,$mean_aLRT,$median_aLRT \n";

if($PHYML_tree eq $ERROR)
{
	unlink(@files2remove);
	die "# $progname : sorry, runPHYML3 failed\n";
}
elsif($label_tree)
{
        my ($newick,$labelled_newick) = ('','');

        open(TREE,$PHYML_tree) || die "# $0 : cannot read $PHYML_tree\n";
        while(<TREE>){ $newick .= $_; }
        close(TREE);

        if($newick)
        {
                my %ampFASTA = read_FASTA_sequence($AAfile);
                $labelled_newick = add_labels2newick_tree( $newick, \%ampFASTA );
                open(LABELTREE,">$PHYML_tree") || die "#$0 : cannot rewrite $PHYML_tree\n";
                print LABELTREE $labelled_newick;
                close(LABELTREE); #print "#$newick#$labelled_newick#\n";
        }
        else{ print "# $0 : cannot label tree $PHYML_tree\n"; }
}

#push(@files2remove,$PHYML_stat); # $PHYML_tree

## 5) create a file with PUZZLE_LM paramerers for protein alignments
#print "# $progname running: create_PUZZLE_LM_prot_parameters_file($AA_phylip_file,$PHYML_stat,$PROTTEST_bestAICmodel);\n";
my $PUZZLE_param_file = create_PUZZLE_LM_prot_parameters_file($AA_phylip_file,$PHYML_stat,$PROTTEST_bestAICmodel);
if($PUZZLE_param_file eq $ERROR)
{
	unlink(@files2remove);
	die "# $progname : sorry, create_PUZZLE_LM_prot_parameters_file failed ($AA_phylip_file)\n";
}

push(@files2remove,$PUZZLE_param_file);

## 6) run PUZZLE TREE for Likelihood Mapping (LM) with ML-param-file returned from create_PUZZLE_LM_prot_parameters_file()

if($Puzzle_models) # run puzzle only if modeltest or prottest have been set to run only the subset of puzze-compatible models
{
    #print "# $progname running:run_PUZZLE_LM($AA_phylip_file,$PUZZLE_param_file, $REMOVETMPFILES)\n";
    my ($PUZZLE_dist,$PUZZLE_LM_puzzle,$alpha,$aln_stats,$distance_stats) = run_PUZZLE_LM($AA_phylip_file,$PUZZLE_param_file, $REMOVETMPFILES);
    #print "# run_PUZZLE_LM returned: $PUZZLE_dist,$PUZZLE_LM_puzzle,$alpha,$aln_stats,$distance_stats\n";

    if($PUZZLE_LM_puzzle eq $ERROR)
    {
    	unlink(@files2remove);
	die "# $progname : sorry, runPUZZLE_LM failed ($AA_phylip_file)\n";
    }

    push(@files2remove,$PUZZLE_dist,); # $PUZZLE_LM_puzzle


    if($REMOVETMPFILES){ unlink @files2remove; }

    ## 7) parse LM_PUZZLE file and print the corresponding MySQL.input table
    my $DB_row = parse_LM_PUZZLE_file($PUZZLE_LM_puzzle);
    if($DB_row eq $ERROR)
    {
	unlink(@files2remove);
	die "# $progname : sorry, parse_LM_PUZZLE_file failed ($AAfile)\n";
    } 

    print "# output : subs_model = $PROTTEST_bestAICmodel stats = $DB_row\n";
    print "# alpha = $alpha alignment stats : $aln_stats\n";
    print "# distance stats : $distance_stats\n";
    print "# mean aLRT = $mean_aLRT median aLRT = $median_aLRT\n";
    print "# outfiles : $AA_phylip_file $PHYML_tree $PUZZLE_param_file $PROTTEST_outfile $PUZZLE_dist $PUZZLE_LM_puzzle $PHYML_stat $AA_aln_file\n";
}

#open ( OUT, ">>Table_puzzle_LM_DNA_alphas.input" ) || die "cannot write Table_puzzle_LM_DNA_alphas.input";
#print OUT "$DB_row\n";
#close(OUT);

# 8) Upload data into MySQL DB
#if(!upload_LM_data2DB(%LM_puzzle_data))
#{
#  die "# $progname : sorry, upload_LM_data2DB failed\n";
#} 
