#!/usr/bin/perl -w
# Bruno Contreras-Moreira, Pablo Vinuesa, ...
# 2005-8 UNAM/Mexico, EEAD/CSIC

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

my $progname = "marfil_amps.pl";
my $VERSION = 2; # Oct 2cnd, 2011, upgraded by PV to run with phyml_3.0.1_linux64; 
                 # getPHYML3parameters, run_PHYML3, get_PUZZLE_param_from_PHYML3_stat, run_PUZZLE_LM in phyTools were adapted accordingly 
my $REMOVETMPFILES;
my @files2remove;

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

my(%opts,$DNAfile,$skip_ModelTest,$label_tree);
my($search,$optimizations,$user_tree_file,$random_trees,$no_of_random_trees);
getopts('hcln:m:N:O:S:U:R', \%opts);

if(($opts{'h'})||(scalar(keys(%opts))==0))
{
        print "usage: $progname [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 "-n \t input file with nucleotide sequences .fna\n";
        print "-m \t skip ModelTest, run model provided by user (optional,[GTRG|HKYG|TrNIG|...])\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 "-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 "-l \t label tree branches with FASTA headers    (optional)\n";
        exit;
}

if(defined($opts{'n'})){ $DNAfile = $opts{'n'}; }
else{ die "# $progname : need a valid nucleotide FNA file\n" }
if(defined($opts{'c'})){ $REMOVETMPFILES = 0; }else{ $REMOVETMPFILES = 1 }
if(defined($opts{'m'})){ $skip_ModelTest = $opts{'m'} }
else{ $skip_ModelTest = '' }
if(defined($opts{'S'})){ $search = $opts{'S'}; }else{ $search = 'BEST'; }
if(defined($opts{'O'})){ $optimizations = $opts{'O'}; }else{ $optimizations = 'tlr';} 
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 }

print "# $progname : -n $DNAfile -m $skip_ModelTest -S $search -O $optimizations -U $user_tree_file -R $random_trees -N $no_of_random_trees -l $label_tree -c $REMOVETMPFILES\n";

######################################################
 
# 1) convert DNA alignment to phylip and nexus formats
my $DNA_phylip_file = convert_FAS2PHY($DNAfile,$REMOVETMPFILES);
if($DNA_phylip_file eq $ERROR)
{
	die "# $progname : sorry, FAS2PHY failed ($DNAfile)\n";
}

push(@files2remove,$DNA_phylip_file);

my $DNA_nexus_file = convert_FAS2NXS($DNAfile,$REMOVETMPFILES);
if($DNA_nexus_file eq $ERROR)
{
	unlink(@files2remove);
	die "# $progname : sorry, FAS2NXS failed ($DNAfile)\n";
}

push(@files2remove,$DNA_nexus_file);

my ($MODELTEST_bestAICmodel,$MODELTEST_file);
unless($skip_ModelTest)
{
	# 2) run MODELTEST3.7 using the PAUP output of -lnL scores generated by running MODELBLOCK3
	($MODELTEST_bestAICmodel,$MODELTEST_file) = run_MODELTEST($DNA_nexus_file,$REMOVETMPFILES);
	if($MODELTEST_bestAICmodel eq $ERROR)
	{
		unlink(@files2remove);
		die "# $progname : sorry, run_MODELTEST failed ($DNAfile)\n";
	}
	#system("rm modelfit.log model.scores"); # cleanup
}
else
{
	# 2.1) use model provided by user
	$MODELTEST_bestAICmodel = $skip_ModelTest;
	$MODELTEST_file = '-';
}

# 3) select PHYML parameters for this model
my $PHYML3params = getPHYML3parameters($MODELTEST_bestAICmodel);
if($PHYML3params eq $ERROR)
{
	unlink(@files2remove);
	die "# $progname : sorry, getPHYML3parameters failed with model $MODELTEST_bestAICmodel ($DNAfile)\n";
}

print "# The PHYML3params:\n$PHYML3params\n";

# 4) run PHYML with this parameters
#$phylip_file,$PHYMLparams,$clean,$search,$optimizations,$user_tree_file,$random_trees,$no_of_random_trees
my ($PHYML_tree,$PHYML_stat,$mean_aLRT,$median_aLRT) = run_PHYML3($DNA_phylip_file,$PHYML3params,$REMOVETMPFILES,$search,$optimizations,$user_tree_file,$random_trees,$no_of_random_trees);
#print "# got from run_PHYML3: [$PHYML_tree,$PHYML_stat,$mean_aLRT,$median_aLRT]\n"; 

if($PHYML_tree eq $ERROR)
{
	unlink(@files2remove);
	die "# $progname : sorry, runPHYML3 failed ($DNAfile)\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($DNAfile);	
		$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_tree,$PHYML_stat);

# 5) parse PHYML stats in order to create params file to run PUZZLETREE
# print "# $progname is runnning: get_PUZZLE_param_from_PHYML3_stat($PHYML_stat)\n";
my $PUZZLE_param_file = get_PUZZLE_param_from_PHYML3_stat($PHYML_stat);
if($PUZZLE_param_file eq $ERROR)
{
	unlink(@files2remove);
	die "# $progname : sorry, get_PUZZLE_param_from_PHYML3_stat failed ($DNAfile)\n";
}

#push(@files2remove,$PUZZLE_param_file);

# 6) run PUZZLE TREE for Likelihood Mapping (LM) with ML-param-file returned from get_PUZZLE_param_from_PHYML_stat()
my ($PUZZLE_dist,$PUZZLE_LM_puzzle,$alpha,$aln_stats,$distance_stats) = run_PUZZLE_LM($DNA_phylip_file,$PUZZLE_param_file, $REMOVETMPFILES);
if($PUZZLE_LM_puzzle eq $ERROR)
{
	unlink(@files2remove);
	die "# $progname : sorry, runPUZZLE_LM failed ($DNAfile)\n";
}

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

# 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 ($DNAfile)\n";
}

print "# output : subs_model = $MODELTEST_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 : $DNA_phylip_file $PHYML_tree $PUZZLE_param_file $MODELTEST_file $PUZZLE_dist $PUZZLE_LM_puzzle $PHYML_stat $DNAfile\n";

if($REMOVETMPFILES)
{
	unlink($DNA_nexus_file);
}
