# Bruno Contreras-Moreira, Pablo Vinuesa
# 2005-12 CCG/UNAM, Mexico, EEAD/CSIC, Zaragoza, Spain
# This is a collection of subroutines used in our projects,
# including primers4clades and get_homologues.pl


package phyTools;
require Exporter;

@ISA = qw(Exporter);
@EXPORT = qw(
set_phyTools_env $ERROR read_FASTA_sequence feature_is_installed check_installed_features
extract_gene_name extract_CDSs_from_genbank extract_intergenic_from_genbank extract_features_from_genbank extract_genes_from_genbank
find_taxa_FASTA_headers find_taxa_FASTA_array_headers read_FASTA_file_array NAME SEQ 
run_muscle run_muscle3831 run_stable_py run_mafft_einsi run_mablock run_codehop convert_FAS2PHY convert_FAS2NXS run_MODELTEST
run_MrModeltest get_DNA_AA_seq_pairs check_aligned_FASTA_sequences convert_AAaln2NTaln
getPHYMLparameters getPHYML3parameters getPHYMLparameters_prot getPHYML3_parameters_prot getPHYMLparameters_prot_100boot run_PHYML run_PHYML3
get_PUZZLE_param_from_PHYML_stat get_PUZZLE_param_from_PHYML3_stat run_PUZZLE_LM run_PUZZLE_DIST run_PUZZLE_DIST_HKYG run_NEIGHBOR
run_codeml get_mean_boot4WAG_dist_tree get_mean_aLRT4WAGG_tree run_PHYML_WAG_G parse_LM_PUZZLE_file
select_amplicons select_intergene_amplicons calc_amplicon_length_aa amplicon_too_similar
cut_amplicon_aln_files cut_intergen_amplicon_aln_files run_PROTTEST run_PROTTEST2 create_PUZZLE_LM_prot_parameters_file
put_fullname_in_topology filter_aln_rate4site get_GO_terms4gi translate_GO_terms find_exon_boundaries_FASTA_headers
find_introns_FASTA_headers adjust_introns_alignment translate_dna2prot convert_FASTA_dna2prot run_BRANCHCLUST
cluster_sequences select_codon_tables calc_correlation_coef calc_oligo_degeneracy rank_codon_usage_tables
check_primer_match check_primer_match_excluded extract_codon_freqs_from_fasta2spsumfile add_labels2newick_tree 
find_redundant_taxa_FASTA_CDHIT run_CDHIT 
plot_pangenome_heatmap
);

$phyTools::ERROR = "phyError";
$phyTools::NERROR = "-999";

use strict;
use File::Basename;
use Cwd;

# bioperl installation path if not standard
use lib "/home/contrera/soft/marfil/CVS_code/bioperl-1.5.2_102/";

eval{ require Bio::Seq } || die "phyTools: failed to import Bio::Seq, please install Bioperl (please check manual)\n@_";
eval{ require Bio::SeqIO } || die "phyTools: failed to import Bio::SeqIO, please install Bioperl (please check manual)\n@_";

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

my $MARFILPATH = "/path/to/marfil/";

$phyTools::ERROR = "phyError";
$phyTools::NERROR = "-999";

# constant as indices to access read_FASTA_file_ arrays
use constant NAME => 0; 
use constant SEQ  => 1;
use constant IDENTICALS => 2; 
use constant NUMBER_IDENTICALS  => 3;

my %feature_output = ( 
# default output of binaries when corrected installed, add as required for 
# check_installed_features
'EXE_BLASTP'=>'BLAST','EXE_BLASTN'=>'BLAST','EXE_FORMATDB'=>'database',
'EXE_MCL'=>'Usage:','EXE_MULTIPARANOID'=>'args','EXE_R'=>'R is free software',
'EXE_MAKEHASH'=>'t open output file','EXE_READBLAST'=>'Error opening','EXE_COGLSE'=>'Usage:',
'EXE_HMMPFAM'=>'Usage:','EXE_COGTRI'=>'Error opening', #'EXE_COGTRI'=>'Runs',
'EXE_INPARA'=>'usage','EXE_ORTHO'=>'usage','EXE_HOMOL'=>'usage','EXE_SPLITBLAST'=>'Usage',
'EXE_SPLITHMMPFAM'=>'Usage'
);


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

set_phyTools_env();
#foreach my $key (keys %ENV)
#{
#    print "# $key => $ENV{$key}\n";
#
#} exit;

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

sub set_phyTools_env
{
   if( ! defined($ENV{'MARFIL_MISSING_BINARIES'}) ) { $ENV{'MARFIL_MISSING_BINARIES'} =  ''; }
   if( ! defined($ENV{'MARFIL'}) ) { $ENV{'MARFIL'} =  $MARFILPATH .'/'; }
   if( ! defined($ENV{'BIN'}) ) { $ENV{'BIN'} =  $ENV{'MARFIL'} . 'bin/'; }	   
   if( ! defined($ENV{'TMP_DIR'}) ){ $ENV{'TMP_DIR'} = "/tmp/" } #data/TempData/marfil/" }

   if( ! defined($ENV{'BLIMPS_DIR'}) ) {$ENV{'BLIMPS_DIR'} =  $ENV{'BIN'} . "/BLOCKS3.7_software/blimps3.7/blimps-3.7/";}
   if( ! defined($ENV{'BLIMPS_BIN'}) ) {$ENV{'BLIMPS_BIN'} =  "$ENV{'BLIMPS_DIR'}/bin/linux/";}
   if( ! defined($ENV{'BLIMPS_DOCS'}) ) {$ENV{'BLIMPS_DOCS'} = "$ENV{'BLIMPS_DIR'}/docs/";}
   if( ! defined($ENV{"FILE_MODELBLOCK"}) ) { $ENV{"FILE_MODELBLOCK"} = $ENV{'BIN'} . "/mymodelblock3"; }
   if( ! defined($ENV{"GB_SPSUM_DIR"}) ) { $ENV{"GB_SPSUM_DIR"} = $ENV{'BIN'} . "/CUD/"; }
   if( ! defined($ENV{"GBBCT_SPSUM"}) ) { $ENV{"GBBCT_SPSUM"} = $ENV{"GB_SPSUM_DIR"} . "/gbbct.spsum"; }
   
   # $ENV{"GB_UNIVERSAL_CODON_TABLES"} generated with bin/CUD/get_universal_tables.pl
   if( ! defined($ENV{"GB_UNIVERSAL_CODON_TABLES"}) ) { $ENV{"GB_UNIVERSAL_CODON_TABLES"} = "10090.mitochondrion:mitochondrion Mus musculus,gbrod.spsum,10090:Mus musculus,gbrod.spsum,102126:Pseudanabaena sp. PCC 6903,gbbct.spsum,10360:Human herpesvirus 5,gbvrl.spsum,10710:Enterobacteria phage lambda,gbphg.spsum,11676:Human immunodeficiency virus 1,gbvrl.spsum,10566:Human papillomavirus,gbvrl.spsum,118965:Bacillus subtilis,gbbct.spsum,188937:Methanosarcina acetivorans C2A,gbbct.spsum,224911:Bradyrhizobium japonicum USDA 110,gbbct.spsum,239758:SARS coronavirus Taiwan TC3,gbvrl.spsum,11908:Human T-lymphotropic virus,gbvrl.spsum,273057:Sulfolobus solfataricus P2,gbbct.spsum,311553.chloroplast:chloroplast Oryza sativa,gbpln.spsum,311553.mitochondrion:mitochondrion Oryza sativa,gbpln.spsum,311553:Oryza sativa,gbpln.spsum,32608:Human adenovirus type 2,gbvrl.spsum,3702.chloroplast:chloroplast Arabidopsis thaliana,gbpln.spsum,3702.mitochondrion:mitochondrion Arabidopsis thaliana,gbpln.spsum,3702:Arabidopsis thaliana,gbpln.spsum,5270:Ustilago maydis,gbpln.spsum,273507:Phanerochaete chrysosporium,gbpln.spsum,57667:Simian-Human immunodeficiency virus,gbvrl.spsum,6239.mitochondrion:mitochondrion Caenorhabditis elegans,gbinv.spsum,6239:Caenorhabditis elegans,gbinv.spsum,7227.mitochondrion:mitochondrion Drosophila melanogaster,gbinv.spsum,7227:Drosophila melanogaster,gbinv.spsum,83331:Mycobacterium tuberculosis CDC1551,gbbct.spsum,83333:Escherichia coli K12,gbbct.spsum,208524.mitochondrion:mitochondrion Gallus gallus,gbvrt.spsum,208526:Gallus gallus,gbvrt.spsum,9606.mitochondrion:mitochondrion Homo sapiens,gbpri.spsum,9606:Homo sapiens,gbpri.spsum,9615.mitochondrion:mitochondrion Canis familiaris,gbmam.spsum,9615:Canis familiaris,gbmam.spsum,9913.mitochondrion:mitochondrion Bos taurus,gbmam.spsum,297284:Bos taurus,gbmam.spsum" } 
   
   #if( ! defined($ENV{"DBSERVER"}) ) { $ENV{"DBSERVER"} = "132.248.34.219"; }
   #if( ! defined($ENV{"DBNAME"}) ) { $ENV{"DBNAME"} = "marfil"; }
   #if( ! defined($ENV{"DBARCH"}) ) { $ENV{"DBARCH"} = "mysql"; }
   #if( ! defined($ENV{"DBUSER"}) ) { $ENV{"DBUSER"} = "marfil"; }
   #if( ! defined($ENV{"DBPASS"}) ) { $ENV{"DBPASS"} = "m4rf1l"; } # maya pswd m4rf1lm4y4
   #if( ! defined($ENV{"BLASTNDBPATH"}) ) { $ENV{"BLASTNDBPATH"} = $ENV{'MARFIL'}."moreno_blastnDB/"; }
   #if( ! defined($ENV{"IMG_GENOMES_FASTA_DIR"}) ) { $ENV{"IMG_GENOMES_FASTA_DIR"} = $ENV{'MARFIL'}."IMG_GENOMES/"; } 
   if( ! defined($ENV{"GO_DIR"}) ) { $ENV{"GO_DIR"} = $ENV{'MARFIL'}."GO/GO_ANNOTATIONS/"; }
   if( ! defined($ENV{"GO_DICTIONARY"}) ) { $ENV{"GO_DICTIONARY"} = $ENV{'MARFIL'}."GO/GO.terms_and_ids"; }   
   if( ! defined($ENV{"PUZZLE_PROT_DIST_PARAMS_FILE"}) ) { $ENV{"PUZZLE_PROT_DIST_PARAMS_FILE"} = $ENV{'MARFIL'}."CVS_code/marfil_data/puzzle_prot_dist_params.txt"; } # auto model  
   if( ! defined($ENV{"PUZZLE_PROT_WAG_DIST_PARAMS_FILE"}) ) { $ENV{"PUZZLE_PROT_WAG_DIST_PARAMS_FILE"} = $ENV{'MARFIL'}."CVS_code/marfil_data/puzzle_prot_WAG_dist_params.txt"; } # hemos visto que WAG es el mas chingon
   if( ! defined($ENV{'PUZZLE_DNA_HKY_DIST_PARAMS_FILE'}) ) { $ENV{'PUZZLE_DNA_HKY_DIST_PARAMS_FILE'} = $ENV{'MARFIL'}."CVS_code/marfil_data/puzzle_dna_HKY_dist_params.txt"; } 
   if( ! defined($ENV{"NEIGHBOR_PARAMS_FILE"}) ) { $ENV{"NEIGHBOR_PARAMS_FILE"} = $ENV{'MARFIL'}."CVS_code/marfil_data/NEIGHBOR_params.txt"; }   
   if( ! defined($ENV{"CODEML_CTRL_FILES_DIR"}) ) { $ENV{"CODEML_CTRL_FILES_DIR"} = $ENV{'MARFIL'}."CVS_code/marfil_data/"; }   
   if( ! defined($ENV{"PFAMDB"}) ) { $ENV{"PFAMDB"} = $ENV{'MARFIL'}."db/Pfam-A.hmm"; } # HMMER3
 
   #wget ftp://ftp.kazusa.or.jp/pub/codon/current/gbbct.spsum # current GB release version GB151 of Jan06
   # if( ! defined($ENV{"FILE_MrMODELBLOCK"}) ) { $ENV{"FILE_MrMODELBLOCK"} = "/home/vinuesa/Software_downloads/stow/MrModeltest2.2/MrModelblock"; }
   
   ## generic executables and directories
   # if( ! defined($ENV{"EXE_PYTHON"}) ){ $ENV{"EXE_PYTHON"} = "/usr/bin/python"; }
   if( ! defined($ENV{"EXE_MUSCLE"}) ){ $ENV{"EXE_MUSCLE"} = $ENV{'MARFIL'}."bin/muscle"; }
   if( ! defined($ENV{"EXE_MUSCLE3831"}) ){ $ENV{"EXE_MUSCLE3831"} = $ENV{'MARFIL'}."bin/muscle3831"; }            # <=== PV Sep. 28th 2011
   if( ! defined($ENV{'EXE_STABLEPY'}) ){ $ENV{'EXE_STABLEPY'} = $ENV{'MARFIL'}."bin/stable_muscle.py"; }          # <=== PV Sep. 28th 2011
   if( ! defined($ENV{'EXE_MAFFT_EINSI'}) ){ $ENV{'EXE_MAFFT_EINSI'} = $ENV{'MARFIL'}.'bin/mafft-6.814-without-extensions/binaries/einsi'; } # pablo Jul 15th, 2010   
   if( ! defined($ENV{"EXE_CLUSTAL"}) ){ $ENV{"EXE_CLUSTAL"} = $ENV{'MARFIL'}."bin/clustalw1.83/clustalw"; }
   if( ! defined($ENV{"EXE_CLUSTALO"}) ){ $ENV{"EXE_CLUSTALO"} = $ENV{'MARFIL'}."bin/clustalo-1.2.0-x86_64"; }

   if( ! defined($ENV{"EXE_MABLOCK"}) ){ $ENV{"EXE_MABLOCK"} =  "$ENV{'BLIMPS_BIN'}/mablock"; }
   if( ! defined($ENV{"EXE_MAKE_SPSUM"}) ){ $ENV{"EXE_MAKE_SPSUM"} = $ENV{'MARFIL'}."CVS_code/make_spsum_codon_file.pl"; }
   #if( ! defined($ENV{"EXE_CODUSE"}) ){ $ENV{"EXE_CODUSE"} = "$blimps_bin/coduse"; }
   if( ! defined($ENV{"EXE_CODEHOP"}) ){ $ENV{"EXE_CODEHOP"} = "$ENV{'BLIMPS_BIN'}/codehop"; }
   if( ! defined($ENV{"EXE_PHYML"}) ){ $ENV{"EXE_PHYML"} = $ENV{'MARFIL'}."bin/phyml_alrt_linux"; } # changed to alart_linux 32 bit on 18th June 07 PV; originally phyml_linux
   if( ! defined($ENV{"EXE_PHYML3"}) ){ $ENV{"EXE_PHYML3"} = $ENV{'MARFIL'}."bin/PhyML-3.1_linux64"; } # PV May2014
   #if( ! defined($ENV{"EXE_PHYML"}) ){ $ENV{"EXE_PHYML"} = '/usr/local/biotools/bin/phyml_3.0.1_linux64'; }#
   #if( ! defined($ENV{"EXE_PAUP"}) ){ $ENV{"EXE_PAUP"} = $ENV{'MARFIL'}.bin/paup }
   if( ! defined($ENV{"EXE_PAUP"}) ){ $ENV{"EXE_PAUP"} = $ENV{'MARFIL'}."/bin/paup"; }
   if( ! defined($ENV{"EXE_PUZZLE"}) ){ $ENV{"EXE_PUZZLE"} = $ENV{'MARFIL'}."bin/puzzle"; } #
   if( ! defined($ENV{"EXE_MODELTEST"}) ){ $ENV{"EXE_MODELTEST"} = $ENV{'MARFIL'}."bin/modeltest3.7"; }
   #if( ! defined($ENV{"EXE_PROTTEST"}) ){ $ENV{"EXE_PROTTEST"} = $ENV{'MARFIL'}."bin/ProtTest1.3/runProtTest"; }
   if( ! defined($ENV{"EXE_PROTTEST2"}) ){ $ENV{"EXE_PROTTEST2"} = $ENV{'MARFIL'}."bin/ProtTest2.4/runProtTest"; } # <=== PV Sep. 28th 2011
   if( ! defined($ENV{"EXE_CODEML"}) ){ $ENV{"EXE_CODEML"} = $ENV{'MARFIL'}."bin/paml3.15/bin/codeml"; }
   if( ! defined($ENV{"PROTTEST_DIR"}) ){ $ENV{"PROTTEST_DIR"} = $ENV{'MARFIL'}."bin/ProtTest1.3/"; }
   if( ! defined($ENV{"PROTTEST2_DIR"}) ){ $ENV{"PROTTEST2_DIR"} = $ENV{'MARFIL'}."bin/ProtTest2.4/"; }            # <=== PV Sep. 28th 2011                
   if( ! defined($ENV{"EXE_MAT2PNG"}) ){ $ENV{"EXE_MAT2PNG"} = $ENV{'MARFIL'}."bin/mat2png/matrix2png"; }
   if( ! defined($ENV{"EXE_AU_MAKERMT"}) ){ $ENV{"EXE_AU_MAKERMT"} = $ENV{'MARFIL'}."bin/consel01i/bin/makermt"; }
   if( ! defined($ENV{"EXE_AU_CONSEL"}) ){ $ENV{"EXE_AU_CONSEL"} = $ENV{'MARFIL'}."bin/consel01i/bin/consel"; }
   if( ! defined($ENV{"EXE_AU_CATPV"}) ){ $ENV{"EXE_AU_CATPV"} = $ENV{'MARFIL'}."bin/consel01i/bin/catpv"; }

   if( ! defined($ENV{"EXE_AU_MAKERMT2"}) ){ $ENV{"EXE_AU_MAKERMT2"} = $ENV{'MARFIL'}."bin/consel020/bin/makermt"; }
   if( ! defined($ENV{"EXE_AU_CONSEL2"}) ){ $ENV{"EXE_AU_CONSEL2"} = $ENV{'MARFIL'}."bin/consel020/bin/consel"; }
   if( ! defined($ENV{"EXE_AU_CATPV2"}) ){ $ENV{"EXE_AU_CATPV2"} = $ENV{'MARFIL'}."bin/consel020/bin/catpv"; }

   if( ! defined($ENV{"EXE_CHECKPRIMERS"}) ){ $ENV{"EXE_CHECKPRIMERS"} = $ENV{'MARFIL'}."CVS_code/check_primers.py"; }
   
   # blast+
if( ! defined($ENV{"BLAST_PATH"}) ){ $ENV{"BLAST_PATH"} = $ENV{'MARFIL'}."bin/ncbi-blast-2.2.24+/bin/"; }
if( ! defined($ENV{"BLAST_PATH"}) ){ $ENV{"BLAST_PATH"} = $ENV{'MARFIL'}."bin/ncbi-blast-2.2.24+/bin/"; }
   if( ! defined($ENV{"EXE_BLASTP"}) ){ $ENV{"EXE_BLASTP"} = $ENV{"BLAST_PATH"}."blastp"; }
   if( ! defined($ENV{"EXE_BLASTN"}) ){ $ENV{"EXE_BLASTN"} = $ENV{"BLAST_PATH"}."blastn"; }
   if( ! defined($ENV{"EXE_FORMATDB"}) ){ $ENV{"EXE_FORMATDB"} = $ENV{"BLAST_PATH"}."makeblastdb"; }
 
   # old blast <== reactivated and edited by PV June 4th, 2011
   #if( ! defined($ENV{"BLAST_PATH_OLD"}) ){ $ENV{"BLAST_PATH_OLD"} = $ENV{'MARFIL'}."bin/blast-2.2.18/bin/"; }
   #$ENV{"EXE_BLASTP"} = $ENV{"BLAST_PATH_OLD"}."blastall "; 
   #$ENV{"EXE_FORMATDB"} = $ENV{"BLAST_PATH_OLD"}."formatdb "; 
   #$ENV{"EXE_FASTACMD"} = $ENV{"BLAST_PATH_OLD"}."fastacmd "; 
   
   if( ! defined($ENV{"EXE_SPLITBLAST"}) ){ $ENV{"EXE_SPLITBLAST"} = $ENV{'MARFIL'}."forge/_split_blast.pl"; }
   if( ! defined($ENV{"EXE_SPLITHMMPFAM"}) ){ $ENV{"EXE_SPLITHMMPFAM"} = $ENV{'MARFIL'}."forge/_split_hmmscan.pl"; }
 
   if( ! defined($ENV{"EXE_RATE4SITE"}) ){ $ENV{"EXE_RATE4SITE"} = $ENV{'MARFIL'}."bin/rate4site.exe"; }
   if( ! defined($ENV{"EXE_NEIGHBOR"}) ){ $ENV{"EXE_NEIGHBOR"} = $ENV{'MARFIL'}."bin/phylip3.66/exe/neighbor_marfil"; }
   if( ! defined($ENV{"EXE_SEQBOOT"}) ){ $ENV{"EXE_SEQBOOT"} = $ENV{'MARFIL'}."bin/phylip3.66/exe/seqboot"; }
   if( ! defined($ENV{"EXE_CONSENSE"}) ){ $ENV{"EXE_CONSENSE"} = $ENV{'MARFIL'}."bin/phylip3.66/exe/consense"; }
   if( ! defined($ENV{"EXE_CODEML"}) ){ $ENV{"EXE_CODEML"} = $ENV{'MARFIL'}."bin/paml3.15/bin/codeml"; }
   if( ! defined($ENV{"EXE_MARFILFACT"}) ){ $ENV{"EXE_MARFILFACT"} = $ENV{'MARFIL'}."bin/AutoFACT/scripts/marfilFACT.pl"; }
   if( ! defined($ENV{"EXE_MCL"}) ){ $ENV{"EXE_MCL"} = $ENV{'MARFIL'}."/bin/mcl-02-063/shmcl/mcl"; }
   if( ! defined($ENV{"EXE_MULTIPARANOID"}) ){ $ENV{"EXE_MULTIPARANOID"} = $ENV{'MARFIL'}."/bin/marfil_multiparanoid.pl"; }
   if( ! defined($ENV{"EXE_INPARA"}) ){ $ENV{"EXE_INPARA"} = $ENV{'MARFIL'}."forge/_cluster_makeInparalog.pl"; }
   if( ! defined($ENV{"EXE_ORTHO"}) ){ $ENV{"EXE_ORTHO"} = $ENV{'MARFIL'}."forge/_cluster_makeOrtholog.pl"; }
   if( ! defined($ENV{"EXE_HOMOL"}) ){ $ENV{"EXE_HOMOL"} = $ENV{'MARFIL'}."forge/_cluster_makeHomolog.pl"; }
   if( ! defined($ENV{"EXE_CDHIT"}) ){ $ENV{"EXE_CDHIT"} = $ENV{'MARFIL'}."/bin/cd-hit/cd-hit "; }
   if( ! defined($ENV{"EXE_HMMPFAM"}) ){ $ENV{"EXE_HMMPFAM"} = $ENV{'MARFIL'}."/bin/hmmer-3.0/src/hmmscan64 --noali --acc --cut_ga "; } # HMMER3
   if( ! defined($ENV{"EXE_MAKEHASH"}) ){ $ENV{"EXE_MAKEHASH"} = $ENV{'MARFIL'}."/bin/COGsoft/COGmakehash/COGmakehash "; }
   if( ! defined($ENV{"EXE_READBLAST"}) ){ $ENV{"EXE_READBLAST"} = $ENV{'MARFIL'}."/bin/COGsoft/COGreadblast/COGreadblast "; }
   if( ! defined($ENV{"EXE_COGLSE"}) ){ $ENV{"EXE_COGLSE"} = $ENV{'MARFIL'}."/bin/COGsoft/COGlse/COGlse "; }
   if( ! defined($ENV{"EXE_COGTRI"}) ){ $ENV{"EXE_COGTRI"} = $ENV{'MARFIL'}."/bin/COGsoft/COGtriangles/COGtriangles "; }
   #if( ! defined($ENV{"EXE_COGTRI"}) ){ $ENV{"EXE_COGTRI"} = "perl ".$ENV{'MARFIL'}."/bin/COGsoft/COGtriangles/COGtriangles.reformat.pl "; }

   if( ! defined($ENV{"EXE_R"}) ){ $ENV{"EXE_R"} = "R " } # expected to be in the path
 
   ## default parameters also required in primer4clades web server
   if( ! defined($ENV{"DEF_CODEHOP_T"}) ){ $ENV{"DEF_CODEHOP_T"} = 35; } # C centigrades default temp for 
   if( ! defined($ENV{"DEF_QUALITY_PHYLO_EVAL"}) ){ $ENV{"DEF_QUALITY_PHYLO_EVAL"} = 70; } # required for 

   ## user-defined parameters
   if( ! defined($ENV{"GBKFEATURES"}) ){ $ENV{"GBKFEATURES"} = 'rRNA,tRNA'; } # used by extract_features_from_genbank
   if( ! defined($ENV{"MAXPRIMERHEADER"}) ){ $ENV{"MAXPRIMERHEADER"} = 100; } # 25 5May09; max length for header in cut_amplicons primer blocks (option -s) 
    #$ENV{'MARFIL'}bin/BLOCKS3.7_software/blimps3.7/blimps-3.7/examples/codehop/my_codehop.csh
   if( ! defined($ENV{"MAX_PRIMER_TEMP_DIFF"}) ){ $ENV{"MAX_PRIMER_TEMP_DIFF"} = 3; } # C centigrades Codehop, 2, 3
   if( ! defined($ENV{"MAX_PRIMER_SCORE_DIFF"}) ){ $ENV{"MAX_PRIMER_SCORE_DIFF"} = 4.0; } # 3.0, 4.0 codehops score, 6.0
   if( ! defined($ENV{"MAX_PRIMER_DEGEN_DIFF"}) ){ $ENV{"MAX_PRIMER_DEGEN_DIFF"} = 2; } #2, 3 degeneracy in log(2) scale 
   if( ! defined($ENV{"MAX_PRIMER_LENGTH_DIFF"}) ){ $ENV{"MAX_PRIMER_LENGTH_DIFF"} = 3; } # 3, 4 nt difference in total primer length  
   if( ! defined($ENV{"MAX_AMPLICON_LENGTH"}) ){ $ENV{"MAX_AMPLICON_LENGTH"} = 1000; } # nt en cut amplicons
   if( ! defined($ENV{"MIN_AMPLICON_LENGTH"}) ){ $ENV{"MIN_AMPLICON_LENGTH"} = 450; } # nt en cut amplicons
   if( ! defined($ENV{"MAX_AMPLICON_COVER"}) ){ $ENV{"MAX_AMPLICON_COVER"} = 0.1; } # 0.2 proportion of coverage when comparing 2 amplicons of same CUT
   if( ! defined($ENV{"MIN_PROP_PRIMER_MATCHES"}) ){ $ENV{"MIN_PROP_PRIMER_MATCHES"} = 0.80; } #0.80, 0.85, 0.65 proportion of matches over the primer length when deciding whether a degenerate primers matches well its genomic target (only to determine primer/amplicon coordinates on the genome in BLAST analysis)
   if( ! defined($ENV{"MIN_CONGR_PVALUE"}) ){ $ENV{"MIN_CONGR_PVALUE"} = 0.05; }
   if( ! defined($ENV{"HIGH_CONGR_PVALUE"}) ){ $ENV{"HIGH_CONGR_PVALUE"} = 0.25; } 
   if( ! defined($ENV{"MIN_FREQ_CONSENSUS"}) ){ $ENV{"MIN_FREQ_CONSENSUS"} = 0.7; } 
   #if( ! defined($ENV{"MEGABLAST_EVAL"}) ){ $ENV{"MEGABLAST_EVAL"} = 0.001; }
   if( ! defined($ENV{"GENOMIC_MAP_UNIT"}) ){ $ENV{"GENOMIC_MAP_UNIT"} = 10000; } #25000 pb or 25 kb
   #if( ! defined($ENV{"BRANCHCLUST_FRACTION"}) ){ $ENV{"BRANCHCLUST_FRACTION"} = 0.8; } #proportion of taxa
   #if( ! defined($ENV{"BRANCHCLUST_SIZE"}) ){ $ENV{"BRANCHCLUST_SIZE"} = 10; } #number of taxa
   # not used anymore if( ! defined($ENV{"NRTABLECHISQP"}) ){ $ENV{"NRTABLECHISQP"} = 0.05; } # p=0.05,0.01,0.001 
   if( ! defined($ENV{"MIN_PAIR_DIST_REDUND"}) ){ $ENV{"MIN_PAIR_DIST_REDUND"} = 0.005; } # cluster_sequences
   if( ! defined($ENV{"MAX_PAIR_DIST"}) ){ $ENV{"MAX_PAIR_DIST"} = 2.21; } # cluster_sequences(),Fasta2clusters.pl: 2.2 seems to be ok for cellullases
   if( ! defined($ENV{"MAX_AMP_MARGIN_DEV"}) ){ $ENV{"MAX_AMP_MARGIN_DEV"} = 50; } # it will move +- $ENV{"MAX_AMP_MARGIN_DEV amino positions in order to check for similar amplicons
   if( ! defined($ENV{"MAX_INTAMP_MARGIN_DEV"}) ){ $ENV{"MAX_INTAMP_MARGIN_DEV"} = 30; } # same as above, for intergenic amplicons, in nts
   if( ! defined($ENV{"ABS_MAX_SITE_RATE"}) ){ $ENV{"ABS_MAX_SITE_RATE"} = 3.0; } 
   # $ENV{"ABS_MAX_SITE_RATE"}  absolute max. site rates for AA alns as calculated by rate4site.exe in the output r4sOrig.res, values should range between 2.5 and 3.5

   # user-defined parameters for evaluate_primers.pl, relative to primer quality
   if( ! defined($ENV{"MAX_DEGEN"}) ){ $ENV{"MAX_DEGEN"} = 17; } # 13 # can be made more relaxed 9, 13, 19, 33, 49 65 or 129 or strengthened 13, 17
   if( ! defined($ENV{"MAX_RELAXED_DEGEN"}) ){ $ENV{"MAX_RELAXED_DEGEN"} = 25; } # 25 # can be 17, 25, 33, 65, 73, 121
   if( ! defined($ENV{"MAX_FULLDEGEN"}) ){ $ENV{"MAX_FULLDEGEN"} = 129;} # 129  # can be stringent 49, 65, 129, 193, 257 or relaxed 289, 513 577, 850, 1153, 1853, 4609
   if( ! defined($ENV{"MAX_HAIRPINPOT"}) ){ $ENV{"MAX_HAIRPINPOT"} = 0.58; } # 0.58, all 0.51 or 0.56 values can be relaxec to 0.61 or vice vs.
   if( ! defined($ENV{"MAX_SELFPOT"}) ){ $ENV{"MAX_SELFPOT"} = 0.58 } # 0.58, 0.61
   if( ! defined($ENV{"MAX_F_R_CROSSPOT"}) ){ $ENV{"MAX_F_R_CROSSPOT"} = 0.56; } # 0.56
   if( ! defined($ENV{"MAX_NO_OVERLAPS"}) ){ $ENV{"MAX_NO_OVERLAPS"} = 10; }
   if( ! defined($ENV{"PHYML_WAG_BOOTSTRAP_FACTOR"}) ){ $ENV{"PHYML_WAG_BOOTSTRAP_FACTOR"} = 100; }
   
   ## user-defined DB parameters that refer to columns in the amplicons table of the DB
   if( ! defined($ENV{"MAXCHISQFAIL"}) ){ $ENV{"MAXCHISQFAIL"} = 1 } # edit as necessary (1 is min)
   if( ! defined($ENV{"MINPERFRQ"}) ){ $ENV{"MINPERFRQ"} = 80.0; } # edit as necessary	

   ## server parameters
   if( ! defined($ENV{"MAXPHYLOEVALS"}) ){ $ENV{"MAXPHYLOEVALS"} = 5; } # Fasta2primers -l
   if( ! defined($ENV{"MAXRUNTIME"}) ){ $ENV{"MAXRUNTIME"} = 400; } # seconds, Fasta2primers -
   if( ! defined($ENV{"TRAILINGZEROS"}) ){ $ENV{"TRAILINGZEROS"} = 3; } # trailing zeros for sequence identifiers, produced by read_FASTA_sequences, (3 for max 999 sequences)
}

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

sub check_installed_features
{
	# check all needed binaries and data sources required by functions here
        # fills $ENV{"MISSING_BINARIES"} with missing binaries
	my (@to_be_checked) = @_;

        my ($check_summary,$output) = ("\nChecking required binaries and data sources, all set in phyTools.pm :\n");

        foreach my $bin (@to_be_checked)
        {
                $check_summary .= sprintf("%18s : ",$bin);
                if($ENV{$bin})
		{ 
			if('EXE_MAFFT_EINSI' eq $bin){ $output = `$ENV{$bin} -h 2>&1 ` }
			elsif('EXE_R' eq $bin){ $output = `$ENV{$bin} --version 2>&1 ` }
			else{ $output = `$ENV{$bin} 2>&1 ` }
			if(!$output){ $output = '' }
		}
                
		if($output =~ /$feature_output{$bin}/)
                {
                        if($bin eq 'EXE_HMMPFAM') # special check for hammer (Pfam hmm library)
                        {
                                if(!-s $ENV{'PFAMDB'})
                                {
                                        $output = "Pfam-A.hmm file needs to be installed (ftp://ftp.sanger.ac.uk/pub/databases/Pfam/)\n";
                                        $ENV{"MARFIL_MISSING_BINARIES"} .= "$bin,";
                                }
                                elsif(!-s $ENV{'PFAMDB'}.'.h3m')
                                {
                                        $output = "Pfam-A.hmm file should be prepared with hmmpress (same folder as $ENV{'EXE_HMMPFAM'})\n";
                                        $ENV{"MARFIL_MISSING_BINARIES"} .= "$bin,";
                                }
                                else{ $output = "OK ($ENV{$bin} $ENV{'PFAMDB'})" }
                        }
                        else{ $output = "OK (path:$ENV{$bin})" }
                }
		else
                {
                        $ENV{"MARFIL_MISSING_BINARIES"} .= "$bin,";
                        $output = " needs to be installed (wrong path:$ENV{$bin})";
                }
                $check_summary .= $output."\n";
        }#die $ENV{"MARFIL_MISSING_BINARIES"};

        return $check_summary;
}

sub feature_is_installed
{
	# checks whether a passed software can be used before calling it
        my ($feature) = @_;

        my $env_missing = $ENV{"MARFIL_MISSING_BINARIES"}; 
        if($feature eq 'BLAST')
        {
                if($env_missing =~ /BLASTP/ || $env_missing =~ /BLASTN/ || $env_missing =~ /FORMATDB/){ return 0 }
        }
        elsif($feature eq 'PFAM')
        {
                if($env_missing =~ /PFAM/){ return 0 }
        }
        elsif($feature eq 'OMCL')
        {
                if($env_missing =~ /MCL/){ return 0 }
        }
        elsif($feature eq 'PRND')
        {
                if($env_missing =~ /PARANOID/){ return 0 }
        }
        elsif($feature eq 'COGS')
        {
                if($env_missing =~ /HASH/ || $env_missing =~ /COG/ || $env_missing =~ /READBLAST/){ return 0 }
        }
	elsif($feature eq 'MAFFT')
        {
                if($env_missing =~ /MAFFT/){ return 0 }
        }
	elsif($feature eq 'R')
        {
                if($env_missing =~ /EXE_R/){ return 0 }
        }

        return 1;
}

# needs a pangenome_matrix.tab file as single argument to run
# Pablo May2011
sub plot_pangenome_heatmap
{
    my ($pangenome_matrix_file) = @_;
    my ($matrix_basename) = (split /\.tab/, $pangenome_matrix_file)[0];
    my $pangenome_heatmap_output_png = $matrix_basename . '_heatmap.png';
    my $pangenome_heatmap_output_pdf = $matrix_basename . '_heatmap.pdf';

    open (R, '|R --no-save') or die "can't open R: $!\n";
    print R<<EOF;

    require(graphics); require(grDevices)
    table <- read.table("$pangenome_matrix_file", header=TRUE)
    x <- as.matrix(table)
    rc <- rainbow(nrow(x), start=0, end=.3)
    cc <- rainbow(ncol(x), start=0, end=.3)
    pdf(file="$pangenome_heatmap_output_pdf")
    # hv <- heatmap(x, col = cm.colors(256), scale="column", RowSideColors = rc, ColSideColors = cc, margins=c(5,10), xlab = "genomes", ylab="genes", main = "heatmap(<Mtcars data>, ..., scale = \"column\")")
    # Error: evaluaciÃ³n anidada demasiado profunda; recursiÃ³n infinita options(expressions= )?
    #heatmap(x, Colv = NA, col = cm.colors(256), scale="column",RowSideColors = rc, margins=c(8,20),  xlab = "genes", ylab="genomes", main ="heatmap genes presence/absence")
    heatmap(x, Colv = NA,  scale="column", margins=c(8,18),  xlab = "gene families", ylab="genomes", main ="Pangenome heatmap")
    png(file="$pangenome_heatmap_output_png")
    heatmap(x, Colv = NA,  scale="column", margins=c(8,18),  xlab = "gene families", ylab="genomes", main ="Pangenome heatmap")
    dev.off()
    q()

EOF
    close R;
}


sub run_muscle
{
   my($infile,$clean) = @_;

    my @split = split(/\//,$infile); # mejor usar File::Basename
    my $basefile = $split[$#split];
    my @file_name_parts;
    @file_name_parts = split(/\./,$basefile);
	
	
   my $tmpfile1 = (split(/\./,$basefile))[0] . "_tmp1.mus";
   my $tmpfile2 = (split(/\./,$basefile))[0] . "_tmp2.mus";
   my $outfile  = (split(/\./,$basefile))[0] . "_aln.$file_name_parts[$#file_name_parts]"; #i.e. file extension (added by PV April08)

   # run muscle as we want
   system(" $ENV{'EXE_MUSCLE'} < $infile > $tmpfile1 -stable -quiet ");
   system(" $ENV{'EXE_MUSCLE'} < $tmpfile1 > $tmpfile2 -stable -refine -quiet ");
   system(" $ENV{'EXE_MUSCLE'} < $tmpfile2 > $outfile -stable -refine -quiet ");

   # clean file
   if($clean==1){ system(" rm -f $tmpfile1 $tmpfile2 "); }

   # check the output file is there
   if(! -s $outfile){ return $phyTools::ERROR; }
   else{ return $outfile; }
}
#--------------------------------------------------------------------------------------#

# PV Jun 2014 f2p_v2
sub run_clustalo
{
    # outfmt='clu|fa|phy'
    my($infile,$outfmt) = @_;
    my @split = split(/\//,$infile); # mejor usar File::Basename
    my $basefile = $split[$#split];
    my @file_name_parts;
    my ($basename, $ext) = split(/\./,$basefile, 2);
    my $outfile = $basename . '_cluo' . ".$ext";
    
    # run clustalo
    system(" $ENV{'EXE_CLUSTALO'} -i $infile -o $outfile --outfmt=$outfmt")    
}
#--------------------------------------------------------------------------------------#

sub run_muscle3831    # PV, Sept 28th, 2011; also set $ENV{'EXE_MUSCLE3831'}   
{
   my($infile,$clean) = @_;

    my @split = split(/\//,$infile); # mejor usar File::Basename
    my $basefile = $split[$#split];
    my @file_name_parts;
    @file_name_parts = split(/\./,$basefile);
	
	
   my $tmpfile1 = (split(/\./,$basefile))[0] . "_tmp1.mus";
   my $tmpfile2 = (split(/\./,$basefile))[0] . "_tmp2.mus";
   my $outfile  = (split(/\./,$basefile))[0] . "_aln.$file_name_parts[$#file_name_parts]"; #i.e. file extension (added by PV April08)

   # run muscle with x2 refinement rounds after the primary alignment; note, -stable is not available anymore in muscle 3.8.31
   system(" $ENV{'EXE_MUSCLE3831'} < $infile > $tmpfile1 -quiet ");
   system(" $ENV{'EXE_MUSCLE3831'} < $tmpfile1 > $tmpfile2 -refine -quiet ");
   system(" $ENV{'EXE_MUSCLE3831'} < $tmpfile2 > $outfile -refine -quiet ");
   
   #  run_stable_py($file,$base.'_musAln.'.$ext);
   run_stable_py($infile,$outfile);
   
   # clean file
   if($clean==1){ system(" rm -f $tmpfile1 $tmpfile2 "); }

   # check the output file is there
   if(! -s $outfile){ return $phyTools::ERROR; }
   else{ return $outfile; }
}
#--------------------------------------------------------------------------------------#
     
sub run_stable_py     #  PV, Sept 28th, 2011; also set $ENV{'EXE_STABLEPY'}; to supplement the -stable option missing in muscle3.8.31
{
    my ($input_file,$aligned_file) = @_;
    my ($base,$ext) = (split /\./, $aligned_file)[0,1];
    my $stable_file = $base . '_stable.' . $ext;
    
    # auxiliary script stable_muscle.py, http://drive5.com/muscle/stable.html
    # NOTE: the script doesn't contain a shebang line; need to provide full-path to the executable
    system("python $ENV{'EXE_STABLEPY'} $input_file $aligned_file > $stable_file; cp $stable_file $aligned_file; rm $stable_file;"); 
}



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

sub run_mafft_einsi
{
    # PubMed=18372315
    # use to align intergenic regions
    my($infile) = @_;

    my @split = split(/\//,$infile); # mejor usar File::Basename
    my $basefile = $split[$#split];
    my @file_name_parts = split(/\./,$basefile);
	
    my $outfile  = (split(/\./,$basefile))[0] . "_aln.$file_name_parts[$#file_name_parts]"; #i.e. file extension (added by PV April08)

   # run mafft-einsi for INTERGENES, using Iterative refinement methods using WSP and consistency scores
   # mafft-einsi is applicable to a difficult problem such as RNA polymerase, 
   # which has several conserved motifs embedded in long unalignable regions. 
   # As E-INS-i has the minimum assumption of the three methods (i.e. einsi, linsi, ginsi) 
   # this is recommended if the nature of sequences to be aligned is not clear. 
   # Note that E-INS-i assumes that the arrangement of the conserved motifs is shared by all sequences.

   system(" $ENV{'EXE_MAFFT_EINSI'} --quiet $infile > $outfile ");
 
   # check the output file is there
   if(! -s $outfile){ return $phyTools::ERROR; }
   else{ return $outfile; }
}

sub run_NEIGHBOR
{
	# runs phylip neighbor from a distance matrix file and renames ouput according to $infile
	# as in/out files are always named the same cannot run more than one process simultaneosuly
	# in the same directory
	
	my($infile,$clean) = @_;

	my @split = split(/\//,$infile);
	my $basefile = $split[$#split];
	
   	my $tmp_infile  = dirname($infile) . "/infile";	
   	my $tmp_outfile = dirname($infile) . "/outfile"; 
   	my $tmp_outtree = dirname($infile) . "/outtree";
   	my $outfile     = (split(/\./,$basefile))[0] . "_nj.ph"; 
	my $outgramfile = (split(/\./,$basefile))[0] . "_nj_txt.graph"; 
   	my ($treeline,$treegram,$readOK) = ('','',0);
	
	# 1) copy input file to tmp input
	system("cp $infile $tmp_infile");

	# 2) run neighbor
	open(NEIGH,"cat $ENV{'NEIGHBOR_PARAMS_FILE'} | $ENV{'EXE_NEIGHBOR'} |") || die "# run_NEIGHBOR: cannot run $ENV{'EXE_NEIGHBOR'}\n";
	while(<NEIGH>)
	{
		if(/^# END marfil/){ last }
		elsif($readOK){ $treegram .= $_; }
		elsif(/^# BEGIN marfil/){ $readOK = 1 }
	}
	close(NEIGH);
	
	open(GRAM,">$outgramfile") || die "# run_NEIGHBOR: cannot create $outgramfile\n";
	print GRAM $treegram;
	close(GRAM);
	
			
	# 3) open the $tmp_outtree to read lines and get rid of \n and \s, and write the single treeline to file
     	open(OUTFILE, ">$tmp_outfile") || die "run_NEIGHBOR cannot write tmp_outfile $tmp_outfile: $!\n";
	open(TMP, $tmp_outtree) || die "run_NEIGHBOR cannot open outfile $outfile: $!\n";
    	while(<TMP>)
	{
     	  	$treeline .= $_;       # both operations, concat and s/\n//g, are required to get the single treeline!
     	  	$treeline =~ s/\n//g;
     	}
	close(TMP);
     	print OUTFILE $treeline;  #print "##########PHYTOOLS run_NEIGHBOR edited outfile is: $treeline\n\n";
     
     	# 4) rename output
	#system("cp $tmp_outtree $outfile"); 
        system("cp $tmp_outfile $outfile");
	
	# 5) clean file
   	if($clean==1){ system(" rm -f $tmp_infile $tmp_outfile $tmp_outtree"); }

  	# 6) check the output file is there
   	if(! -s $outfile){ return $phyTools::ERROR; }
   	else{ return ($outfile,$outgramfile); }
}


sub run_CDHIT
{
	# input: 1) .faa file 2) fraction of sequency identity cutoff ie 0.9
	my ($faafile,$outfile,$cutoff) = @_;

	my ($cdfile1,$cdfile2,$cdfile3) = ($faafile.'.ch',$faafile.'.ch.clstr',$faafile.'.ch.bak.clstr'); # ch for CdHit
	my $n_of_nr_seqs = 0;	

        open(CDHIT,"$ENV{'EXE_CDHIT'} -i $faafile -o $cdfile1 -c $cutoff |")
        || die "# find_redundancy_FASTA_CDHIT : cannot run $ENV{'EXE_CDHIT'} -i $faafile -o $cdfile1 -c $cutoff\n";
        while(<CDHIT>)
        {
        	#print;
		if(/(\d+) clusters/){ $n_of_nr_seqs = $1 }
        }
        close(CDHIT);

	if(!-s $cdfile1){ die "# run_CDHIT : failed, cannot find $cdfile1\n"; }
	else{ rename($cdfile1,$outfile) }

	unlink($cdfile2,$cdfile3);

	return $n_of_nr_seqs;
}

sub find_redundant_taxa_FASTA_CDHIT
{
	# input: 1) .faa file 2) fraction of sequency identity cutoff ie 0.9
	my ($faafile,$cutoff) = @_;

	my ($infile,$cdfile1,$cdfile2,$cdfile3) = ($faafile.'.short',$faafile.'.ch',$faafile.'.ch.clstr',$faafile.'.ch.bak.clstr'); # ch for CdHit
	my (%redundant_taxa,@redundant_strings,$header_string,$taxon,$seq);	

	# first make short version of FASTA headers to avoid cd-hit errors such as:
	# file reading error, possible a description line is too long
	my %fasta = read_FASTA_sequence($faafile);

	open(IN,">$infile") || die "# find_redundant_taxa_FASTA_CDHIT : cannot create $infile\n";
	foreach $seq (sort {$a<=>$b} keys(%fasta))
	{
		if($fasta{$seq}{'NAME'} =~ /(\[.*?\])/){ $header_string = ">$seq $1\n"; }
		else { $header_string = $fasta{$seq}{'NAME'} }
		print IN $header_string . $fasta{$seq}{'SEQ'} . "\n";
	}
	close(IN); 

	open(CDHIT,"$ENV{'EXE_CDHIT'} -i $infile -o $cdfile1 -c $cutoff |") 
	|| die "# find_redundancy_FASTA_CDHIT : cannot run $ENV{'EXE_CDHIT'} -i $infile -o $cdfile1 -c $cutoff\n";	
	while(<CDHIT>)
	{
		#print;
	}
	close(CDHIT);

	# check redundant sequences at $cdfile3
	if(!-s $cdfile3){ die "# find_redundancy_FASTA_CDHIT : cannot read $cdfile3, probably cd-hit failed!\n"; }
	open(OUT,$cdfile3);
	while(<OUT>)
	{
		#3       111aa, >GI:91802055 |[Nit... at 70%i
		#print;
		if(/\d+aa, >(.*?) at \d+%/)
		{
			$header_string = $1; 
			$header_string =~ s/\.\.\.//;
			$header_string = (split(/[\[|\(]/,$header_string))[0]; # to avoid cases such as GI:74420406 |[Nit ...
			push(@redundant_strings,$header_string); #print "# $header_string\n";	
		}
	}
	close(OUT);

	if(@redundant_strings)
	{	
		%fasta = read_FASTA_sequence($infile);
		foreach $seq (keys(%fasta))
		{
			$taxon = '';
			foreach $header_string (@redundant_strings)
			{
				if($fasta{$seq}{'NAME'} =~ $header_string)
				{
					if($fasta{$seq}{'NAME'} =~ /(\[.*?\])/)
					{ 
						$taxon = $1; 
						$redundant_taxa{$taxon} = $seq; #print "# $taxon $seq\n";
						last;
					}
				}
			}
		}		
	}

	unlink($infile,$cdfile1,$cdfile2,$cdfile3);

	return %redundant_taxa;
}

sub run_BRANCHCLUST
{
	# runs branchclust in order to cluster tree leafs
	# if $filter_paralogs == 1 reads families, otherwise it reads clusters
	 
	my($treefile, $ref_taxa, $filter_paralogs, $cluster_size, $clean) = @_;
	
	my (%clusters,$taxa,$member,$n_of_taxa,$n_of_clusters,$readOK,$pattern,$tree_string);
	
   	my $taxafile     = dirname($treefile) . "/gi_numbers.out";	
   	my $outclusters  = dirname($treefile) . "/clusters.out";
   	my $outlog       = dirname($treefile) . "/cluster.log";
	my $outfamilies  = dirname($treefile) . "/families.list";
	
	# 1) write taxa file
	open(TAXA,">$taxafile") || die "# run_BRANCHCLUST: cannot write to $taxafile $!\n";
	foreach $taxa (keys (%{$ref_taxa}))
	{
		print TAXA "$taxa $ref_taxa->{$taxa}{'SIZE'} | ";
		foreach $member (@{$ref_taxa->{$taxa}{'MEMBERS'}}){ printf TAXA "$member\t"; }
		print TAXA "\n";
	}
	close(TAXA);
	
	# 2) run branchclust
	system("$ENV{'EXE_BRANCHCLUST'} $treefile $cluster_size > /dev/null");
	
	# 3) parse clusters file
	if($filter_paralogs) { $pattern = "------------ FAMILY ------------"; }
	else { $pattern = "------------ CLUSTER -----------"; }

	$n_of_clusters = $readOK = 0;
	open(CLUSTERS,$outclusters) || die "# run_BRANCHCLUST: cannot read $outclusters\n";
	while(<CLUSTERS>)
	{
		if($readOK)

		{
			my @cluster = split;
			$clusters{$n_of_clusters} = [ @cluster ] ;
			$readOK = 0
		}
		elsif(/^$pattern/) 
		{
			$n_of_clusters++;
			$readOK = 1;
		}
	}
	close(CLUSTERS);

	# 4) label tree leafs with cluster numbers
	$tree_string = '';
	open(TREE,$treefile) || die "# run_BRANCHCLUST: cannot read $treefile\n";
	while(<TREE>){ $tree_string .= $_; }
	close(TREE);	

	foreach $n_of_clusters (keys(%clusters))
	{
		foreach $member (@{$clusters{$n_of_clusters}})
		{
			$tree_string =~ s/$member/$member\___$n_of_clusters/; #print "#$member#\n";
		}
	}

	open(TREE,">$treefile") || die "# run_BRANCHclust: cannot rewrite $treefile\n";
	print TREE $tree_string;
	close(TREE);

	# 5) clean if required
	if($clean==1){ system(" rm -f $taxafile $outclusters $outlog $outfamilies"); }

	return %clusters;
}

sub run_mablock
{
   my($infile,$clean) = @_;

	my @split = split(/\//,$infile);
	my $basefile = $split[$#split];
	
   my $filename = (split(/\./,$basefile))[0];	
   my $outfile  = (split(/\./,$basefile))[0] . ".blks";
   my $logfile  = (split(/\./,$basefile))[0] . "_log"; 
   my $seqsfile = (split(/\./,$basefile))[0] . ".seqs"; 
	
   #print "my infile in run_mablock is $infile\n";
	#print "my outfile in run_mablock is $outfile\n";
   #print "my logfile in run_mablock is $logfile\n"; 
   
   # run mablocks as we want it
   system("$ENV{'EXE_MABLOCK'} $infile $filename B 10 50 > $logfile");
   
   # clean file
   if($clean==1){ system(" rm -f $logfile $seqsfile"); }

   # check the output file is there
   if(! -s $outfile){ return $phyTools::ERROR; }
   else { return $outfile; }
}

sub extract_codon_freqs_from_fasta2spsumfile
{
	# assumes introns have been removed Nov2008
	my ($ref_hash_FASTA,$label,$spsum_filename,$taxon) = @_; 
 
 	my @codons = qw/CGA CGC CGG CGT AGA AGG CTA CTC CTG CTT TTA TTG TCA TCC TCG TCT AGC AGT ACA ACC ACG ACT CCA CCC CCG CCT GCA GCC GCG GCT GGA GGC GGG GGT GTA GTC GTG GTT AAA AAG AAC AAT CAA CAG CAC CAT GAA GAG GAC GAT TAC TAT TGC TGT TTC TTT ATA ATC ATT ATG TGG TAA TAG TGA/;

	my ($codonFreq,$no_of_seqs,$no_of_codons,%codonTable,@triplets,$sequence,$triplet) = ('',0,0);

	foreach my $seq (keys(%{$ref_hash_FASTA}))
        {
		next if($taxon && $ref_hash_FASTA->{$seq}{'NAME'} !~ /$taxon/);

		$no_of_seqs++;
		$sequence = $ref_hash_FASTA->{$seq}{'SEQ'};
		chomp($sequence);
		# read the sequence until the end progressing three positions each time and count each triplet that has been found
		for(my $i=0; $i < length($sequence); $i +=3)
		{
			$triplet = uc(substr($sequence,$i,3));
			push(@triplets, $triplet);
			$codonTable{$triplet}++ ;
			$no_of_codons = scalar(@triplets);
		}
	}
        
	foreach my $codon (@codons)
	{
		if(grep(/$codon/,@triplets))			
		{
			$codonFreq .= $codonTable{$codon} . " ";
		}
		else{ $codonFreq .= "0 "; }
	}
	
	# print outfile 
	open(SPSUM,">$spsum_filename") || die "# extract_codon_freqs_from_fasta2spsumfile : cannot create $spsum_filename\n";
	print SPSUM "$label: $no_of_codons\n$codonFreq\n";
	close SPSUM;
}

sub run_codehop
{
	# ojo, codehop tiene que correr en el directorio examples a huevo
	# user table puede ser más de una tabla separada por comas: tabla1,tabla2,...

   	my($infile, $ref_AAFasta, $user_table, $user_codehopT, $user_CUT_path ) = @_; #CUT=codon usage table
  
	if(!$user_CUT_path){ $user_CUT_path = $ENV{'BLIMPS_DOCS'} }
 
   	my (@codon_files,$table,$short_table,$tablename,$outfile,@outfiles); 
   
	my @split = split(/\//,$infile);
   	my $basefile = $split[$#split];
   	my $default_codon_usage_file = "$ENV{'BLIMPS_DOCS'}/atumC58.codon.use"; #does not work any more
	my $codehopT = $ENV{'DEF_CODEHOP_T'};
	if($user_codehopT){ $codehopT = $user_codehopT }
   
   
   	# 1) infer preferred codon usage tables for each sequence in $infile
	if($user_table)
	{
		$user_table =~ s/ //g; 
		#print "# run_codehop: user table(s) used = $user_table\n"; 
		# not needed anymore, since a ranking of tables will be printed out
		push(@codon_files,split(/,/,$user_table));	
	}
	else
	{ 
		my %codontables = select_codon_tables($ref_AAFasta);
		if(!keys(%codontables))
		{
			print "# run_codehop: no codon usage table was selected for $infile,".
				" default = $default_codon_usage_file\n";
			push(@codon_files,$default_codon_usage_file);
		}
		else
		{
			foreach $table (keys(%codontables))
			{
				# 1.1) run my_coduse with this table and get table file name
				#print "#$table#$codontables{$table}#\n";	
				$tablename = $table;
				$tablename =~ s/\s+|\/|\(|\)/\_/g; # remove chars that cause trouble in the shell
				$tablename =~ s/\.//g;
				$tablename .= ".codon.use"; 
			
				if(-s "$ENV{'GB_SPSUM_DIR'}/$codontables{$table}")
				{
					#if(!-s "$ENV{'BLIMPS_DIR'}/docs/$tablename") # uncomment only if coduse takes too long
					#{
					system("$ENV{'BLIMPS_DIR'}/bin/linux/coduse $ENV{'GB_SPSUM_DIR'}/$codontables{$table} \"$table\" $tablename > /dev/null"); 
					system("mv -f $tablename $ENV{'BLIMPS_DOCS'}");
					#}

					push(@codon_files,$tablename);
				}
			}
		}
   	}

	@codon_files = rank_codon_usage_tables($user_CUT_path,@codon_files);
		
	# 2) remember where we were before running codehop and move to codehop directory
	my $cwdDIR = getcwd();
	chdir "$ENV{'BLIMPS_DIR'}/examples/codehop";
 
 	# 3) run codehop with each selected table
	#print "#$user_CUT_path#\n";
 	foreach $tablename (@codon_files)
	{
 		$short_table = basename($tablename);
		$outfile  = (split(/\./,$basefile))[0] . "__" . $short_table . "__codehops.out";

		if($user_CUT_path ne $ENV{'BLIMPS_DOCS'}) # por alguna razon falla con este directorio y solo con este!!
		{		
			system("./my_codehop_path.csh $cwdDIR/$infile $cwdDIR/$outfile $tablename $codehopT $user_CUT_path");
		}
		else{ system("./my_codehop.csh $cwdDIR/$infile $cwdDIR/$outfile $tablename $codehopT"); }
   		
		if(-s $cwdDIR."/".$outfile)
		{ 
			#print `cat $cwdDIR/$outfile`;
			push(@outfiles,$outfile); 
		}
	}
	
	# 4) return to working directory
	chdir $cwdDIR;

   	return @outfiles; 
}

# sub notas_de_CUD
#{wget ftp://ftp.kazusa.or.jp/pub/codon/current/gbbct.spsum
#../../bin/linux/coduse ../../../../codehops/codehops/codons/GB104/gbbct.spsum "Escherichia colx" kk
#$ENV{'BLIMPS_DIR'} =  "$ENV{'MARFIL'}bin/BLOCKS3.7_software/blimps3.7/blimps-3.7/";
#$blimps_bin = "$ENV{'BLIMPS_DIR'}/bin/linux/";
#$docs = "$ENV{'BLIMPS_DIR'}/docs/";
#system("$bin/codehop $ARGV[0] -C$docs/arab.codon.use -O3 >& $ARGV[1]");
#system("$bin/codehop $ARGV[0] -C$docs/arab.codon.use -O3");
#BLIMPS_DIR=$ENV{'MARFIL'}bin/BLOCKS3.7_software/blimps3.7/blimps-3.7/
#necesitamos un sitio fijo para gbbct.spsum : $ENV{'MARFIL'}bin/
#}

sub rank_codon_usage_tables
{
	# ranks an array of codon usage tables from less to more redundant 
	# distance between tables is computed as cos(angle) between their codon frequencies vectors

	my ($CUT_path,@codon_usage_tables) = @_;

	my (@nr_tables,$table,%tables,%n_of_codons,@table_names);
	my ($t1,$t2,$n_of_tables,$codon,$freq1,$freq2,%cos_angles);
	my ($escalar,$mod1,$mod2,$cos,$short_table);

	# 0) return if only one table
	if(scalar(@codon_usage_tables) == 1)
	{
		push(@nr_tables,$codon_usage_tables[0]);
		return @nr_tables;
	}

	
	# 1) read original tables into arrays
	foreach $table (@codon_usage_tables)
	{
		#sleep(1); system("ls $ENV{'BLIMPS_DOCS'}/$table"); # refresh pinche NFS
		open(TABLE,$CUT_path."/".$table) || die "# rank_codon_usage_tables : cannot read $CUT_path/$table\n";
		while(<TABLE>)
		{
			if(/^\d\./){ push(@{$tables{$table}},(split)[0]); }
			elsif(/: (\d+) codons/){ $n_of_codons{$table} = $1; }
		}
		close(TABLE);
	}
	
	# 2) compare all against all, adding up their pairwise distances
	$n_of_tables = scalar(keys(%tables));
	@table_names = (sort (keys(%n_of_codons)));
	
	for($t1=0;$t1<$n_of_tables;$t1++)
	{
		for($t2=$t1+1;$t2<$n_of_tables;$t2++)
		{
			($escalar,$mod1,$mod2,$cos) = (0,0,0,0);
			for($codon=0;$codon<64;$codon++)
			{
				$freq1 = $tables{$table_names[$t1]}[$codon];
				$freq2 = $tables{$table_names[$t2]}[$codon];
				
				$escalar += $freq1*$freq2;
				if($freq1 > 0.000001){ $mod1 += $freq1**2; }
				if($freq2 > 0.000001){ $mod2 += $freq2**2; }
			}
		
			if($mod1){ $mod1 = sqrt($mod1) }
			if($mod2){ $mod2 = sqrt($mod2) }
	
			if($mod1 && $mod2)
			{ 
				$cos = $escalar/($mod1*$mod2); 
				#printf("# cos(angulo) = %1.4f %s <=> %s\n",$cos,$table_names[$t1],$table_names[$t2]);
	
				$cos_angles{$table_names[$t1]} += $cos;
				$cos_angles{$table_names[$t2]} += $cos;
			}
		}
		
		$cos_angles{$table_names[$t1]} = sprintf("%1.2f",($cos_angles{$table_names[$t1]}/($n_of_tables-1)) );
	}
	
	# 3) rank tables in terms of cos_angles to all other tables, ie, redundant tables should go at the end 	
	foreach $table (sort {$cos_angles{$a}<=>$cos_angles{$b}} (keys(%cos_angles)))
	{
		print "# rank_codon_usage_tables : $table $cos_angles{$table}\n";
		push(@nr_tables,$table);
	}	

	return @nr_tables;	
}	

sub select_intergene_amplicons
{
	# function called from flankingCDS2primers.pl that receives 'left' and 'right' block files, 
	# codehops and strands of two CDS segments that flank an intergenic region
	# OJO: falta poner codigo para intrones, Jul2010
	my ($refH_blockfiles,$refH_primerfiles,$refH_blockstrand) = @_;

	my ($read,%block_coordinates,%block_lengths,%n_of_blocks,%primers,$complement,$offsetnt);
        my ($primerN,$primerC,$offset,$offset_codon,$AA_aln_coord,$AA_aln_coordC,$blockN,$blockC);
        my ($ampl_score,$ampl_length,$diff_temperature,$diff_degen,$diff_score,$diff_length,%aln);
        my (%amplicons,%selected_amplicons,$ampl_length2,$cover,@ampl_coords,$amplicon,$midpoint);
        my ($N,$C,%introns,%aln_introns,$length,$Clength,$intron,$intronFREE,$n_of_intron_primers);

	foreach my $side (keys(%$refH_blockfiles))
	{
		my $BLOCKS_file = $refH_blockfiles->{$side}; #print "#> $side\n";
		my $CODEHOPs_file = $refH_primerfiles->{$side}; #print "# $BLOCKS_file $CODEHOPs_file\n";

		# 1) read BLOCKS file and store coordinates for each block
	        $read=0;
	        open(BLOCKS,$BLOCKS_file) || die "# select_intergene_amplicons : cannot read $BLOCKS_file\n";
	        while(<BLOCKS>)
	        {
			#BL   UNK motif;  width=49; seqs=5; 99.5%=0; strength=0
			#001                (  11) VLTGIEPISIIDEMQRSYLDYAMSVIVSRALPDVRDGLKPVHRRILHAM 9.438889
	                if(/^BL   UNK motif;  width=(\d+)/){ push(@{$block_lengths{$side}},$1); $read=1; }
	                elsif($read && /\w+\s+\(\s*(\d+)/){ push(@{$block_coordinates{$side}},$1); $read=0; }
	        }
	        close(BLOCKS);

		 # 2) parse CODEHOP file and store primers for each block
	        $n_of_blocks{$side}=-1;
	        $complement=$read=$n_of_intron_primers=0;
	        open(CODEHOP,$CODEHOPs_file) || die "# select_intergene_amplicons : cannot read $CODEHOPs_file\n";
	        while(<CODEHOP>)
	        {
			#Processing Complement of Block 39x5215xgB
			#G  F  G  D  S  S  S  N  Q  R  N  T  F  G  Q  N
                        #      ttrttrwbngtTAGACCACTTAGAAAAAGAG -5'  Core: degen=96 len=11  Clamp: score=69, len=20 temp= 55.6
	                if(/CLAMP NEEDS EXTENSION/){ next; }
	                elsif(/^Processing Block/){ $n_of_blocks{$side}++; $complement=0; }
	                elsif(/^Processing Complement/){ $complement=1; }
	                elsif(/(\s*)(\w+) -[3|5]'\s+Core: degen=(\d+) len=(\d+)\s+Clamp: score=(\d+), len=(\d+) temp= (\d+.\d)/)
	                {
                        	($Clength,$length) = ($4,$6);
				$offsetnt = length($1);
        	                $offset = int($offsetnt/3);
	                        $offset_codon = $offsetnt%3;
	                        #$AA_aln_coord = $block_coordinates{$side}[$n_of_blocks{$side}] + $offset - 1; # -1 adds one extra aa
	                        #$AA_aln_coordC = $AA_aln_coord+int(($Clength+$length)/3);

                	        $primers{$side}[$n_of_blocks{$side}]{$2}{'complement'} = $complement;
	                        $primers{$side}[$n_of_blocks{$side}]{$2}{'Cdegen'} = $3;
	                        $primers{$side}[$n_of_blocks{$side}]{$2}{'Clength'} = $4;
	                        $primers{$side}[$n_of_blocks{$side}]{$2}{'score'} = $5;
	                        $primers{$side}[$n_of_blocks{$side}]{$2}{'length'} = $6;
	                        $primers{$side}[$n_of_blocks{$side}]{$2}{'temperature'} = $7;
	                        $primers{$side}[$n_of_blocks{$side}]{$2}{'offset'} = $offset_codon;
				$primers{$side}[$n_of_blocks{$side}]{$2}{'offset_primer'} = $offsetnt + 
					($block_coordinates{$side}[$n_of_blocks{$side}]*3);
	                        #$primers{$side}[$n_of_blocks{$side}]{$2}{'block_length'} = $block_length{$side}[$n_of_blocks{$side}];
	                        #print "$1 $2 $3 $4 $5 $6 $7 $offset $offset_codon $complement\n";
                	}
		}
        	close(CODEHOP); 
	} #print "# $n_of_blocks{'left'} $n_of_blocks{'right'}\n";

	# 3) combine primers from both sides from N to Cterm and exclude bad primer pairs
        for($blockN=0;$blockN<=$n_of_blocks{'left'};$blockN++)
        {
                for($blockC=0;$blockC<=$n_of_blocks{'right'};$blockC++)
                {
                        #print "$blockN $blockC\n";
                        foreach $primerN (keys(%{$primers{'left'}[$blockN]}))
                        {
				# 5->3 complement=0
				# 3<-5 complement=1
				# leftCDS(1) 5->3  ----------- 3<-5 rightCDS(1)
				# leftCDS(-1) 3<-5 ----------- 3<-5 rightCDS(1)
				# leftCDS(1) 5->3  ----------- 5->3 rightCDS(-1)
				# leftCDS(-1) 3<-5 ----------- 5->3 rightCDS(-1)

				next if($primers{'left'}[$blockN]{$primerN}{'complement'}==0 && 
                                        $refH_blockstrand->{'left'} eq '-1');
				next if($primers{'left'}[$blockN]{$primerN}{'complement'}==1 &&
                                        $refH_blockstrand->{'left'} eq '1');
				
                                foreach $primerC (keys(%{$primers{'right'}[$blockC]}))
                                {
                                        next if ($primers{'right'}[$blockC]{$primerC}{'complement'}==0 &&
						$refH_blockstrand->{'right'} eq '1');
					next if ($primers{'right'}[$blockC]{$primerC}{'complement'}==1 &&
                                                $refH_blockstrand->{'right'} eq '-1');

                                        # check other amplicon properties
                                        $ampl_score = 0;

                                        $diff_temperature = abs($primers{'left'}[$blockN]{$primerN}{'temperature'}-
								$primers{'right'}[$blockC]{$primerC}{'temperature'});
                                        next if($diff_temperature > $ENV{"MAX_PRIMER_TEMP_DIFF"});

					$diff_degen = abs((log($primers{'left'}[$blockN]{$primerN}{'Cdegen'})/log(2))-
							(log($primers{'right'}[$blockC]{$primerC}{'Cdegen'})/log(2)));
                                        next if($diff_degen > $ENV{"MAX_PRIMER_DEGEN_DIFF"});

                                        $diff_score = abs($primers{'left'}[$blockN]{$primerN}{'score'}-
							$primers{'right'}[$blockC]{$primerC}{'score'});
                                        next if($diff_score > $ENV{"MAX_PRIMER_SCORE_DIFF"});

                                        $diff_length = abs( ($primers{'left'}[$blockN]{$primerN}{'Clength'}+
								$primers{'left'}[$blockN]{$primerN}{'length'}) - 
								($primers{'right'}[$blockC]{$primerC}{'Clength'}+
								$primers{'right'}[$blockC]{$primerC}{'length'}) ); #print "#$diff_length#\n";
                                        next if($diff_length > $ENV{"MAX_PRIMER_LENGTH_DIFF"});

                                        # attempt to calculate an amplicon score
                                        $ampl_score = $primers{'left'}[$blockN]{$primerN}{'score'} + 
							$primers{'right'}[$blockC]{$primerC}{'score'};

                                        #print "|$ampl_score|$diff_temperature|$diff_degen|$diff_score|$primerN <-> $primerC\n";

                                        # store amplicons with their scores
                                        $amplicons{$ampl_score}{$primerN}{$primerC} = [ $primers{'left'}[$blockN]{$primerN} , 
										$primers{'right'}[$blockC]{$primerC} ];
                                }
                        }
                }
        }					
	return %amplicons;
}

sub select_amplicons
{
	my ($FASTA_aln_file,$BLOCKS_file,$CODEHOPs_file,$ref_seen_amplicons,$length_range) = @_;
	
	my ($read,@block_coordinates,@block_lengths,$n_of_blocks,@primers,$complement);
	my ($primerN,$primerC,$offset,$offset_codon,$AA_aln_coord,$AA_aln_coordC,$blockN,$blockC);
	my ($ampl_score,$ampl_length,$diff_temperature,$diff_degen,$diff_score,$diff_length,%aln);
	my (%amplicons,%selected_amplicons,$ampl_length2,$cover,@ampl_coords,$amplicon,$midpoint);
	my ($N,$C,%introns,%aln_introns,$length,$Clength,$intron,$intronFREE,$n_of_intron_primers);
	

	# -1) check amplicon length parameters 
	my ($min_length,$max_length) = ($ENV{"MIN_AMPLICON_LENGTH"},$ENV{"MAX_AMPLICON_LENGTH"});
	if($length_range)
	{
		($min_length,$max_length) = split(/:/,$length_range);
	}

	# 0) read FASTA sequence
	%aln = read_FASTA_sequence($FASTA_aln_file);
 	%introns = find_introns_FASTA_headers(\%aln,1);
	if(%introns){ %aln_introns = adjust_introns_alignment(\%aln,\%introns); }

	# 1) read BLOCKS file and store coordinates for each block
	$read=0;
	open(BLOCKS,$BLOCKS_file) || die "# select_amplicons : cannot read $BLOCKS_file\n";
	while(<BLOCKS>)
	{
		if(/^BL   UNK motif;  width=(\d+)/){ push(@block_lengths,$1); $read=1; }
		elsif($read && /\w+\s+\(\s*(\d+)/){ push(@block_coordinates,$1); $read=0; }
	}
	close(BLOCKS);
	
	# 2) parse CODEHOP file and store primers for each block
	$n_of_blocks=-1;
	$complement=$read=$n_of_intron_primers=0;
	open(CODEHOP,$CODEHOPs_file) || die "# select_amplicons : cannot read $CODEHOPs_file\n";
	while(<CODEHOP>)
	{
		if(/CLAMP NEEDS EXTENSION/){ next; }
		elsif(/^Processing Block/){ $n_of_blocks++; $complement=0; }
		elsif(/^Processing Complement/){ $complement=1; }
		elsif(/(\s*)(\w+) -[3|5]'\s+Core: degen=(\d+) len=(\d+)\s+Clamp: score=(\d+), len=(\d+) temp= (\d+.\d)/) 
		{
			($Clength,$length) = ($4,$6);
			$offset = int(length($1)/3);
                        $offset_codon = length($1)%3;
                        $AA_aln_coord = $block_coordinates[$n_of_blocks] + $offset - 1; # -1 adds one extra aa
			$AA_aln_coordC = $AA_aln_coord+int(($Clength+$length)/3);

			$intronFREE=-1;
			foreach $intron (keys(%aln_introns)){ if($intron >= $AA_aln_coord && $intron <= $AA_aln_coordC){ $intronFREE=$intron; last }}
			if($intronFREE>-1)
			{
				#print "# select_amplicons : skip primer $2 , contains introns ($AA_aln_coord , $intronFREE , $AA_aln_coordC)\n";	
				$n_of_intron_primers++;
				next;
			}

			$primers[$n_of_blocks]{$2}{'complement'} = $complement;
			$primers[$n_of_blocks]{$2}{'Cdegen'} = $3;
			$primers[$n_of_blocks]{$2}{'Clength'} = $4;
			$primers[$n_of_blocks]{$2}{'score'} = $5;
			$primers[$n_of_blocks]{$2}{'length'} = $6;
			$primers[$n_of_blocks]{$2}{'temperature'} = $7;
			$primers[$n_of_blocks]{$2}{'offset'} = $offset_codon;	
			$primers[$n_of_blocks]{$2}{'AA_align_coord'} = $AA_aln_coord;	
			#print "$1 $2 $3 $4 $5 $6 $7 $offset $offset_codon $AA_aln_coord $complement\n";
		}	
	}
	close(CODEHOP);
	print "# number of intron-spanning primers skipped = $n_of_intron_primers\n" if($n_of_intron_primers);

	# 3) recombine primers from left(N) to right(C) (blocks are sorted that way) 
	for($blockN=0;$blockN<$n_of_blocks;$blockN++)
	{
		for($blockC=$blockN+1;$blockC<=$n_of_blocks;$blockC++)
		{
			#print "$blockN $blockC\n";
			foreach $primerN (keys(%{$primers[$blockN]}))
			{
				next if ($primers[$blockN]{$primerN}{'complement'}==1);

				foreach $primerC (keys(%{$primers[$blockC]}))
				{
					next if ($primers[$blockC]{$primerC}{'complement'}==0);
	
					# check amplicon length			
					$ampl_length = 3 * calc_amplicon_length_aa( $primers[$blockN]{$primerN}{'AA_align_coord'},$primers[$blockC]{$primerC}{'AA_align_coord'}, %aln);
					# correct amplicon length if introns are present!!
					foreach $intron (keys(%aln_introns))
					{ 
						if($intron >= $primers[$blockN]{$primerN}{'AA_align_coord'} && 
							$intron <= $primers[$blockC]{$primerC}{'AA_align_coord'})
						{
							$ampl_length += $aln_introns{$intron};
						}		
					}

					next if($ampl_length < $min_length || $ampl_length > $max_length);
					#print "$ampl_length $min_length $max_length\n";
				
					# check other amplicon properties
					$ampl_score = 0;
					
					$diff_temperature = abs($primers[$blockN]{$primerN}{'temperature'}-$primers[$blockC]{$primerC}{'temperature'});
					next if($diff_temperature > $ENV{"MAX_PRIMER_TEMP_DIFF"});
					
					$diff_degen = abs((log($primers[$blockN]{$primerN}{'Cdegen'})/log(2))-(log($primers[$blockC]{$primerC}{'Cdegen'})/log(2)));
					next if($diff_degen > $ENV{"MAX_PRIMER_DEGEN_DIFF"});
					
					$diff_score = abs($primers[$blockN]{$primerN}{'score'}-$primers[$blockC]{$primerC}{'score'});
					next if($diff_score > $ENV{"MAX_PRIMER_SCORE_DIFF"});
					
					$diff_length = abs( ($primers[$blockN]{$primerN}{'Clength'}+$primers[$blockN]{$primerN}{'length'}) - ($primers[$blockC]{$primerC}{'Clength'}+$primers[$blockC]{$primerC}{'length'}) );#print "#$diff_length#\n";
					next if($diff_length > $ENV{"MAX_PRIMER_LENGTH_DIFF"});


					# attempt to calculate an amplicon score
					$ampl_score = $primers[$blockN]{$primerN}{'score'} + $primers[$blockC]{$primerC}{'score'};
					
					#print "|$ampl_score|$diff_temperature|$diff_degen|$diff_score|$primerN <-> $primerC\n";
					
					# store amplicons with their scores
					$amplicons{$ampl_score}{$primerN}{$primerC} = [ $primers[$blockN]{$primerN} , $primers[$blockC]{$primerC} ];
				}
			}
		}
	}
	
	# 4) select non-overlapping amplicons to ensure maximal gene coverage
	foreach $ampl_score (sort {$b<=>$a} (keys(%amplicons)))
	{
		foreach $primerN (keys(%{$amplicons{$ampl_score}}))
		{
			foreach $primerC (keys(%{$amplicons{$ampl_score}{$primerN}}))
			{
				#print "$ampl_score $primerN $primerC\n";
				
				$N = $amplicons{$ampl_score}{$primerN}{$primerC}[0]{'AA_align_coord'};
				$C = $amplicons{$ampl_score}{$primerN}{$primerC}[1]{'AA_align_coord'} + int(length($primerC)/3);	
				
				# re-check amplicon length as unseen primer pairs can get through (with same score and primer N, diff primerC)
				$ampl_length = 3 * ($C-$N+1);
				next if($ampl_length < $min_length || $ampl_length > $max_length);
			
				if(!@ampl_coords)
				{ 
					next if(amplicon_too_similar($ref_seen_amplicons,$N,$C));
					
					push(@ampl_coords,[$N,$C]); 
					$selected_amplicons{$N}{$ampl_score}{$primerN}{$primerC} = $amplicons{$ampl_score}{$primerN}{$primerC};
					$ref_seen_amplicons->{$N}{$C} = 1;
					print "# amplicon coords: $N $C\n";
				}
				else
				{
					$ampl_length2 = $C - $N + 1;

					foreach $amplicon (@ampl_coords)
					{
						$ampl_length = $amplicon->[1] - $amplicon->[0] + 1;
						
						if($ampl_length < $ampl_length2)
						{
							$midpoint = $amplicon->[0] + ($ampl_length/2);
							if($midpoint < $N || $midpoint > $C){ $cover = 0; }
							elsif($midpoint < $C)
							{
								$cover = ($C - $amplicon->[0]+1)/$ampl_length;
							}
							else
							{
								$cover = ($amplicon->[1] - $N+1)/$ampl_length2;
							}#print "< $cover\n";
						}
						else
						{
							$midpoint = $N + ($ampl_length2/2);
							if($midpoint < $amplicon->[0] || $midpoint > $amplicon->[1]){ $cover = 0; }
							elsif($midpoint < $amplicon->[1])
							{
								$cover = ($amplicon->[1] - $N+1)/$ampl_length2;
							}
							else
							{
								$cover = ($C - $amplicon->[0]+1)/$ampl_length;
							}#print "> $cover\n";
						}
						last if($cover >= $ENV{"MAX_AMPLICON_COVER"});#stop checks here print "# $cover\n";
					}
				
					next if(amplicon_too_similar($ref_seen_amplicons,$N,$C));#print "= $cover\n";				

					if($cover < $ENV{"MAX_AMPLICON_COVER"})
					{
						push(@ampl_coords,[$N,$C]); 
						$selected_amplicons{$N}{$ampl_score}{$primerN}{$primerC} = $amplicons{$ampl_score}{$primerN}{$primerC};
						$ref_seen_amplicons->{$N}{$C} = 1;
						print "# amplicon coords: $N $C\n";
					}
				}
			}
		}	
	} 
	
	return %selected_amplicons;
}

sub amplicon_too_similar
{
	# check whether a new amplicon in the form of N,C coordinates is too similar
	# to amplicons already stored in $ref_seen_amplicons
	# returns 1 in case amplicon N,C found to be too similar, 0 otherwise
	# updated Jul2010 to check degeneracy levels, stored in $ref_seen_amplicons, if required
	
	my ($ref_seen_amplicons,$N,$C,$is_intergenic,$degeneracy) = @_;
	
	my $lim = $ENV{'MAX_AMP_MARGIN_DEV'};	
	if($is_intergenic){ $lim = $ENV{'MAX_INTAMP_MARGIN_DEV'}; }
	my ($pos,$pos2);
	
	if(!$degeneracy)
	{
		if($ref_seen_amplicons->{$N}{$C}){ return 1; }
	
		for($pos=$N-$lim;$pos<$N+$lim;$pos++)
		{
			if($ref_seen_amplicons->{$pos}{$C}){ return 1 }
		}
		
		for($pos=$C-$lim;$pos<$C+$lim;$pos++)
		{
			if($ref_seen_amplicons->{$N}{$pos}){ return 1 }
		}
		
		for($pos=$N-$lim;$pos<$N+$lim;$pos++)
		{
			for($pos2=$C-$lim;$pos2<$C+$lim;$pos2++)
			{	
				if($ref_seen_amplicons->{$pos}{$pos2}){ return 1 }
			}
		}
	}
	else
	{
		if($ref_seen_amplicons->{$N}{$C} && 
			$degeneracy >= $ref_seen_amplicons->{$N}{$C}){ return 1; }
         
                for($pos=$N-$lim;$pos<$N+$lim;$pos++)
                {
                        if($ref_seen_amplicons->{$pos}{$C} &&
				$degeneracy >= $ref_seen_amplicons->{$pos}{$C}){ return 1 }
                }
         
                for($pos=$C-$lim;$pos<$C+$lim;$pos++)
                {
                        if($ref_seen_amplicons->{$N}{$pos} &&
				$degeneracy >= $ref_seen_amplicons->{$N}{$pos}){ return 1 }
                }
         
                for($pos=$N-$lim;$pos<$N+$lim;$pos++)
                {
                        for($pos2=$C-$lim;$pos2<$C+$lim;$pos2++)
                        {
                                if($ref_seen_amplicons->{$pos}{$pos2} &&
					$degeneracy >= $ref_seen_amplicons->{$pos}{$pos2}){ return 1 }
                        }
                }
	}
					
	return 0;
}

sub calc_amplicon_length_aa
{
	my ($Ncoord,$Ccoord,%aln) = @_;
	
	my ($seq,$amplicon,$length,$minlength);
	
	$minlength = $ENV{'MAX_AMPLICON_LENGTH'};
	foreach $seq (keys(%aln))
	{
		$amplicon = substr($aln{$seq}{'SEQ'},$Ncoord,$Ccoord-$Ncoord+1);
		#print $amplicon,"\n";
		$amplicon =~ s/-//g;
		$length = length($amplicon);
		if($length < $minlength){ $minlength = $length; }	
	}
	
	return $minlength;
}


sub cut_intergen_amplicon_aln_files
{
	my ($refH_amplicons,$refH_sequences,$refH_strands,$refH_DNA_aln,$refH_lengthflank,$rootfilename) = @_;
	
	my (@amplicon_file_names,@primers_file_names,$pos,$N,$C,$score,%intergene);
        my (%DNA_aln,$seq,$amplicon,$header,$N_prop_matches,$C_prop_matches,$primerN);
        my ($n_of_amplicons,$n_of_matches,$primerpos,$pair,$left,$right,$revright) = (1);
        my ($DNA_aln_file_name,$dna_filename,$amplicon_match_aln,$primer_filename);	
	my ($N_CDScoord,$C_CDScoord,%alignedCDS,%unalignedCDS,$length,$tmpfile1,$side);
	my ($C_strand,$N_strand,$N_primer_sense,$C_primer_sense,$tmpfile2) = ('1','1');
	my (%FASTA,$seqnozeros,$sequence,$take_complement,$primersC);
	my ($left_intergene_aligned_length,$total_aligned_length,$right_primer_length);	
	my @split = split(/\//,$rootfilename); 
	my $basefile = $split[$#split];
	my @file_name_parts = split(/\./,$basefile);	
	my $root = $file_name_parts[0] . '_aln'; 

	# first mutiply align intergenic sequence
        $tmpfile1 = '_tmp_'.$root.'_'.scalar(keys(%$refH_amplicons)).'_'.'.fna';
	open(INTFNA,">$tmpfile1") || 
		die "# cut_intergen_amplicon_aln_files : cannot create $tmpfile1\n";
        foreach $seq (sort {$a<=>$b} (keys(%$refH_sequences)))
        {
        	print INTFNA ">$seq\n$refH_sequences->{$seq}{'intergene'}\n";
        }
	close(INTFNA);

	$tmpfile2 = run_mafft_einsi( $tmpfile1 );
	if($tmpfile2 eq $phyTools::ERROR){ die "# cut_intergen_amplicon_aln_files : sorry, run_mafft_einsi failed\n"; } 
	%FASTA = read_FASTA_sequence( $tmpfile2 );
        foreach $seq (keys(%FASTA))
        {
 		$seqnozeros = $FASTA{$seq}{'NAME'};
               	$seqnozeros =~ s/^>|\n//g; 
		$intergene{$seqnozeros} = $FASTA{$seq}{'SEQ'};
		
		if($refH_sequences->{$seqnozeros}{'leftSTOP'})
		{
			my $pos=0;
			while($pos<length($intergene{$seqnozeros}))
			{
				last if(substr($intergene{$seqnozeros},$pos,1) ne '-');
				$pos++;
			}
			my $stop = uc(substr($intergene{$seqnozeros},$pos,3));
			my $gaps = '-' x $pos;
			$intergene{$seqnozeros} = $gaps . $stop . 
				substr($intergene{$seqnozeros},$pos+3);
		}	
		if($refH_sequences->{$seqnozeros}{'rightSTOP'})
                {
			my $pos=length($intergene{$seqnozeros})-1;
			while($pos>=0)
			{
				last if(substr($intergene{$seqnozeros},$pos,1) ne '-');
				$pos--;
			}

                        my $stop = uc(substr($intergene{$seqnozeros},$pos-2,3));
			my $gaps = '-' x (length($intergene{$seqnozeros})-1-$pos);
                        $intergene{$seqnozeros} = substr($intergene{$seqnozeros},0,$pos-2) .
				$stop . $gaps;
                }
        }
	unlink($tmpfile1,$tmpfile2);

	# process each pair of primers 
	foreach $score (sort {$b<=>$a} (keys(%$refH_amplicons)))
	{
        	PAIR: foreach $pair (keys(%{$refH_amplicons->{$score}}))
                {
			($left,$right) = split(/,/,$pair); 
			$revright = reverse($right);
			$N_CDScoord = $refH_amplicons->{$score}{$pair}[0]{'CDScoord'};
			$C_CDScoord = $refH_amplicons->{$score}{$pair}[1]{'CDScoord'};
			$N_primer_sense = $refH_amplicons->{$score}{$pair}[0]{'complement'};
			$C_primer_sense = $refH_amplicons->{$score}{$pair}[1]{'complement'};
			$take_complement = 0;

			foreach $seq (keys(%$refH_strands))
                        {			
				if($refH_strands->{$seq}{'left'} eq '-1'){ $N_strand = '-1' }
				if($refH_strands->{$seq}{'right'} eq '-1'){ $C_strand = '-1'; last; }
			}
			
			# put primers in genomic orientation
                        if($N_strand eq '-1')
                        {
				$take_complement = 1; # in our tests 2,6,10,39,47 only left-1 require this
				$left = reverse($left);
                        	if($N_primer_sense){ $left =~ tr/acgtwsrykmbvdhACGT/tgcaswyrmkvbhdTGCA/; }
                        }
                        if($C_strand eq '-1')
                        {
				$right = reverse($right);
				$revright = reverse($revright);
                        	if($C_primer_sense)
				{ 
					$right =~ tr/acgtwsrykmbvdhACGT/tgcaswyrmkvbhdTGCA/; 		
					$revright =~ tr/acgtwsrykmbvdhACGT/tgcaswyrmkvbhdTGCA/;
				}
				
                        }#print "# ($score) $N_CDScoord $left $N_primer_sense $C_CDScoord $right $C_primer_sense\n";

			# check left flank alignment
			my @left_align;
			foreach $seq (keys(%$refH_sequences))
			{
				#print "lf $refH_DNA_aln->{'left'}{$seq}\n"; # full length		
				$length = $refH_amplicons->{$score}{$pair}[0]{'Clength'} +
					$refH_amplicons->{$score}{$pair}[0]{'length'};

				$alignedCDS{'left'}{$seq} = substr($refH_DNA_aln->{'left'}{$seq},$N_CDScoord,$length);
				$amplicon = "$alignedCDS{'left'}{$seq}\n"; #print "la $amplicon\n";
				
				$unalignedCDS{'left'}{$seq} = substr($refH_DNA_aln->{'left'}{$seq},$N_CDScoord+$length);
				#print "lu $unalignedCDS{'left'}{$seq}\n";		

				# check primer genomic matching
				($n_of_matches, $amplicon_match_aln) = calc_matches($amplicon,$left,0,$take_complement); 	
				$N_prop_matches = sprintf("%1.1f",$n_of_matches / length($left));                                     
                                if($n_of_matches < sprintf("%.0f",length($left)*$ENV{"MIN_PROP_PRIMER_MATCHES"}))
	                        {
                               		print "# cut_intergen_amplicon_aln_files : left primer in amplicon $n_of_amplicons ".
						"matches poorly $N_prop_matches\n";# ($amplicon,$left)\n";
					next PAIR;
				}
				if(scalar(@left_align) == 0)
				{
					my $codehop_primer = $left;	
					$codehop_primer =~ tr/acgtACGT/tgcaTGCA/ if($take_complement); 
					push(@left_align,"$codehop_primer 5'->3'"); 
				}  
                                push(@left_align,"$amplicon_match_aln $refH_sequences->{$seq}{'NAME'}");
			}
			
			# check right flank alignment
			my @right_align;
                        foreach $seq (keys(%$refH_sequences))
                        {
				#print "rf $refH_DNA_aln->{'right'}{$seq}\n"; # full length               
                                $length = $refH_amplicons->{$score}{$pair}[1]{'Clength'} +
                                        $refH_amplicons->{$score}{$pair}[1]{'length'};

				$alignedCDS{'right'}{$seq} = substr($refH_DNA_aln->{'right'}{$seq},$C_CDScoord,$length);
                                $amplicon = reverse($alignedCDS{'right'}{$seq})."\n"; #print "ra $alignedCDS{'right'}{$seq}\n";
                                
                                $unalignedCDS{'right'}{$seq} = substr($refH_DNA_aln->{'right'}{$seq},0,$C_CDScoord);
                                #print "ru $unalignedCDS{'right'}{$seq}\n";  
				
				# check primer genomic matching, taking complement of primer
				($n_of_matches, $amplicon_match_aln) = calc_matches($amplicon,$revright,0,1);        
				$C_prop_matches = sprintf("%1.1f",$n_of_matches / length($revright));
                                if($n_of_matches < sprintf("%.0f",length($revright)*$ENV{"MIN_PROP_PRIMER_MATCHES"}))
                                {
                                        print "# cut_intergen_amplicon_aln_files : right primer in amplicon $n_of_amplicons ".
                                                "matches poorly $C_prop_matches\n";# ($amplicon,$revright)\n";
                                        next PAIR;
                                }
                                if(scalar(@right_align) == 0)
				{
					my $codehop_primer = $revright;
					$codehop_primer =~ tr/acgtACGT/tgcaTGCA/;
					push(@right_align,"$codehop_primer 5'->3'"); 
				}  
                                push(@right_align,"$amplicon_match_aln $refH_sequences->{$seq}{'NAME'}");
				$right_primer_length = $length;
                        }	

			# concatenate unaligned segments to intergenic sequences
			$DNA_aln_file_name = $root .'_'. $n_of_amplicons . '__'.
				$refH_amplicons->{$score}{$pair}[0]{'table'}.'-'.
				$refH_amplicons->{$score}{$pair}[1]{'table'}.'.dna_amp'; #print "# $DNA_aln_file_name\n";
			open(DNAAMP,">$DNA_aln_file_name") || 
				die "# cut_intergen_amplicon_aln_files : cannot create $DNA_aln_file_name\n";
			foreach $seq (keys(%$refH_sequences))
                        {
				$header = $refH_sequences->{$seq}{'NAME'}; 
				chomp($header);
				$header .= " #$N_CDScoord|$C_CDScoord|";
                                $header .=  "$refH_amplicons->{$score}{$pair}[0]{'complement'}|";
                                $header .=  "$refH_amplicons->{$score}{$pair}[1]{'complement'}|";
                                $header .=  "$refH_amplicons->{$score}{$pair}[0]{'offset'}|";
                                $header .=  "$refH_amplicons->{$score}{$pair}[1]{'offset'}|";
                                $header .=  "$refH_amplicons->{$score}{$pair}[0]{'temperature'}|";
                                $header .=  "$refH_amplicons->{$score}{$pair}[1]{'temperature'}|";
                                $header .=  "$left|$right|$N_prop_matches|$C_prop_matches";

				$sequence = $unalignedCDS{'left'}{$seq}.$intergene{$seq}.$unalignedCDS{'right'}{$seq};
				$left_intergene_aligned_length = length( $refH_DNA_aln->{'left'}{$seq}.$intergene{$seq} );
				$total_aligned_length = length( $refH_DNA_aln->{'left'}{$seq} . 
					$intergene{$seq} . $refH_DNA_aln->{'right'}{$seq} );
				print DNAAMP "$header\n$sequence\n";
                        }
                	close(DNAAMP);
			push(@amplicon_file_names,$DNA_aln_file_name);
			
			if($n_of_amplicons == 1)
			{
				print "# reference intergene coordinates: 1 .. $total_aligned_length ".
				"flank_sizes: $refH_lengthflank->{'left'} $refH_lengthflank->{'right'} ".
				"strands: $N_strand $C_strand ".
				"\n\n";
			}

			# calculate amplicon boundary coordinates over aligned nt sequences
			$N = $N_CDScoord;
			$C = $left_intergene_aligned_length + $C_CDScoord + $right_primer_length;
			print "# amplicon coords: $N $C\n";
				
			# print primer alignment blocks
                        $left_align[0] .= " N $N $C (aligned nucleotides)\n";
			$right_align[0] .= " C $N $C (aligned nucleotides)\n";
			$primer_filename = $root .'_'. $n_of_amplicons . '__'. 
                                $refH_amplicons->{$score}{$pair}[0]{'table'}.'-'.
                                $refH_amplicons->{$score}{$pair}[1]{'table'}.'.primers';
                        open(PRIMERS,">$primer_filename") 
                                || die "# cut_intergen_amplicon_aln_files : cannot create $primer_filename\n";
                        print PRIMERS "\n## Amplicon $root $n_of_amplicons:\n";
                        foreach $amplicon_match_aln (check_primer_match($root,$n_of_amplicons,'N',@left_align))
                        {
                                print PRIMERS $amplicon_match_aln;
                        }
                        print PRIMERS "\n";
			foreach $amplicon_match_aln (check_primer_match($root,$n_of_amplicons,'C',@right_align))
                        {
                                print PRIMERS $amplicon_match_aln;
                        }
                        close(PRIMERS);

			push(@primers_file_names,$primer_filename);	
			$n_of_amplicons++;
                }
        }	

	return (\@primers_file_names,\@amplicon_file_names);
}


sub cut_amplicon_aln_files
{
	my ($DNA_aln_file,$AA_aln_file,$codon_table_name,$print_amplicons,%amplicons) = @_;
	my (@amplicon_file_names,$pos,$N,$C,$score,$primerN,$primerC);
	my (%AA_aln,%DNA_aln,$seq,$amplicon,$header,$prop_matches);
	my ($n_of_amplicons,$n_of_matches,$primerpos,$primersOK);
	my ($AA_aln_file_name,$DNA_aln_file_name,$aa_filename,$dna_filename);
	my ($amplicon_match_aln,@primers_file_names,$primer_filename);
	
       	$codon_table_name = (split(/__/,$codon_table_name))[1];
	$codon_table_name =~ s/\.codon\.use//;
	
	# 1) read FASTA sequences
	%AA_aln  = read_FASTA_sequence($AA_aln_file);
	%DNA_aln = read_FASTA_sequence($DNA_aln_file);
	
	# 2) define filenames for each amplicon file
	$AA_aln_file_name  = (split(/\./,$AA_aln_file))[0];
	$DNA_aln_file_name = (split(/\./,$DNA_aln_file))[0];
	
	#print "$AA_aln_file\n$DNA_aln_file\n";
	
	# 3) loop through selected amplicons cutting aligned sequences 
	$n_of_amplicons = 0;
	foreach $pos (sort {$a<=>$b} (keys(%amplicons)))
	{
		foreach $score (keys(%{$amplicons{$pos}}))
		{
			foreach $primerN (keys(%{$amplicons{$pos}{$score}}))
			{
				foreach $primerC (keys(%{$amplicons{$pos}{$score}{$primerN}}))
				{
					my (@AAFASTA,@DNAFASTA);
					my (@N_align,@C_align);
					$primersOK=1;
					
					$N = $amplicons{$pos}{$score}{$primerN}{$primerC}[0]{'AA_align_coord'};
					$C = $amplicons{$pos}{$score}{$primerN}{$primerC}[1]{'AA_align_coord'};			
					$n_of_amplicons++;
					#print "$N $C $n_of_amplicons\n";

					# 3.1) store AA amplicons in FASTA format
					foreach $seq (sort {$a<=>$b} (keys(%AA_aln)))
					{
						$amplicon = substr($AA_aln{$seq}{'SEQ'},$N,$C-$N+2+(length($primerC)/3));
						#$amplicon =~ s/-//g;
						$header = $AA_aln{$seq}{'NAME'}; chomp $header;
						$header .= " #$N|$C|";
						$header .=  "$amplicons{$pos}{$score}{$primerN}{$primerC}[0]{'complement'}|";
						$header .=  "$amplicons{$pos}{$score}{$primerN}{$primerC}[1]{'complement'}|";
						$header .=  "$amplicons{$pos}{$score}{$primerN}{$primerC}[0]{'offset'}|";
						$header .=  "$amplicons{$pos}{$score}{$primerN}{$primerC}[1]{'offset'}|";
						$header .=  "$amplicons{$pos}{$score}{$primerN}{$primerC}[0]{'temperature'}|";
						$header .=  "$amplicons{$pos}{$score}{$primerN}{$primerC}[1]{'temperature'}|";
						$header .=  "$primerN|$primerC";			
						push(@AAFASTA,"$header\n$amplicon\n");
					}
					
					# 3.2) store DNA amplicons in FASTA format
					foreach $seq (sort {$a<=>$b} (keys(%DNA_aln)))
					{
						$amplicon = substr($DNA_aln{$seq}{'SEQ'},$N*3,($C*3)-($N*3)+4+length($primerC));
						#$amplicon =~ s/-//g;
						$header = $DNA_aln{$seq}{'NAME'}; chomp $header;
						$header .= " #$N|$C|";
						$header .=  "$amplicons{$pos}{$score}{$primerN}{$primerC}[0]{'complement'}|";
						$header .=  "$amplicons{$pos}{$score}{$primerN}{$primerC}[1]{'complement'}|";
						$header .=  "$amplicons{$pos}{$score}{$primerN}{$primerC}[0]{'offset'}|";
						$header .=  "$amplicons{$pos}{$score}{$primerN}{$primerC}[1]{'offset'}|";
						$header .=  "$amplicons{$pos}{$score}{$primerN}{$primerC}[0]{'temperature'}|";
						$header .=  "$amplicons{$pos}{$score}{$primerN}{$primerC}[1]{'temperature'}|";
						$header .=  "$primerN|$primerC|";			
						
						($n_of_matches, $amplicon_match_aln) = 
							calc_matches($amplicon,$primerN,3+$amplicons{$pos}{$score}{$primerN}{$primerC}[0]{'offset'},
							$amplicons{$pos}{$score}{$primerN}{$primerC}[0]{'complement'});
						
						if(scalar(@N_align) == 0)
						{
							#push(@N_align,"$primerN 5'->3' N $N $C $amplicons{$pos}{$score}{$primerN}{$primerC}[0]{'temperature'}\n");
							push(@N_align,"$primerN 5'->3' N $N $C (aligned residues)\n");
						}
						push(@N_align,"$amplicon_match_aln $DNA_aln{$seq}{'NAME'}");
						
						$prop_matches = sprintf("%1.1f",$n_of_matches / length($primerN));																	
																										
						if($n_of_matches < sprintf("%.0f",length($primerN)*$ENV{"MIN_PROP_PRIMER_MATCHES"}))
						{
							print "# cut_amplicon_aln_files : primer N in amplicon $n_of_amplicons matches poorly $prop_matches\n";
							$primersOK=0;
						}									
						$header .= "$prop_matches|";
					
						$primerpos = ($C*3)-($N*3)+3;
						($n_of_matches,$amplicon_match_aln) = 
							calc_matches($amplicon,$primerC,
								$primerpos+$amplicons{$pos}{$score}{$primerN}{$primerC}[1]{'offset'},
								$amplicons{$pos}{$score}{$primerN}{$primerC}[1]{'complement'});
						
						if(scalar(@C_align) == 0)
						{
							#push(@C_align,reverse($primerC)." 5'->3' C $N $C $amplicons{$pos}{$score}{$primerN}{$primerC}[1]{'temperature'}\n");
							push(@C_align,reverse($primerC)." 5'->3' C $N $C (aligned residues)\n");
						}
						$amplicon_match_aln =~ tr/atgcATGC/tacgTACG/; 
						push(@C_align,reverse($amplicon_match_aln)." $DNA_aln{$seq}{'NAME'}");
						
						$prop_matches = sprintf("%1.1f",$n_of_matches / length($primerC));
						
						if($n_of_matches < sprintf("%.0f",length($primerC)*$ENV{"MIN_PROP_PRIMER_MATCHES"}))
						{
							#printf("#%d#%.0f#\n", $n_of_matches, length($primerC)*$ENV{"MIN_PROP_PRIMER_MATCHES"});
							print "# cut_amplicon_aln_files : primer C in amplicon $n_of_amplicons matches poorly $prop_matches\n";
							$primersOK=0;
						}	
						$header .= "$prop_matches|";
													
						push(@DNAFASTA,"$header\n$amplicon\n");
					}
				
					if($primersOK)
					{
						# 3.3) print primer alignment blocks
						$primer_filename = $AA_aln_file_name ."_" . $n_of_amplicons . "__" . $codon_table_name . ".primers";
						open(PRIMERS,">$primer_filename") || die "# cut_amplicon_aln_files 3.3 : cannot create $primer_filename\n";
						
						print PRIMERS "\n## Amplicon $AA_aln_file_name $n_of_amplicons:\n";
						
						foreach $amplicon_match_aln (check_primer_match($AA_aln_file_name,$n_of_amplicons,'N',@N_align))
						{
							print PRIMERS $amplicon_match_aln;
						}
						print PRIMERS "\n";
						foreach $amplicon_match_aln (check_primer_match($AA_aln_file_name,$n_of_amplicons,'C',@C_align))
						{
							print PRIMERS $amplicon_match_aln;
						}
						close(PRIMERS); #print "\n###################################################################################\n";
						
						push(@primers_file_names,$primer_filename);
						
						# 3.4) write AA & DNA files and store their filenames
						if($print_amplicons)
						{
							$aa_filename = $AA_aln_file_name ."_".$n_of_amplicons. "__" . $codon_table_name .".aa_amp";
							open(AAFAS,">$aa_filename") || die "# cut_amplicon_aln_files : cannot create $aa_filename\n";
							print AAFAS @AAFASTA;
							close(AAFAS);
					
							$dna_filename = $DNA_aln_file_name ."_".$n_of_amplicons. "__" . $codon_table_name .".dna_amp";
							open(DNAFAS,">$dna_filename") || die "# cut_amplicon_aln_files : cannot create $dna_filename\n";
							print DNAFAS @DNAFASTA;
							close(DNAFAS);
									
							push(@amplicon_file_names,$aa_filename ."|". $dna_filename);
						}
					}
				}
			}
		}	
	}
	
	return (\@primers_file_names,\@amplicon_file_names);
}	
	
sub check_primer_match
{
	my ($name,$ampnumber,$NC,@aligned_matches) = @_; #print "# @aligned_matches\n";
	my ($primer,$dnaseq,$n_of_sequences,$length,$s,$p,$domOK,$degen_start);
	my ($nt,$dnt,$primernt,%freq,$qualityline,$quality,$degen);
	my ($corrected_primer,$corrected_degen_primer,$corrected_relaxed_primer);
	
	$n_of_sequences = scalar(@aligned_matches) - 1;
	$primer = (split(/ /,$aligned_matches[0]))[0];
	$length = length($primer);
	$qualityline = $corrected_primer = $corrected_degen_primer = $corrected_relaxed_primer = '';
	
	for($p=0;$p<$length;$p++)
	{
		if(substr($primer,$p,1) =~ /[a-z]/){ $degen_start = $p; last; } 
	}
	
	# 1) first compute alignment stats
	for($p=0;$p<$length;$p++)
	{
		$freq{$p}{'A'} = $freq{$p}{'T'} = $freq{$p}{'C'} = $freq{$p}{'G'} = 0; 
	}
	
	for($p=0;$p<$length;$p++)
	{
		for($s=1;$s<$n_of_sequences+1;$s++)
		{
			$nt = uc(substr($aligned_matches[$s],$p,1));
			$freq{$p}{$nt}++;
		}
	}
	
	# 2) check each column $p from left to right
	for($p=0;$p<$length;$p++)
	{
		$primernt = uc(substr($primer,$p,1));
		$nt = $primernt;
		$quality = '.';
		$domOK = '';
	
		# 100% conserved columns
		if($freq{$p}{'A'} == $n_of_sequences || $freq{$p}{'T'} == $n_of_sequences ||
			$freq{$p}{'C'} == $n_of_sequences || $freq{$p}{'G'} == $n_of_sequences)
		{
			#if(!$aligned_matches[1]){ print @aligned_matches; exit; }
			if(uc(substr($aligned_matches[1],$p,1)) ne $primernt)
			{ 
				$quality = '!'; 
				$dnt = lc(substr($aligned_matches[1],$p,1));
				$nt = uc(substr($aligned_matches[1],$p,1)); 
				if($p >= $degen_start)
				{ 
					$nt = lc(substr($aligned_matches[1],$p,1)); 
				}
			}
		} 
		elsif($freq{$p}{'C'}/$n_of_sequences> $ENV{"MIN_FREQ_CONSENSUS"}){ $domOK = 'C' }
		elsif($freq{$p}{'G'}/$n_of_sequences> $ENV{"MIN_FREQ_CONSENSUS"}){ $domOK = 'G' }
		elsif($freq{$p}{'A'}/$n_of_sequences> $ENV{"MIN_FREQ_CONSENSUS"}){ $domOK = 'A' }
		elsif($freq{$p}{'T'}/$n_of_sequences> $ENV{"MIN_FREQ_CONSENSUS"}){ $domOK = 'T' }
		
		# dominated columns (with freq > $ENV{"MIN_FREQ_CONSENSUS"})
		if($domOK ne '' && $p < $degen_start)
		{
			if($primernt ne $domOK)
			{
				$quality = '?'; 
				if($p >= $degen_start){ $nt = $domOK; }
				$dnt = $domOK;
			}	
		}
		
		if($freq{$p}{'C'} > 0 && $freq{$p}{'T'} > 0 && $freq{$p}{'A'} == 0 && $freq{$p}{'G'} == 0)
		{
			if($primernt ne 'Y'){ $quality = '?' }
			if($p >= $degen_start){ $nt = 'y'; }
			$dnt = 'y';
		}
		elsif($freq{$p}{'C'} == 0 && $freq{$p}{'T'} == 0 && $freq{$p}{'A'} > 0 && $freq{$p}{'G'} > 0)
		{
			if($primernt ne 'R'){ $quality = '?' }
			if($p >= $degen_start){ $nt = 'r'; }
			$dnt = 'r';
		}
		elsif($freq{$p}{'C'} > 0 && $freq{$p}{'T'} == 0 && $freq{$p}{'A'} > 0 && $freq{$p}{'G'} == 0)
		{
			if($primernt ne 'M'){ $quality = '?' }
			if($p >= $degen_start){ $nt = 'm'; }
			$dnt = 'm';	
		}
		elsif($freq{$p}{'C'} == 0 && $freq{$p}{'T'} > 0 && $freq{$p}{'A'} > 0 && $freq{$p}{'G'} == 0)
		{
			if($primernt ne 'W'){ $quality = '?' }
			if($p >= $degen_start){ $nt = 'w'; }
			$dnt = 'w';
		}
		elsif($freq{$p}{'C'} > 0 && $freq{$p}{'T'} == 0 && $freq{$p}{'A'} == 0 && $freq{$p}{'G'} > 0)
		{
			if($primernt ne 'S'){ $quality = '?' }
			if($p >= $degen_start){ $nt = 's'; }
			$dnt = 's';
		}
		elsif($freq{$p}{'C'} == 0 && $freq{$p}{'T'} > 0 && $freq{$p}{'A'} == 0 && $freq{$p}{'G'} > 0)
		{
			if($primernt ne 'K'){ $quality = '?' }
			if($p >= $degen_start){ $nt = 'k'; }
			$dnt = 'k';
		}
		elsif($freq{$p}{'C'} == 0 && $freq{$p}{'T'} > 0 && $freq{$p}{'A'} > 0 && $freq{$p}{'G'} > 0)
		{
			if($primernt ne 'D'){ $quality = '?' }
			if($p >= $degen_start){ $nt = 'd'; }
			$dnt = 'd';
		}
		elsif($freq{$p}{'C'} > 0 && $freq{$p}{'T'} > 0 && $freq{$p}{'A'} > 0 && $freq{$p}{'G'} == 0)
		{
			if($primernt ne 'H'){ $quality = '?' }
			if($p >= $degen_start){ $nt = 'h'; }
			$dnt = 'h';
		}
		elsif($freq{$p}{'C'} > 0 && $freq{$p}{'T'} > 0 && $freq{$p}{'A'} == 0 && $freq{$p}{'G'} > 0)
		{
			if($primernt ne 'B'){ $quality = '?' }
			if($p >= $degen_start){ $nt = 'b'; }
			$dnt = 'b';
		}
		elsif($freq{$p}{'C'} > 0 && $freq{$p}{'T'} == 0 && $freq{$p}{'A'} > 0 && $freq{$p}{'G'} > 0)
		{
			if($primernt ne 'V'){ $quality = '?' }
			if($p >= $degen_start){ $nt = 'v'; }
			$dnt = 'v';
		}
		elsif($freq{$p}{'C'} > 0 && $freq{$p}{'T'} > 0 && $freq{$p}{'A'} > 0 && $freq{$p}{'G'} > 0)
		{
			if($primernt ne 'N'){ $quality = '?' }
			if($p >= $degen_start){ $nt = 'n'; }
			$dnt = 'n';
		}
		elsif($quality ne '!')
		{
			#print "hola |$nt|$dnt|\n";
			if($p >= $degen_start){ $nt = lc($nt); }
			$dnt = lc($nt);
		}
		
		$qualityline .= $quality;
		$corrected_primer .= $nt;
		$corrected_degen_primer .= $dnt;
	}
	
	# calculate relaxed primer (Pablo's brilliant idea of 22032007)

	$corrected_relaxed_primer = $corrected_primer;
	$corrected_relaxed_primer =~ s/[AGCT]//g;

	$p = $degen_start-1;
	while($p>=0) 
        {
		$dnt = substr($corrected_degen_primer,$p,1);
                $nt = substr($corrected_primer,$p,1);

		$degen = calc_oligo_degeneracy($dnt . $corrected_relaxed_primer);
		if($degen < $ENV{"MAX_RELAXED_DEGEN"})
		{
			$corrected_relaxed_primer = $dnt . $corrected_relaxed_primer;
			$p--;
		}
		else
		{ 
			$corrected_relaxed_primer = substr($corrected_primer,0,$p+1) . $corrected_relaxed_primer;
			last;
		}
	}
	
	#print "#$corrected_primer\n";
	#print "#$corrected_degen_primer\n";
	#print "#$corrected_relaxed_primer\n";	
	#exit;

	$name =~ s/my_//;
	$name =~ s/\_aln/\_amp/;
	if($NC eq 'C')
	{
		$corrected_primer .= " codeh_corr $name$ampnumber\_C" . (split(/ /,$aligned_matches[0]))[4] . "\n"; 
		$corrected_relaxed_primer .= " relax_corr $name$ampnumber\_C" . (split(/ /,$aligned_matches[0]))[4] . "\n";
		$corrected_degen_primer .= " degen_corr $name$ampnumber\_C" . (split(/ /,$aligned_matches[0]))[4] . "\n";
	}
	else
	{
		$corrected_primer .= " codeh_corr $name$ampnumber\_N" . (split(/ /,$aligned_matches[0]))[3] . "\n";
		$corrected_relaxed_primer .= " relax_corr $name$ampnumber\_N" . (split(/ /,$aligned_matches[0]))[3] . "\n"; 
		$corrected_degen_primer .= " degen_corr $name$ampnumber\_N" . (split(/ /,$aligned_matches[0]))[3] . "\n";
	}
	
	push(@aligned_matches,"$qualityline\n");	
	push(@aligned_matches,$corrected_primer);
	push(@aligned_matches,$corrected_relaxed_primer);
	push(@aligned_matches,$corrected_degen_primer);	

	return @aligned_matches;
}


sub check_primer_match_excluded
{
	my ($name,$ampnumber,$NC,$ref_aligned_matches,$ref_excluded_matches) = @_; 
	my ($primer,$dnaseq,$n_of_sequences,$length,$s,$p,$domOK,$degen_start);
	my ($nt,$dnt,$primernt,%freq,$qualityline,$quality,$degen);
	my ($corrected_primer,$corrected_degen_primer,$corrected_relaxed_primer);
	my @aligned_matches = @{$ref_aligned_matches};	

	$n_of_sequences = scalar(@aligned_matches) - 1;
	$primer = (split(/ /,$aligned_matches[0]))[0];
	$length = length($primer);
	$qualityline = $corrected_primer = $corrected_degen_primer = $corrected_relaxed_primer = '';
	
	for($p=0;$p<$length;$p++)
	{
		if(substr($primer,$p,1) =~ /[a-z]/){ $degen_start = $p; last; } 
	}
	
	# 1) first compute alignment stats
	for($p=0;$p<$length;$p++)
	{
		$freq{$p}{'A'} = $freq{$p}{'T'} = $freq{$p}{'C'} = $freq{$p}{'G'} = 0; 
	}
	
	for($p=0;$p<$length;$p++)
	{
		for($s=1;$s<$n_of_sequences+1;$s++)
		{
			$nt = uc(substr($aligned_matches[$s],$p,1));
			$freq{$p}{$nt}++;
		}
	}
	
	# 2) check each column $p from left to right
	for($p=0;$p<$length;$p++)
	{
		$primernt = uc(substr($primer,$p,1));
		$nt = $primernt;
		$quality = '.';
		$domOK = '';
	
		# 100% conserved columns
		if($freq{$p}{'A'} == $n_of_sequences || $freq{$p}{'T'} == $n_of_sequences ||
			$freq{$p}{'C'} == $n_of_sequences || $freq{$p}{'G'} == $n_of_sequences)
		{
			if(uc(substr($aligned_matches[1],$p,1)) ne $primernt)
			{ 
				$quality = '!'; 
				$dnt = lc(substr($aligned_matches[1],$p,1));
				$nt = uc(substr($aligned_matches[1],$p,1)); 
				if($p >= $degen_start)
				{ 
					$nt = lc(substr($aligned_matches[1],$p,1)); 
				}
			}
		} 
		elsif($freq{$p}{'C'}/$n_of_sequences> $ENV{"MIN_FREQ_CONSENSUS"}){ $domOK = 'C' }
		elsif($freq{$p}{'G'}/$n_of_sequences> $ENV{"MIN_FREQ_CONSENSUS"}){ $domOK = 'G' }
		elsif($freq{$p}{'A'}/$n_of_sequences> $ENV{"MIN_FREQ_CONSENSUS"}){ $domOK = 'A' }
		elsif($freq{$p}{'T'}/$n_of_sequences> $ENV{"MIN_FREQ_CONSENSUS"}){ $domOK = 'T' }
		
		# dominated columns (with freq > $ENV{"MIN_FREQ_CONSENSUS"})
		if($domOK ne '' && $p < $degen_start)
		{
			if($primernt ne $domOK)
			{
				$quality = '?'; 
				if($p >= $degen_start){ $nt = $domOK; }
				$dnt = $domOK;
			}	
		}
		
		if($freq{$p}{'C'} > 0 && $freq{$p}{'T'} > 0 && $freq{$p}{'A'} == 0 && $freq{$p}{'G'} == 0)
		{
			if($primernt ne 'Y'){ $quality = '?' }
			if($p >= $degen_start){ $nt = 'y'; }
			$dnt = 'y';
		}
		elsif($freq{$p}{'C'} == 0 && $freq{$p}{'T'} == 0 && $freq{$p}{'A'} > 0 && $freq{$p}{'G'} > 0)
		{
			if($primernt ne 'R'){ $quality = '?' }
			if($p >= $degen_start){ $nt = 'r'; }
			$dnt = 'r';
		}
		elsif($freq{$p}{'C'} > 0 && $freq{$p}{'T'} == 0 && $freq{$p}{'A'} > 0 && $freq{$p}{'G'} == 0)
		{
			if($primernt ne 'M'){ $quality = '?' }
			if($p >= $degen_start){ $nt = 'm'; }
			$dnt = 'm';	
		}
		elsif($freq{$p}{'C'} == 0 && $freq{$p}{'T'} > 0 && $freq{$p}{'A'} > 0 && $freq{$p}{'G'} == 0)
		{
			if($primernt ne 'W'){ $quality = '?' }
			if($p >= $degen_start){ $nt = 'w'; }
			$dnt = 'w';
		}
		elsif($freq{$p}{'C'} > 0 && $freq{$p}{'T'} == 0 && $freq{$p}{'A'} == 0 && $freq{$p}{'G'} > 0)
		{
			if($primernt ne 'S'){ $quality = '?' }
			if($p >= $degen_start){ $nt = 's'; }
			$dnt = 's';
		}
		elsif($freq{$p}{'C'} == 0 && $freq{$p}{'T'} > 0 && $freq{$p}{'A'} == 0 && $freq{$p}{'G'} > 0)
		{
			if($primernt ne 'K'){ $quality = '?' }
			if($p >= $degen_start){ $nt = 'k'; }
			$dnt = 'k';
		}
		elsif($freq{$p}{'C'} == 0 && $freq{$p}{'T'} > 0 && $freq{$p}{'A'} > 0 && $freq{$p}{'G'} > 0)
		{
			if($primernt ne 'D'){ $quality = '?' }
			if($p >= $degen_start){ $nt = 'd'; }
			$dnt = 'd';
		}
		elsif($freq{$p}{'C'} > 0 && $freq{$p}{'T'} > 0 && $freq{$p}{'A'} > 0 && $freq{$p}{'G'} == 0)
		{
			if($primernt ne 'H'){ $quality = '?' }
			if($p >= $degen_start){ $nt = 'h'; }
			$dnt = 'h';
		}
		elsif($freq{$p}{'C'} > 0 && $freq{$p}{'T'} > 0 && $freq{$p}{'A'} == 0 && $freq{$p}{'G'} > 0)
		{
			if($primernt ne 'B'){ $quality = '?' }
			if($p >= $degen_start){ $nt = 'b'; }
			$dnt = 'b';
		}
		elsif($freq{$p}{'C'} > 0 && $freq{$p}{'T'} == 0 && $freq{$p}{'A'} > 0 && $freq{$p}{'G'} > 0)
		{
			if($primernt ne 'V'){ $quality = '?' }
			if($p >= $degen_start){ $nt = 'v'; }
			$dnt = 'v';
		}
		elsif($freq{$p}{'C'} > 0 && $freq{$p}{'T'} > 0 && $freq{$p}{'A'} > 0 && $freq{$p}{'G'} > 0)
		{
			if($primernt ne 'N'){ $quality = '?' }
			if($p >= $degen_start){ $nt = 'n'; }
			$dnt = 'n';
		}
		elsif($quality ne '!')
		{
			#print "hola |$nt|$dnt|\n";
			if($p >= $degen_start){ $nt = lc($nt); }
			$dnt = lc($nt);
		}
		
		$qualityline .= $quality;
		$corrected_primer .= $nt;
		$corrected_degen_primer .= $dnt;
	}
	
	# calculate relaxed primer (Pablo's brilliant idea of 22032007)

	$corrected_relaxed_primer = $corrected_primer;
	$corrected_relaxed_primer =~ s/[AGCT]//g;

	$p = $degen_start-1;
	while($p>=0) 
        {
		$dnt = substr($corrected_degen_primer,$p,1);
                $nt = substr($corrected_primer,$p,1);

		$degen = calc_oligo_degeneracy($dnt . $corrected_relaxed_primer);
		if($degen < $ENV{"MAX_RELAXED_DEGEN"})
		{
			$corrected_relaxed_primer = $dnt . $corrected_relaxed_primer;
			$p--;
		}
		else
		{ 
			$corrected_relaxed_primer = substr($corrected_primer,0,$p+1) . $corrected_relaxed_primer;
			last;
		}
	}
	
	$name =~ s/my_//;
	$name =~ s/\_aln/\_amp/;
	if($NC eq 'C')
	{
		$corrected_primer .= " codeh_corr $name$ampnumber\_C" . (split(/ /,$aligned_matches[0]))[4] . "\n"; 
		$corrected_relaxed_primer .= " relax_corr $name$ampnumber\_C" . (split(/ /,$aligned_matches[0]))[4] . "\n";
		$corrected_degen_primer .= " degen_corr $name$ampnumber\_C" . (split(/ /,$aligned_matches[0]))[4] . "\n";
	}
	else
	{
		$corrected_primer .= " codeh_corr $name$ampnumber\_N" . (split(/ /,$aligned_matches[0]))[3] . "\n";
		$corrected_relaxed_primer .= " relax_corr $name$ampnumber\_N" . (split(/ /,$aligned_matches[0]))[3] . "\n"; 
		$corrected_degen_primer .= " degen_corr $name$ampnumber\_N" . (split(/ /,$aligned_matches[0]))[3] . "\n";
	}
	
	push(@aligned_matches,"$qualityline\n");	
	push(@aligned_matches,$corrected_primer);
	push(@aligned_matches,$corrected_relaxed_primer);
	push(@aligned_matches,$corrected_degen_primer);	

	## enumerate mismatches between corrected  primer and excluded sequences
	foreach my $primer (@{$ref_excluded_matches})
	{
		push(@aligned_matches,calc_2oligo_mismatches($corrected_relaxed_primer,$primer,1));
	}

	return @aligned_matches;
}

sub calc_2oligo_mismatches
{
	my ($ref_oligo,$oligo,$only_upper) = @_;

	my ($pos,$report,$taxon,$length,$refnt,$refntorig,$repnt,$nt) = (0,'','');
	my ($n_of_mismatches5,$n_of_mismatches3) = (0,0);

	$ref_oligo = (split(/ /,$ref_oligo))[0];
	($oligo,$taxon) = split(/ /,$oligo,2);	
	$length = length($ref_oligo);
	if($only_upper)
	{ 
		while(substr($ref_oligo,$pos,1) =~ /[A-Z]/)
		{ 
			$pos++; 
			$report .= '-';
		} 
	}
	
	while($pos < $length)
        {
                $nt = substr($oligo,$pos,1);
		$refntorig = substr($ref_oligo,$pos,1);
		$refnt = lc($refntorig);

		if($refnt eq $nt || $refnt eq 'n'){ $report .= '-' }
		else
		{
			$repnt = '';
			if($refnt eq 'r' && ($nt ne 'a' && $nt ne 'g')){ $repnt = $nt; }
			elsif($refnt eq 'y' && ($nt ne 'c' && $nt ne 't')){ $repnt = $nt; }
			elsif($refnt eq 's' && ($nt ne 'c' && $nt ne 'g')){ $repnt = $nt; }
			elsif($refnt eq 'w' && ($nt ne 'a' && $nt ne 't')){ $repnt = $nt; }
			elsif($refnt eq 'm' && ($nt ne 'a' && $nt ne 'c')){ $repnt = $nt; }
			elsif($refnt eq 'k' && ($nt ne 't' && $nt ne 'g')){ $repnt = $nt; }
			elsif($refnt eq 'v' && $nt eq 't'){ $repnt = $nt; }
			elsif($refnt eq 'h' && $nt eq 'g'){ $repnt = $nt; }
			elsif($refnt eq 'd' && $nt eq 'c'){ $repnt = $nt; }
			elsif($refnt eq 'b' && $nt eq 'a'){ $repnt = $nt; }
                	else{ $report .= '-' }
		
			if($repnt)
			{
				if($refntorig =~ m/[A-Z]/){ $n_of_mismatches5++; $repnt = uc($repnt); }
				else{ $n_of_mismatches3++; }
				$report .= $repnt;
			}
		}
                $pos++;
        }

	#return sprintf("%s %d %d %s",$report,$n_of_mismatches5,$n_of_mismatches3,$taxon);
	return sprintf("%s %s",$report,$taxon);
}

sub calc_oligo_degeneracy
{
	my ($oligo) = @_;

	my ($pos,$nt,$length,$degen) = (0,'',length($oligo),1);

	$oligo = lc($oligo);

	while($pos < $length)
	{
		$nt = substr($oligo,$pos,1);
		if($nt =~ /[ryswmk]/){ $degen *= 2 }
		elsif($nt =~ /[vhdb]/){ $degen *= 3 }
		elsif($nt =~ /n/){ $degen *= 4 }
		$pos++;
	}

	return $degen;
}

sub calc_matches
{
	my ($amplicon,$primer,$pos,$complement) = @_;
	my ($n_of_matches,$length) = (0,0);
	my ($b,$p,$dna);
	
	# 1) locate hybridization position
	$amplicon = lc(substr($amplicon,$pos,length($primer)));
	$primer = lc($primer);
	
	#print "> $complement\n> $amplicon\n> $primer\n";
	
	# 2) calculate matchings
	if($complement)
	{ 
		$primer =~ tr/atgc/tacg/; 
		#print "> $primer\n";
	}
	
	$length = length($primer);
	for($b=0;$b<$length;$b++)
	{
		$p = substr($primer,$b,1);
		$dna = substr($amplicon,$b,1);
		
		if($p eq 'n'){ $n_of_matches++; }
		elsif($p eq 'v' && ($dna eq 'g' || $dna eq 'a' || $dna eq 'c')){ $n_of_matches++; }
		elsif($p eq 'b' && ($dna eq 'g' || $dna eq 't' || $dna eq 'c')){ $n_of_matches++; }
		elsif($p eq 'h' && ($dna eq 'a' || $dna eq 't' || $dna eq 'c')){ $n_of_matches++; }
		elsif($p eq 'd' && ($dna eq 'g' || $dna eq 'a' || $dna eq 't')){ $n_of_matches++; }
		elsif($p eq 'k' && ($dna eq 'g' || $dna eq 't')){ $n_of_matches++; }
		elsif($p eq 's' && ($dna eq 'g' || $dna eq 'c')){ $n_of_matches++; }
		elsif($p eq 'w' && ($dna eq 'a' || $dna eq 't')){ $n_of_matches++; }
		elsif($p eq 'm' && ($dna eq 'a' || $dna eq 'c')){ $n_of_matches++; }
		elsif($p eq 'y' && ($dna eq 'c' || $dna eq 't')){ $n_of_matches++; }
		elsif($p eq 'r' && ($dna eq 'a' || $dna eq 'g')){ $n_of_matches++; }
		elsif($p eq $dna) { $n_of_matches++; }
	}
	
	return ($n_of_matches,lc($amplicon));
}	 
																
sub run_revtrans
{
   my($infile_DNA,$infile_AA) = @_;
	
   my $outfile = (split(/\./,$infile_DNA))[0] . "_aln.fas";

   # run revtrans as we want OJO: igual conviene a?dir alguna opci? sobre el c?igo gen?ico 11
   system("$ENV{'EXE_PYTHON'} $ENV{'EXE_REVTRANS'} -match pos $infile_DNA $infile_AA > $outfile ");

   # check the output file is there
   if(! -s $outfile){ return $phyTools::ERROR; }
   else{ return $outfile; }
}

sub convert_FAS2PHY
{
    my($infile,$clean) = @_;

    my $tmpfile = (split(/\./,$infile))[0] . "_tmp_convertFAS2PHY";
    my $outfile = (split(/\./,$infile))[0] . ".phy";

    #print "#$infile#$outfile#\n";

    # run clustal to convert these formats
    system("$ENV{'EXE_CLUSTAL'} -infile=$infile -convert -output=phylip > $tmpfile");

    # clean file
    if($clean==1){ system(" rm -f $tmpfile"); }

    # check the output file is there
    if(! -s $outfile){ return $phyTools::ERROR; }
    else{ return $outfile; }
}

sub convert_FAS2NXS
{
    my($infile,$clean) = @_;

    my $tmpfile = (split(/\./,$infile))[0] . "_tmp_convertFAS2PHY";
    my $outfile = (split(/\./,$infile))[0] . ".nxs";

    # run clustal to convert these formats
    system("$ENV{'EXE_CLUSTAL'} -infile=$infile -convert -output=nexus > $tmpfile");

    # clean file
    if($clean==1){ system(" rm -f $tmpfile"); }

    # check the output file is there
    if(! -s $outfile){ return $phyTools::ERROR; }
    else{ return $outfile; }
}

sub run_MODELTEST
{
    my($infile,$clean) = @_;
 
    my $tmpfile = (split(/\./,$infile))[0] . "_tmp_PaupModeltest";
    my $outfile = (split(/\./,$infile))[0] . "_modeltest.out";
    my $bestAICmodel = $phyTools::ERROR;
    my $file_specific_modelblock = (split(/\./,$infile))[0] . "_modelblock";    
    my $file_specific_model_log = (split(/\./,$infile))[0] . "_model_log";    
    my $model_scores_file = (split(/\./,$infile))[0] . "_model.scores";    

    open (MODELBLOCK, "> $file_specific_modelblock") || die "cannot write $file_specific_modelblock ";
    open (TEMPLATE, "$ENV{'FILE_MODELBLOCK'}") || die "cannot open $ENV{'FILE_MODELBLOCK'} ";  
    while(<TEMPLATE>)
    {
	$_ =~ s/model\.scores/$model_scores_file/g;
	$_ =~ s/modelfit\.log/$file_specific_model_log/g;
	print MODELBLOCK; 
    } 
    close(TEMPLATE);	    

    # run PAUP* with the file-specific modelblocks for model selection
    
    if($ENV{GDMSESSION} && $ENV{GDMSESSION} eq 'gnome') # May2015 test var exists before using it
    { 
        $ENV{'EXE_PAUP'}='TERM=gnome ' . "$ENV{EXE_PAUP}";
	#print "# executing: cat $file_specific_modelblock | $ENV{EXE_PAUP} -n $infile > $tmpfile\n";
	system ("cat $file_specific_modelblock | $ENV{EXE_PAUP} -n $infile > $tmpfile"); 
    }
    else
    {
          system("cat $file_specific_modelblock | $ENV{'EXE_PAUP'} -n $infile > $tmpfile");
    }



    #unless(-s $model_scores_file){ sleep(1) }
    system("$ENV{'EXE_MODELTEST'} < $model_scores_file > $outfile");
    open(OUT,$outfile) || die "# run_MODELTEST : cannot read $outfile\n";
	 while(<OUT>)
	 {
	 	
		if(/^Likelihood settings from best-fit model \((.+)\) selected by AIC/)
		{
			$bestAICmodel = $1;
			$bestAICmodel =~ s/\+//g;
			#print "bestAICmodel $bestAICmodel\n";
		}
	 }
	 close(OUT);
	     
    # clean file
    #if($clean==1){ system(" rm -f $tmpfile $file_specific_model_log $model_scores_file $file_specific_modelblock"); }
    # added for use in Fasta2primers : $outfile
    if($clean==1){ system(" rm -f $tmpfile $file_specific_model_log $model_scores_file $file_specific_modelblock"); }
    
    return ($bestAICmodel,$outfile); 
}


sub convert_AAaln2NTaln
{
    my ($DNAfile,$AA_aln_file) = @_;

	 my @split = split(/\//,$DNAfile);
	 my $basefile = $split[$#split];

    my $outfile = (split(/\./,$basefile))[0] . "_aln.fna"; 
	 my ($seq,$sequence,$aa,$n_of_aa,$aligned_nts);
	 my ($align_positions,$FASTA,@length,$l);
	 
    # 1) read AA aligned fasta sequences
	 my %AAalign = read_FASTA_sequence($AA_aln_file);

    # 2) read NT fasta sequence
	 my %NTsequences = read_FASTA_sequence($DNAfile);
	 
    # 3) replicate alignment	 
	 foreach $seq (sort {$a<=>$b} (keys(%NTsequences)))
	 {
		#print "mira $seq $NTsequences{$seq}{'NAME'}\n";print "$NTsequences{$seq}{'SEQ'}\n$AAalign{$seq}{'SEQ'}\n";

		  $FASTA .= $NTsequences{$seq}{'NAME'};
		  
		  # 3.1) dissect codons
		  my @codons;
		  $sequence = $NTsequences{$seq}{'SEQ'};
		#print STDERR "#$sequence#\n";  
		while($sequence)
		  {
		   	push(@codons,substr($sequence,0,3));
			if(length($sequence) > 2){ $sequence = substr($sequence,3); }
			else{ die "# convert_AAaln2NTaln: sequence $NTsequences{$seq}{'NAME'} not divisible by 3\n"; }
		  }
		  
		  # 3.2) loop through aminoacids in $AAalign{$seq}
		  $sequence = $AAalign{$seq}{'SEQ'};
		  $n_of_aa = length($sequence);
		  $align_positions=0;
	 	  $aligned_nts = '';
		  for($aa=0;$aa<$n_of_aa;$aa++)
		  {
		   	if(substr($sequence,$aa,1) ne '-')
				{
					 $aligned_nts .= $codons[$align_positions];
					 $align_positions++;
				}
				else
				{
					 $aligned_nts .= '---';
				}
		  } 	
	 
	 	  $l = length($aligned_nts);
	 	  if(!grep(/$l/,@length)){ push(@length,$l);  }
		  
		  $FASTA .= $aligned_nts."\n";
	}
	 
	 if(scalar(@length) > 1)
	 {
	 	  print "# convert_AAaln2NTaln : warning, sequence lengths differ\n";
	 	  return $phyTools::ERROR;
	 }
	 
    # write the output FASTA file
	 open(FASTA,">$outfile") || die "# convert_AAaln2NTaln : cannot write to $outfile\n";
	 print FASTA $FASTA;
	 close FASTA;
    
	 return $outfile;
}

sub read_FASTA_sequence
{
	# if $remove_gap_col assumes FASTA file is actually a multiple alignment
	# in FASTA format
	# if $skipbadCDSs ignores nt sequences with length%3 > 0

    	my ( $infile, $remove_gap_cols, $skipbadCDSs, $skipidentical ) = @_;
	 
	 my (%FASTA,$name,$seq,$n_of_sequences,$length,$maxlength,$pos,$isgap,$seqid);
	 
	 $n_of_sequences = $maxlength = 0;
	 open(FASTA,$infile) || die "# read_FASTA_sequence: cannot read $infile $!:\n";
	 while(<FASTA>)
	 {
	 	  if(/^\>/)
		  {
		   	$name = $_; 
			$n_of_sequences++;
			$seqid = sprintf("%0$ENV{'TRAILINGZEROS'}d",$n_of_sequences); #print "#$seqid#\n";
			$FASTA{$seqid}{'NAME'} = $name;
		  }    	 	   	
		  else
		  {
		   	$FASTA{$seqid}{'SEQ'} .= $_;
			$FASTA{$seqid}{'SEQ'} =~ s/[\s|\n]//g;
				
			$length = length($FASTA{$seqid}{'SEQ'});
			if($length > $maxlength){ $maxlength = $length; }
		  }
	 }
	 close(FASTA);

	if($skipbadCDSs)
	{
		# check length is divisible by 3 and that no inframe STOP codons exist
		my %goodCDSs_FASTA;
		foreach $seq (keys(%FASTA))
                {
			my $ntseq = Bio::Seq->new( -display_id => 'tmp', -seq => $FASTA{$seq}{'SEQ'} );
                        my $protseq = $ntseq->translate()->seq(); chop $protseq;
			if($protseq =~ /\*/)
                        {
                                print "# read_FASTA_sequence : skipped CDS sequence (inframe STOP codon) $FASTA{$seq}{'NAME'}\n";
                                next;
                        }
		
                	if($FASTA{$seq}{'SEQ'} !~ /QWERYIPSDFHKLMNV/ && length($FASTA{$seq}{'SEQ'})%3)
			{
				my $choppednts = 0;
				while(length($FASTA{$seq}{'SEQ'})%3){ chop $FASTA{$seq}{'SEQ'}; $choppednts++; }
				print "# read_FASTA_sequence : chopped CDS sequence (3' $choppednts nts) $FASTA{$seq}{'NAME'}\n";
			}
			
			#if(length($protseq) < 100){ print "mira: $protseq\n"; }
			$goodCDSs_FASTA{$seq} = $FASTA{$seq};
		}
		%FASTA = %goodCDSs_FASTA;
	}

	if($skipidentical)
	{
		my (%nrFASTA,%identical,$seq2);
                foreach $seq (keys(%FASTA))
                {
			$FASTA{$seq}{'IDENTICALS'} = '';
			next if($identical{$seq});
			foreach $seq2 (keys(%FASTA))
			{
				next if($seq == $seq2);
				if($FASTA{$seq}{'SEQ'} eq $FASTA{$seq2}{'SEQ'} ||        # same length
					$FASTA{$seq}{'SEQ'} =~ /$FASTA{$seq2}{'SEQ'}/ || # different length 
					$FASTA{$seq2}{'SEQ'} =~ /$FASTA{$seq}{'SEQ'}/)
				{
					$identical{$seq2} = $seq;
					$FASTA{$seq}{'IDENTICALS'} .= "$seq2,"; #print "mira $seq $seq2\n";
					$FASTA{$seq}{'NUMBER_IDENTICALS'}++;
				}
			}
		}

		foreach $seq (keys(%FASTA))
                {
			if($identical{$seq})
			{
				print "# read_FASTA_sequence : skipped sequence identical to $identical{$seq}: $FASTA{$seq}{'NAME'}\n";
				next;
			}
			$nrFASTA{$seq} = $FASTA{$seq};

			# keep track of identical sequences
			if($FASTA{$seq}{'NUMBER_IDENTICALS'})
			{ 
				chomp($nrFASTA{$seq}{'NAME'});
				$nrFASTA{$seq}{'NAME'} .= " | identical sequences=$FASTA{$seq}{'NUMBER_IDENTICALS'}\n";
			}
		}
		%FASTA = %nrFASTA;
	}

	#if($remove_gap_cols && $remove_gap_cols == 1) 
	if($remove_gap_cols)
	{
		for($pos=0;$pos < $length; $pos++)
		{	
			$isgap=0;
			foreach $seq (keys(%FASTA))
			{
				if(substr($FASTA{$seq}{'SEQ'},$pos,1))
				{
					if(substr($FASTA{$seq}{'SEQ'},$pos,1) eq '-'){ $isgap=1; last; }
				}
				else{ $isgap=1; last }
			}
			
			if($isgap == 0)
			{
				foreach $seq (keys(%FASTA))
				{
					$FASTA{$seq}{'NOGAPSEQ'} .= substr($FASTA{$seq}{'SEQ'},$pos,1);
				}
			}
		}
		
		
		foreach $seq (keys(%FASTA))
		{
			$FASTA{$seq}{'SEQ'} = $FASTA{$seq}{'NOGAPSEQ'}
		}
	}

	return %FASTA;
}

# added in 2011 as arrays take less RAM and are faster to browse:
# http://bioinfoperl.blogspot.com/2011/02/estructuras-de-datos-basicas-en-perl.html
sub read_FASTA_file_array
{
	# in FASTA format
	# returns a reference to a 2D array for 4 secondary keys: NAME,SEQ,IDENTICALS,NUMBER_IDENTICALS
	# first valid index (first sequence) is '0'

    	my ( $infile, $skipamino, $skipidentical ) = @_;
	 
	my (@FASTA,$name,$seq,$n_of_sequences);
	 
	$n_of_sequences = -1; 
	open(FASTA,"<$infile") || die "# read_FASTA_sequence: cannot read $infile $!:\n";
	while(<FASTA>)
	{
            next if(/^$/);
            next if(/^#/);
	    if(/^\>(.*?)[\n\r]/)
            {
		$n_of_sequences++; # first sequence ID is 0
	   	$name = $1; #print "# $n_of_sequences $name\n"; 			
		$FASTA[$n_of_sequences][NAME] = $name;
            }    	 	   	
            elsif($n_of_sequences>-1)
            {
		$_ =~ s/[\s|\n|\-|\.]//g;
		$FASTA[$n_of_sequences][SEQ] .= $_; # conserve case
            }
	}
	close(FASTA);

	if($n_of_sequences == -1){ return \@FASTA }
		
	if($skipamino || $skipidentical)
	{
		my ($n_of_nr_sequences,@nrFASTA,%identical) = 0;
		
		if($skipidentical)
		{
			foreach $seq ( 0 .. $n_of_sequences-1 )
			{
				$FASTA[$seq][IDENTICALS] = ''; 
				next if($identical{$seq});
				foreach my $seq2 ( $seq + 1 .. $n_of_sequences )
				{
					if($FASTA[$seq][SEQ] eq $FASTA[$seq2][SEQ] ||        # same length
						$FASTA[$seq][SEQ] =~ /$FASTA[$seq2][SEQ]/ || # different length 
						$FASTA[$seq2][SEQ] =~ /$FASTA[$seq][SEQ]/)
					{
						$identical{$seq2} = $seq;
						$FASTA[$seq][IDENTICALS] .= "$seq2,"; #print "mira $seq $FASTA[$seq][IDENTICALS]\n";
						$FASTA[$seq][NUMBER_IDENTICALS]++;
					}
				}
			}
		}	

		foreach $seq ( 0 .. $n_of_sequences )
                {
			# skip protein sequences
			if($FASTA[$seq][SEQ] =~ /QWERTYIPSDFHKLVNM/i)
			{
				print "# read_FASTA_sequence : skipped amino acid sequence ($FASTA[$seq][NAME])\n";
				next;
			}
			
			if($identical{$seq})
			{
				print "# read_FASTA_sequence : skipped identical sequence: ".
					"$FASTA[$seq][NAME]\n";
				next;
			}
			
			$nrFASTA[$n_of_nr_sequences] = $FASTA[$seq]; 
			
			# keep track of identical sequences
			if($FASTA[$seq][NUMBER_IDENTICALS])
			{ 
				chomp($nrFASTA[$n_of_nr_sequences][NAME]);
				$nrFASTA[$n_of_nr_sequences][NAME] .=
					" | identical sequences=$nrFASTA[$seq][NUMBER_IDENTICALS]";
			}	
			$n_of_nr_sequences++;
		}
		
		return \@nrFASTA;
	}
	else{ return \@FASTA }
}

sub check_aligned_FASTA_sequences
{
	# receives a FASTA hash generated by read_FASTA_sequences
	# returns 1 if sequences are aligned
	
	my (%sequences) = @_;
	my $length = 0;

	foreach my $seq (sort {$a<=>$b} (keys(%sequences)))
	{
         	if(!$length){ $length = length($sequences{$seq}{'SEQ'}); }
		else{ if(length($sequences{$seq}{'SEQ'}) != $length){ return 0 } }
	}

	return 1;
}

sub getPHYMLparameters
{
	my ($model) = @_;
	
	# borrowed from MrAIC
	my (%MODELnr,%MODELcmd,$m);
	
	$MODELnr{0}	=	"JC69";
	$MODELnr{1}	=	"JC69I";
	$MODELnr{2}	=	"JC69G";
	$MODELnr{3}	=	"JC69IG";
	$MODELnr{4}	=	"F81";
	$MODELnr{5}	=	"F81I";
	$MODELnr{6}	=	"F81G";
	$MODELnr{7}	=	"F81IG";
	$MODELnr{8}	=	"K80";
	$MODELnr{9}	=	"K80I";
	$MODELnr{10}	=	"K80G";
	$MODELnr{11}	=	"K80IG";
	$MODELnr{12}	=	"HKY";
	$MODELnr{13}	=	"HKYI";
	$MODELnr{14}	=	"HKYG";
	$MODELnr{15}	=	"HKYIG";
	$MODELnr{16}	=	"TrNef";
	$MODELnr{17}	=	"TrNefI";
	$MODELnr{18}	=	"TrNefG";
	$MODELnr{19}	=	"TrNefIG";
	$MODELnr{20}	=	"TrN";	# aka TN93
	$MODELnr{21}	=	"TrNI";
	$MODELnr{22}	=	"TrNG";
	$MODELnr{23}	=	"TrNIG";
	
	$MODELnr{24}	=	"K3P|K81";	# aka K81
	$MODELnr{25}	=	"K3PI|K81I";
	$MODELnr{26}	=	"K3PG|K81G";
	$MODELnr{27}	=	"K3PIG|K81IG";
	$MODELnr{28}	=	"K3Puf|K81uf";
	$MODELnr{29}	=	"K3PufI|K81ufI";
	$MODELnr{30}	=	"K3PufG|K81ufG";
	$MODELnr{31}	=	"K3PufIG|K81ugIG";
	
	$MODELnr{32}	=	"TIMef";
	$MODELnr{33}	=	"TIMefI";
	$MODELnr{34}	=	"TIMefG";
	$MODELnr{35}	=	"TIMefIG";
	$MODELnr{36}	=	"TIM";
	$MODELnr{37}	=	"TIMI";
	$MODELnr{38}	=	"TIMG";
	$MODELnr{39}	=	"TIMIG";
	$MODELnr{40}	=	"TVMef";
	$MODELnr{41}	=	"TVMefI";
	$MODELnr{42}	=	"TVMefG";
	$MODELnr{43}	=	"TVMefIG";
	$MODELnr{44}	=	"TVM";
	$MODELnr{45}	=	"TVMI";
	$MODELnr{46}	=	"TVMG";
	$MODELnr{47}	=	"TVMIG";
	$MODELnr{48}	=	"SYM";
	$MODELnr{49}	=	"SYMI";
	$MODELnr{50}	=	"SYMG";
	$MODELnr{51}	=	"SYMIG";
	$MODELnr{52}	=	"GTR";
	$MODELnr{53}	=	"GTRI";
	$MODELnr{54}	=	"GTRG";
	$MODELnr{55}	=	"GTRIG";

	$MODELcmd{0}	=	"M\nM\nM\nM\nM\nY\n";	# JC69
	$MODELcmd{1}	=	"M\nM\nM\nM\nM\nV\nY\nY\n";	# JC69I
	$MODELcmd{2}	=	"M\nM\nM\nM\nM\nR\nA\nY\nY\n";	# JC69G
	$MODELcmd{3}	=	"M\nM\nM\nM\nM\nR\nA\nY\nV\nY\nY\n";	# JC69IG
	$MODELcmd{4}	=	"M\nM\nM\nM\nM\nM\nM\nY\n";	# F81
	$MODELcmd{5}	=	"M\nM\nM\nM\nM\nM\nM\nV\nY\nY\n";	# F81I
	$MODELcmd{6}	=	"M\nM\nM\nM\nM\nM\nM\nR\nA\nY\nY\n";	# F81G
	$MODELcmd{7}	=	"M\nM\nM\nM\nM\nM\nM\nR\nA\nY\nV\nY\nY\n";	# F81IG
	$MODELcmd{8}	=	"M\nM\nM\nM\nM\nM\nT\nY\nY\n";	# K80
	$MODELcmd{9}	=	"M\nM\nM\nM\nM\nM\nT\nY\nV\nY\nY\n";	# K80I
	$MODELcmd{10}	=	"M\nM\nM\nM\nM\nM\nT\nY\nR\nA\nY\nY\n";	# K80G
	$MODELcmd{11}	=	"M\nM\nM\nM\nM\nM\nT\nY\nR\nA\nY\nV\nY\nY\n";	# K80IG
	$MODELcmd{12}	=	"T\nY\nE\nY\n";	# HKY
	$MODELcmd{13}	=	"T\nY\nE\nV\nY\nY\n";	# HKYI
	$MODELcmd{14}	=	"T\nY\nE\nR\nA\nY\nY\n";	# HKYG
	$MODELcmd{15}	=	"T\nY\nE\nR\nA\nY\nV\nY\nY\n";	# HKYIG
	$MODELcmd{16}	=	"M\nM\nM\nM\nF\n0.25\n0.25\n0.25\n0.25\nK\n010020\n1.00\n1.00\nW\nY\n";	# TrNef
	$MODELcmd{17}	=	"M\nM\nM\nM\nF\n0.25\n0.25\n0.25\n0.25\nK\n010020\n1.00\n1.00\nW\nV\nY\nY\n";	# TrNefI
	$MODELcmd{18}	=	"M\nM\nM\nM\nF\n0.25\n0.25\n0.25\n0.25\nK\n010020\n1.00\n1.00\nW\nR\nA\nY\nY\n";	# TrNefG
	$MODELcmd{19}	=	"M\nM\nM\nM\nF\n0.25\n0.25\n0.25\n0.25\nK\n010020\n1.00\n1.00\nW\nR\nA\nY\nV\nY\nY\n";	# TrNefIG
	$MODELcmd{20}	=	"M\nM\nE\nT\nY\nY\n";	# TrN
	$MODELcmd{21}	=	"M\nM\nE\nT\nY\nV\nY\nY\n";	# TrNI
	$MODELcmd{22}	=	"M\nM\nE\nT\nY\nR\nA\nY\nY\n";	# TrNG
	$MODELcmd{23}	=	"M\nM\nE\nT\nY\nR\nA\nY\nV\nY\nY\n";	# TrNIG
	$MODELcmd{24}	=	"M\nM\nM\nM\nF\n0.25\n0.25\n0.25\n0.25\nK\n012210\n1.00\n1.00\nW\nY\n";	# K3P
	$MODELcmd{25}	=	"M\nM\nM\nM\nF\n0.25\n0.25\n0.25\n0.25\nK\n012210\n1.00\n1.00\nW\nV\nY\nY\n";	# K3PI
	$MODELcmd{26}	=	"M\nM\nM\nM\nF\n0.25\n0.25\n0.25\n0.25\nK\n012210\n1.00\n1.00\nW\nR\nA\nY\nY\n";	# K3PG
	$MODELcmd{27}	=	"M\nM\nM\nM\nF\n0.25\n0.25\n0.25\n0.25\nK\n012210\n1.00\n1.00\nW\nR\nA\nY\nV\nY\nY\n";	# K3PIG
	$MODELcmd{28}	=	"M\nM\nM\nM\nK\n012210\n1.00\n1.00\nW\nY\n";	# K3Puf
	$MODELcmd{29}	=	"M\nM\nM\nM\nK\n012210\n1.00\n1.00\nW\nV\nY\nY\n";	# K3PufI
	$MODELcmd{30}	=	"M\nM\nM\nM\nK\n012210\n1.00\n1.00\nW\nR\nA\nY\nY\n";	# K3PufG
	$MODELcmd{31}	=	"M\nM\nM\nM\nK\n012210\n1.00\n1.00\nW\nR\nA\nY\nV\nY\nY\n";	# K3PufIG
	$MODELcmd{32}	=	"M\nM\nM\nM\nF\n0.25\n0.25\n0.25\n0.25\nK\n012230\n1.00\n1.00\n1.00\nW\nY\n";	# TIMef
	$MODELcmd{33}	=	"M\nM\nM\nM\nF\n0.25\n0.25\n0.25\n0.25\nK\n012230\n1.00\n1.00\n1.00\nW\nV\nY\nY\n";	# TIMefI
	$MODELcmd{34}	=	"M\nM\nM\nM\nF\n0.25\n0.25\n0.25\n0.25\nK\n012230\n1.00\n1.00\n1.00\nW\nR\nA\nY\nY\n";	# TIMefG
	$MODELcmd{35}	=	"M\nM\nM\nM\nF\n0.25\n0.25\n0.25\n0.25\nK\n012230\n1.00\n1.00\n1.00\nW\nR\nA\nY\nV\nY\nY\n";	# TIMefIG
	$MODELcmd{36}	=	"M\nM\nM\nM\nK\n012230\n1.00\n1.00\n1.00\nW\nY\n";	# TIM
	$MODELcmd{37}	=	"M\nM\nM\nM\nK\n012230\n1.00\n1.00\n1.00\nW\nV\nY\nY\n";	# TIMI
	$MODELcmd{38}	=	"M\nM\nM\nM\nK\n012230\n1.00\n1.00\n1.00\nW\nR\nA\nY\nY\n";	# TIMG
	$MODELcmd{39}	=	"M\nM\nM\nM\nK\n012230\n1.00\n1.00\n1.00\nW\nR\nA\nY\nV\nY\nY\n";	# TIMIG
	$MODELcmd{40}	=	"M\nM\nM\nM\nF\n0.25\n0.25\n0.25\n0.25\nK\n012314\n1.00\n1.00\n1.00\n1.00\nW\nY\n";	# TVMef
	$MODELcmd{41}	=	"M\nM\nM\nM\nF\n0.25\n0.25\n0.25\n0.25\nK\n012314\n1.00\n1.00\n1.00\n1.00\nW\nV\nY\nY\n";	# TVMefI
	$MODELcmd{42}	=	"M\nM\nM\nM\nF\n0.25\n0.25\n0.25\n0.25\nK\n012314\n1.00\n1.00\n1.00\n1.00\nW\nR\nA\nY\nY\n";	# TVMefG
	$MODELcmd{43}	=	"M\nM\nM\nM\nF\n0.25\n0.25\n0.25\n0.25\nK\n012314\n1.00\n1.00\n1.00\n1.00\nW\nR\nA\nY\nV\nY\nY\n";	# TVMefIG
	$MODELcmd{44}	=	"M\nM\nM\nM\nK\n012314\n1.00\n1.00\n1.00\n1.00\nW\nY\n";	# TVM
	$MODELcmd{45}	=	"M\nM\nM\nM\nK\n012314\n1.00\n1.00\n1.00\n1.00\nW\nV\nY\nY\n";	# TVMI
	$MODELcmd{46}	=	"M\nM\nM\nM\nK\n012314\n1.00\n1.00\n1.00\n1.00\nW\nR\nA\nY\nY\n";	# TVMG
	$MODELcmd{47}	=	"M\nM\nM\nM\nK\n012314\n1.00\n1.00\n1.00\n1.00\nW\nR\nA\nY\nV\nY\nY\n";	# TVMIG
	$MODELcmd{48}	=	"M\nM\nM\nM\nF\n0.25\n0.25\n0.25\n0.25\nK\n012345\n1.00\n1.00\n1.00\n1.00\n1.00\nW\nY\n";	# SYM
	$MODELcmd{49}	=	"M\nM\nM\nM\nF\n0.25\n0.25\n0.25\n0.25\nK\n012345\n1.00\n1.00\n1.00\n1.00\n1.00\nW\nV\nY\nY\n";	# SYMI
	$MODELcmd{50}	=	"M\nM\nM\nM\nF\n0.25\n0.25\n0.25\n0.25\nK\n012345\n1.00\n1.00\n1.00\n1.00\n1.00\nW\nR\nA\nY\nY\n";	# SYMG
	$MODELcmd{51}	=	"M\nM\nM\nM\nF\n0.25\n0.25\n0.25\n0.25\nK\n012345\n1.00\n1.00\n1.00\n1.00\n1.00\nW\nR\nA\nY\nV\nY\nY\n";	# SYMIG
	$MODELcmd{52}	=	"M\nM\nM\nE\nY\n";	# GTR
	$MODELcmd{53}	=	"M\nM\nM\nE\nV\nY\nY\n";	# GTRI
	$MODELcmd{54}	=	"M\nM\nM\nE\nR\nA\nY\nY\n";	# GTRG
	$MODELcmd{55}	=	"M\nM\nM\nE\nR\nA\nY\nV\nY\nY\n";	# GTRIG

	foreach $m (sort {$a<=>$b} (keys(%MODELnr)))
	{
		#print "#$MODELnr{$m}#$model#\n";
		if($MODELnr{$m} =~ /$model/)
		{
			return $MODELcmd{$m};
		}
	}

	return $phyTools::ERROR;
}


sub getPHYML3parameters
{
	my ($model) = @_;
	
	my (%MODELnr,%MODELcmd,$m);
	
	$MODELnr{0}	=	"JC69";
	$MODELnr{1}	=	"JC69I";
	$MODELnr{2}	=	"JC69G";
	$MODELnr{3}	=	"JC69IG";
	$MODELnr{4}	=	"F81";
	$MODELnr{5}	=	"F81I";
	$MODELnr{6}	=	"F81G";
	$MODELnr{7}	=	"F81IG";
	$MODELnr{8}	=	"K80";
	$MODELnr{9}	=	"K80I";
	$MODELnr{10}	=	"K80G";
	$MODELnr{11}	=	"K80IG";
	$MODELnr{12}	=	"HKY";
	$MODELnr{13}	=	"HKYI";
	$MODELnr{14}	=	"HKYG";
	$MODELnr{15}	=	"HKYIG";
	$MODELnr{16}	=	"TrNef";
	$MODELnr{17}	=	"TrNefI";
	$MODELnr{18}	=	"TrNefG";
	$MODELnr{19}	=	"TrNefIG";
	$MODELnr{20}	=	"TrN";	# aka TN93
	$MODELnr{21}	=	"TrNI";
	$MODELnr{22}	=	"TrNG";
	$MODELnr{23}	=	"TrNIG";
	$MODELnr{24}	=	"K3P";	# aka K81
	$MODELnr{25}	=	"K3PI";
	$MODELnr{26}	=	"K3PG";
	$MODELnr{27}	=	"K3PIG";
	$MODELnr{28}	=	"K3Puf";
	$MODELnr{29}	=	"K3PufI";
	$MODELnr{30}	=	"K3PufG";
	$MODELnr{31}	=	"K3PufIG";
	$MODELnr{32}	=	"TIMef";
	$MODELnr{33}	=	"TIMefI";
	$MODELnr{34}	=	"TIMefG";
	$MODELnr{35}	=	"TIMefIG";
	$MODELnr{36}	=	"TIM";
	$MODELnr{37}	=	"TIMI";
	$MODELnr{38}	=	"TIMG";
	$MODELnr{39}	=	"TIMIG";
	$MODELnr{40}	=	"TVMef";
	$MODELnr{41}	=	"TVMefI";
	$MODELnr{42}	=	"TVMefG";
	$MODELnr{43}	=	"TVMefIG";
	$MODELnr{44}	=	"TVM";
	$MODELnr{45}	=	"TVMI";
	$MODELnr{46}	=	"TVMG";
	$MODELnr{47}	=	"TVMIG";
	$MODELnr{48}	=	"SYM";
	$MODELnr{49}	=	"SYMI";
	$MODELnr{50}	=	"SYMG";
	$MODELnr{51}	=	"SYMIG";
	$MODELnr{52}	=	"GTR";
	$MODELnr{53}	=	"GTRI";
	$MODELnr{54}	=	"GTRG";
	$MODELnr{55}	=	"GTRIG";
	
	#May2015 Bruno
	if($model =~ /K81/){ $model =~ s/K81/K3P/ } # https://groups.google.com/forum/#!topic/beast-users/6SwVGjJjCYs

	$MODELcmd{0}	=	"-m JC69 -c 1 -v 0";	        # JC69
	$MODELcmd{1}	=	"-m JC69 -c 1 -v e";	        # JC69I
	$MODELcmd{2}	=	"-m JC69 -c 4 -a e -v 0";	# JC69G
	$MODELcmd{3}	=	"-m JC69 -c 4 -a e -v e";	# JC69IG
	$MODELcmd{4}	=       "-m F81 -c 1 -v 0";             #  F81
	$MODELcmd{5}	=       "-m F81 -c 1 -v e";             # F81I
	$MODELcmd{6}	=       "-m F81 -c 4 -a e -v 0";        # F81G
	$MODELcmd{7}	=       "-m F81 -c 4 -a e -v e";        # F81IG
	$MODELcmd{8}	=	"-m K80 -c 1 -v 0";	        # K2P
	$MODELcmd{9}	=	"-m K80 -c 1 -v e";	        # K2PI
	$MODELcmd{10}	=	"-m K80 -c 4 -a e -v 0";        # K2PG
	$MODELcmd{11}	=	"-m K80 -c 4 -a e -v e";        # K2PIG
	$MODELcmd{12}	=	"-m HKY85 -c 1 -v 0";		# HKY
	$MODELcmd{13}	=	"-m HKY85 -c 1 -v e";		# HKYI
	$MODELcmd{14}	=	"-m HKY85 -c 4 -a e -v 0";	# HKYG
	$MODELcmd{15}	=	"-m HKY85 -c 4 -a e -v e";	# HKYIG
	$MODELcmd{16}	=	"-f -fA 0.25 -fC 0.25 -fG 0.25 -fT 0.25 -m 010020 -c 1 -v 0";	   # TrNef
	$MODELcmd{17}	=	"-f -fA 0.25 -fC 0.25 -fG 0.25 -fT 0.25 -m 010020 -c 1 -v e";	   # TrNefI
	$MODELcmd{18}	=	"-f -fA 0.25 -fC 0.25 -fG 0.25 -fT 0.25 -m 010020 -c 1 -a e -v 0"; # TrNefG
	$MODELcmd{19}	=	"-f -fA 0.25 -fC 0.25 -fG 0.25 -fT 0.25 -m 010020 -c 1 -a e -v e"; # TrNefIG
	$MODELcmd{20}	=	"-m TN93 -c 1 -v 0";	       # TrN
	$MODELcmd{21}	=	"-m TN93 -c 1 -v e";	       # TrNI
	$MODELcmd{22}	=	"-m TN93 -c 4 -a e -v 0";      # TrNG
	$MODELcmd{23}	=	"-m TN93 -c 4 -a e -v e";      # TrNIG
	$MODELcmd{24}	=	"-f -fA 0.25 -fC 0.25 -fG 0.25 -fT 0.25 -m 012210 -c 1 -v 0";	   # K3P
	$MODELcmd{25}	=	"-f -fA 0.25 -fC 0.25 -fG 0.25 -fT 0.25 -m 012210 -c 1 -v e";	   # K3PI
	$MODELcmd{26}	=	"-f -fA 0.25 -fC 0.25 -fG 0.25 -fT 0.25 -m 012210 -c 1 -a e -v 0"; # K3PG
	$MODELcmd{27}	=	"-f -fA 0.25 -fC 0.25 -fG 0.25 -fT 0.25 -m 012210 -c 1 -a e -v e"; # K3PIG
	#May2015 Bruno due to phyml warining: "Please define the data type (nt or aa) before setting the -f option"
	#http://phyml.googlecode.com/svn-history/r671/trunk/src/cl.c
	$MODELcmd{28}   =       "-d nt -f e -m 012210 -c 1 -v 0";          # K3Puf
        $MODELcmd{29}   =       "-d nt -f e -m 012210 -c 1 -v e";          # K3PufI
        $MODELcmd{30}   =       "-d nt -f e -m 012210 -c 1 -a e -v 0";     # K3PufG
        $MODELcmd{31}   =       "-d nt -f e -m 012210 -c 1 -a e -v e";     # K3PufIG
	#$MODELcmd{28}	=	"-f e -m 012210 -c 1 -v 0";          # K3Puf
	#$MODELcmd{29}	=	"-f e -m 012210 -c 1 -v e";	     # K3PufI
	#$MODELcmd{30}	=	"-f e -m 012210 -c 1 -a e -v 0";     # K3PufG
	#$MODELcmd{31}	=	"-f e -m 012210 -c 1 -a e -v e";     # K3PufIG
	$MODELcmd{32}	=	"-f -fA 0.25 -fC 0.25 -fG 0.25 -fT 0.25 -m 012230 -c 1 -v 0";	    # TIMef
	$MODELcmd{33}	=	"-f -fA 0.25 -fC 0.25 -fG 0.25 -fT 0.25 -m 012230 -c 1 -v e";	    # TIMefI
	$MODELcmd{34}	=	"-f -fA 0.25 -fC 0.25 -fG 0.25 -fT 0.25 -m 012230 -c 1 -a e -v 0";  # TIMefG
	$MODELcmd{35}	=	"-f -fA 0.25 -fC 0.25 -fG 0.25 -fT 0.25 -m 012230 -c 1 -a e -v e";  # TIMefIG
	$MODELcmd{36}   =       "-d nt -f e -m 012230 -c 1 -v 0";          # TIM
        $MODELcmd{37}   =       "-d nt -f e -m 012230 -c 1 -v e";          # TIMI
        $MODELcmd{38}   =       "-d nt -f e -m 012230 -c 1 -a e -v 0";     # TIMG
        $MODELcmd{39}   =       "-d nt -f e -m 012230 -c 1 -a e -v e";     # TIMIG
	#$MODELcmd{36}	=	"-f e -m 012230 -c 1 -v 0";	     # TIM
	#$MODELcmd{37}	=	"-f e -m 012230 -c 1 -v e";	     # TIMI
	#$MODELcmd{38}	=	"-f e -m 012230 -c 1 -a e -v 0";     # TIMG
	#$MODELcmd{39}	=	"-f e -m 012230 -c 1 -a e -v e";     # TIMIG
	$MODELcmd{40}	=	"-f -fA 0.25 -fC 0.25 -fG 0.25 -fT 0.25 -m 012314 -c 1 -v 0";	    # TVMef
	$MODELcmd{41}	=	"-f -fA 0.25 -fC 0.25 -fG 0.25 -fT 0.25 -m 012314 -c 1 -v e";	    # TVMefI
	$MODELcmd{42}	=	"-f -fA 0.25 -fC 0.25 -fG 0.25 -fT 0.25 -m 012314 -c 1 -a e -v 0";  # TVMefG
	$MODELcmd{43}	=	"-f -fA 0.25 -fC 0.25 -fG 0.25 -fT 0.25 -m 012314 -c 1 -a e -v e";  # TVMefIG
	$MODELcmd{44}   =       "-d nt -f e -m 012314 -c 1 -v 0";          # TVM
        $MODELcmd{45}   =       "-d nt -f e -m 012314 -c 1 -v e";          # TVMI
        $MODELcmd{46}   =       "-d nt -f e -m 012314 -c 1 -a e -v 0";     # TVMG
        $MODELcmd{47}   =       "-d nt -f e -m 012314 -c 1 -a e -v e";     # TVMIG	
	#$MODELcmd{44}	=	"-f e -m 012314 -c 1 -v 0";	     # TVM
	#$MODELcmd{45}	=	"-f e -m 012314 -c 1 -v e";	     # TVMI
	#$MODELcmd{46}	=	"-f e -m 012314 -c 1 -a e -v 0";     # TVMG
	#$MODELcmd{47}	=	"-f e -m 012314 -c 1 -a e -v e";     # TVMIG
	$MODELcmd{48}	=	"-f -fA 0.25 -fC 0.25 -fG 0.25 -fT 0.25 -m 012345 -c 1 -v 0";	     # SYM
	$MODELcmd{49}	=	"-f -fA 0.25 -fC 0.25 -fG 0.25 -fT 0.25 -m 012345 -c 1 -v e";	     # SYMI
	$MODELcmd{50}	=	"-f -fA 0.25 -fC 0.25 -fG 0.25 -fT 0.25 -m 012345 -c 1 -a e -v 0";   # SYMG
	$MODELcmd{51}	=	"-f -fA 0.25 -fC 0.25 -fG 0.25 -fT 0.25 -m 012345 -c 1 -a e -v e";   # SYMIG
	$MODELcmd{52}	=	"-m GTR -c 1 -v 0";	   # GTR
	$MODELcmd{53}	=	"-m GTR -c 1 -v e";	   # GTRI
	$MODELcmd{54}	=	"-m GTR -c 4 -a e -v 0";   # GTRG
	$MODELcmd{55}	=	"-m GTR -c 4 -a e -v e";   # GTRIG

	foreach $m (sort {$a<=>$b} (keys(%MODELnr)))
	{
		#print "#$MODELnr{$m}#$model#\n";
		if($MODELnr{$m} =~ /$model/)
		{
			return $MODELcmd{$m};
		}
	}
	return $phyTools::ERROR;
}


sub run_PHYML
{
	my ($DNA_phylip_file,$PHYMLparams,$clean) = @_;

	my $treefile = $DNA_phylip_file . "_phyml_tree.txt";
	my $lkfile =   $DNA_phylip_file . "_phyml_lk.txt";
	my $statfile = $DNA_phylip_file . "_phyml_stat.txt";
	my $cmdfile = (split(/\./,$DNA_phylip_file))[0] . "_tmp_cmd_phyml";
   	my $tmpfile = (split(/\./,$DNA_phylip_file))[0] . "_tmp_log_phyml";
	
	my ($newick,$aLRT,@aLRT_values);
	my ($mean_aLRT,$median_aLRT,$n_of_splits)=(0.0,0.0,0);
	
	# delete preexisting files if necessary, otherwise PHYML la caga
	if(-e $treefile){ system("rm -f $treefile") }
	if(-e $lkfile){ system("rm -f $lkfile") }
	if(-e $statfile){ system("rm -f $statfile") }
	
	# write command file
	open(CMD,">$cmdfile")|| die "# run_PHYML : cannot write $cmdfile\n";
	print CMD "$DNA_phylip_file\n$PHYMLparams\n";
	close (CMD);
	
   	# run PHYML
   	system("$ENV{'EXE_PHYML'} < $cmdfile > $tmpfile");
   
	# parse tree
 	open(TREE,$treefile) || die "# rum_PHYML : cannot read $treefile\n";
	while(<TREE>){ $newick .= $_; }
	close(TREE);

	while($newick =~ /\)(\d+\.\d+):/g)
	{
		$aLRT = $1;
		
		push(@aLRT_values,$aLRT);
		$mean_aLRT += $aLRT;
		$n_of_splits++;
	} 
	
	$mean_aLRT = sprintf("%1.2f",$mean_aLRT/$n_of_splits);
	$median_aLRT = sprintf("%1.2f",calc_median(@aLRT_values));
   
   	# clean file
   	if($clean==1){ system(" rm -f $tmpfile $cmdfile"); }

   # check the output files are there
   if(!(-s $treefile) || !(-s $lkfile) || !(-s $lkfile)){ return ($phyTools::ERROR,$phyTools::ERROR,$phyTools::ERROR); }
   else{ return ($treefile,$lkfile,$statfile,$mean_aLRT,$median_aLRT); }
}


sub run_PHYML3
{
	my ($phylip_file,$PHYMLparams,$clean,$search,$optimizations,$user_tree_file,$random_trees,$no_of_random_trees) = @_;

	my $treefile = $phylip_file . "_phyml_tree.txt";
	my $statfile = $phylip_file . "_phyml_stats.txt";
	my $cmdfile = (split(/\./,$phylip_file))[0] . "_tmp_cmd_phyml";
   	my $tmpfile = (split(/\./,$phylip_file))[0] . "_tmp_log_phyml";
	
	my ($newick,$aLRT,@aLRT_values);
	my ($mean_aLRT,$median_aLRT,$n_of_splits)=(0.0,0.0,0);
	
	# delete preexisting files if necessary, otherwise PHYML la caga
	if(-e $treefile){ system("rm -f $treefile") }
	if(-e $statfile){ system("rm -f $statfile") }
	
	# write command file
	
	my $additional_opts_string = ''; 
	$additional_opts_string .= " -s $search -o $optimizations";
	if($user_tree_file){ $additional_opts_string .= " -u $user_tree_file"; }
	if($random_trees){ $additional_opts_string .= " --rand_start --n_rand_starts $no_of_random_trees"; }
	
	open(CMD,">$cmdfile")|| die "# run_PHYML3 : cannot write $cmdfile\n";
	print CMD "-i $phylip_file $PHYMLparams $additional_opts_string\n";
	close (CMD);
	
	my $cmd_string = "-i $phylip_file $PHYMLparams $additional_opts_string\n";
	print "# will execute $ENV{'EXE_PHYML3'} $cmd_string > $tmpfile";
	
   	# run PHYML3
   	system "$ENV{'EXE_PHYML3'} $cmd_string > $tmpfile";
   
	# parse tree
 	open(TREE,$treefile) || die "# rum_PHYML3 : cannot read $treefile\n";
	while(<TREE>){ $newick .= $_; }
	close(TREE);

	while($newick =~ /\)(\d+\.\d+):/g)
	{
		$aLRT = $1;
		
		push(@aLRT_values,$aLRT);
		$mean_aLRT += $aLRT;
		$n_of_splits++;
	} 
	
	$mean_aLRT = sprintf("%1.2f",$mean_aLRT/$n_of_splits);
	$median_aLRT = sprintf("%1.2f",calc_median(@aLRT_values));
   
   	# clean file
   	if($clean){ system(" rm -f $tmpfile $cmdfile"); }

   # check the output files are there
   if(!(-s $treefile) || !(-s $statfile) ){ return ($phyTools::ERROR,$phyTools::ERROR); }
   else{ return ($treefile,$statfile,$mean_aLRT,$median_aLRT); }
}


sub run_codeml
{
	# codeml needs a codon_aln_file,a tree_file and a codeml_ctrl_file to run; the former 2 are already specified in the ctrl file
	my($codeml_ctrl_file) = @_;
	
	# run codeml
	system("$ENV{'EXE_CODEML'} $codeml_ctrl_file");
	
	# codeml returns:
	my $results_file = 'results.txt';
	my $twoNGdN_file = '2NG.dN';
	my $twoNGdS_file = '2NG.dS';
	my $twoNGtree_file = '2NG.t';
	my $fourfold_nuc_file = '4fold.nuc';
	my $lnf_file = 'lnf';
	my $rst_file = 'rst';
	my $rst1_file = 'rst1';
	my $rub_file = 'rub';
	
	return ($results_file,$twoNGdN_file,$twoNGdS_file,$twoNGtree_file,$fourfold_nuc_file,$lnf_file,$rst_file,$rst1_file,$rub_file);
}

sub get_mean_boot4WAG_dist_tree
{
	# returns the mean bootstrap support reported by phyml using a WAG model

	my ($phylip_file,$ref_FASTA_hash,$clean) = @_;

	my $treefile = $phylip_file . "_phyml_tree.txt";
	my $lkfile =   $phylip_file . "_phyml_lk.txt";
	my $statfile = $phylip_file . "_phyml_stat.txt";
	my $bootstatfile = $phylip_file . "_phyml_boot_stats.txt";
	my $boottreefile = $phylip_file . "_phyml_boot_trees.txt";
   	my $tmpfile = (split(/\./,$phylip_file))[0] . "_tmp_log_phyml";
	
	my (@bootstrap_values,$newick,$labelled_newick,$boot);
	my ($n_of_splits,$mean_boot,$median_boot) = (0,0.0,0.0);

   	# run PHYML
   	system("$ENV{'EXE_PHYML'} $phylip_file 1 i 1 $ENV{'PHYML_WAG_BOOTSTRAP_FACTOR'} WAG 0.0 1 1 BIONJ n n > $tmpfile");
 
 	# parse tree
 	open(TREE,$treefile) || die "# get_mean_boot4WAG_dist_tree : cannot read $treefile\n";
	while(<TREE>){ $newick .= $_; }
	close(TREE);
	
	if(!$newick)
	{
		print "# get_mean_boot4WAG_dist_tree : empty boostrap tree $phylip_file\n";
		return (0.0,0.0,'failed');
	}
	
	# create labelled tree
	$labelled_newick = add_labels2newick_tree( $newick, $ref_FASTA_hash );
	open(LABELTREE,">$treefile") ||die "# get_mean_boot4WAG_dist_tree : cannot create $treefile\n";
	print LABELTREE $labelled_newick;
	close(LABELTREE); #print "#$newick#$labelled_newick#\n";

	while($newick =~ /\)(\d+):/g)
	{
		$boot = $1;
		push(@bootstrap_values,$boot);
		$mean_boot += $boot;
		$n_of_splits++;
	} 
	
	$mean_boot = sprintf("%1.0f",$mean_boot/$n_of_splits);
	$median_boot = sprintf("%1.0f",calc_median(@bootstrap_values));
	
   	# clean files
   	if($clean==1){ system(" rm -f $tmpfile $lkfile $statfile $bootstatfile $boottreefile"); }

   	return ($mean_boot,$median_boot,$treefile); 
}

sub get_mean_aLRT4WAGG_tree
{
	# returns the mean aLRT reported by phyml using a WAG model

	my ($phylip_file,$ref_FASTA_hash,$clean) = @_;

	my (@aLRT_values,$newick,$labelled_newick,$aLRT);
	my ($n_of_splits,$mean_aLRT,$median_aLRT) = (0,0.0,0.0);

	my ($treefile,$lkfile,$statfile) = run_PHYML_WAG_G($phylip_file,$clean);
 
 	# parse tree
 	open(TREE,$treefile) || die "# get_mean_aLRT4WAGG_tree : cannot read $treefile\n";
	while(<TREE>){ $newick .= $_; }
	close(TREE);
	
	# create labelled tree
	$labelled_newick = add_labels2newick_tree( $newick, $ref_FASTA_hash );
	open(LABELTREE,">$treefile") ||die "# get_mean_aLRT4WAGG_tree : cannot create $treefile\n";
	print LABELTREE $labelled_newick;
	close(LABELTREE); #print "#$newick#$labelled_newick#\n";

	while($newick =~ /\)(\d+\.\d+):/g)
	{
		$aLRT = sprintf("%1.2f",$1);
		push(@aLRT_values,$aLRT);
		$mean_aLRT += $aLRT;
		$n_of_splits++;
	} 
	
	$mean_aLRT = sprintf("%1.2f",$mean_aLRT/$n_of_splits);
	$median_aLRT = sprintf("%1.2f",calc_median(@aLRT_values));
	
   	# clean files
   	if($clean==1){ system(" rm -f $lkfile $statfile"); }

   	return ($mean_aLRT,$median_aLRT,$treefile); 	


}

sub calc_median
{
	my (@values) = @_;
	
	@values = sort {$a<=>$b} @values;
	my ($median,$n_of_values) = (0.0,scalar(@values));
	
	if($n_of_values%2)
	{ 
		$median = $values[int($n_of_values/2)]; 
	}
	else
	{ 
		$median = ($values[$n_of_values/2] + $values[($n_of_values/2)-1]) / 2; 
	}

	return $median;
}

sub run_PHYML_WAG_G
{
   my ($phylip_file,$clean) = @_;

   my $treefile = $phylip_file . "_phyml_tree.txt";
   my $lkfile =   $phylip_file . "_phyml_lk.txt";
   my $statfile = $phylip_file . "_phyml_stat.txt";
   my $tmpfile = (split(/\./,$phylip_file))[0] . "_tmp_log_phyml";

   # run PHYML
   system("$ENV{'EXE_PHYML'} $phylip_file 1 i 1 0 WAG 0 4 e BIONJ y y > $tmpfile");
 
   # clean file
   if($clean==1){ system(" rm -f $tmpfile"); }

   # check the output files are there
   if(!(-s $treefile) || !(-s $lkfile) || !(-s $lkfile)){ return ($phyTools::ERROR,$phyTools::ERROR,$phyTools::ERROR); }
   else{ return ($treefile,$lkfile,$statfile); }
}


sub get_PUZZLE_param_from_PHYML_stat
{
        my ($PHYML_stat) = @_;

        my $outfile = (split(/\./,$PHYML_stat))[0] . "_puzzle_params";
        my (@freqs,@rates,$model);
        my ($alpha_parameter,$p_invariants,$ti_tv_ratio,$ti_tv_ratio_pur,$ti_tv_ratio_pyr)=(0,0,0,0,0);

        # parse stat file
        open(STAT,$PHYML_stat) || die "# get_PUZZLE_param_from_PHYML_stat : cannot read $PHYML_stat\n";
        while(<STAT>)
        {
            if(/\. Model of nucleotides substitution : (.+)/)
            {
                $model = $1;
            }
            if(/ - Gamma shape parameter : (.+)/)
            {
                $alpha_parameter = $1;
            }
            if(/\. Proportion of invariant : (.+)/)
            {
                $p_invariants = $1;
            }
            if(/\. Transition\/transversion ratio for purines : (.+)/)
            {
                $ti_tv_ratio_pur = $1;
            }
            if(/\. Transition\/transversion ratio for pyrimidines : (.+)/)
            {
                $ti_tv_ratio_pyr = $1;
            }
            if(/\. Transition\/transversion ratio : (.+)/)
            {
                $ti_tv_ratio = $1;
            }


            if(/  - f\(/)
            {
                push(@freqs,(100*(split)[2]));
            }
            if(/ <-> /)
            {
                push(@rates,(split)[3]);
            }

        }
        close(STAT);

        # write params file
        open(PARAM,">$outfile") || die "# get_PUZZLE_param_from_PHYML_stat : cannot write to $outfile\n";
        print PARAM "b\n";
        print PARAM "f\n";
        print PARAM "$freqs[0]\n";
        print PARAM "$freqs[1]\n";
        print PARAM "$freqs[2]\n";


       if($model eq 'GTR' or $model eq 'custom')
       {
           print PARAM "m\nm\n";
           print PARAM "1\n$rates[0]\n";
           print PARAM "2\n$rates[1]\n";
           print PARAM "3\n$rates[2]\n";
           print PARAM "4\n$rates[3]\n";
           print PARAM "5\n$rates[4]\n";
           #print PARAM "6\n";
       }
       elsif($model eq 'HKY' or $model eq 'K2P' or $model eq 'F81')
       {
           print PARAM "t\n$ti_tv_ratio\n";
       }
       elsif($model eq 'TN93')
       {
           print PARAM "m\n";
       }
        elsif($model eq 'JC69')
       {
           print PARAM "t\n0.5\n";
       }

       if($alpha_parameter and $p_invariants)
       {
           print PARAM "w\nw\nw\n"; ### i\na\n$alpha_parameter\n
       }
       elsif($alpha_parameter)
       {
           print PARAM "w\n"; ### a\n$alpha_parameter\n
       }
       elsif($p_invariants)
       {
           print PARAM "w\nw\n";  ### i\n$p_invariants\n
       }
       else
       {
           print PARAM "w\n";
       }

        print PARAM "y\n";
        close(PARAM);

        # check the output file is there
        if(!@rates and @freqs and ($model ne 'HKY' && $model ne 'K2P' && $model ne 'F81' && $model ne 'JC69' && $model ne 'TN93')){ return $phyTools::ERROR; }
	
        else{ return $outfile; }
}


sub get_PUZZLE_param_from_PHYML3_stat
{
        my ($PHYML_stat) = @_;

        my $outfile = (split(/\./,$PHYML_stat))[0] . "_puzzle_params";
        my (@freqs,@rates,$model);
        my ($alpha_parameter,$p_invariants,$ti_tv_ratio,$ti_tv_ratio_pur,$ti_tv_ratio_pyr)=(0,0,0,0,0);

        # parse stat file
        open(STAT,$PHYML_stat) || die "# get_PUZZLE_param_from_PHYML_stat : cannot read $PHYML_stat\n";
        while(<STAT>)
        {
            if(/\. Model of nucleotides substitution: 	(.+)/)
            {
                $model = $1;
            }
            if(/ - Gamma shape parameter: 		(.+)/)
            {
                $alpha_parameter = $1;
            }
            if(/\. Proportion of invariant:              (.+)/)
            {
                $p_invariants = $1;
            }
            if(/\. Transition\/transversion ratio for purines:                    (.+)/)
            {
                $ti_tv_ratio_pur = $1;
            }
            if(/\. Transition\/transversion ratio for pyrimidines:                        (.+)/)
            {
                $ti_tv_ratio_pyr = $1;
            }
            if(/\. Transition\/transversion ratio: 	(.+)/)
            {
                $ti_tv_ratio = $1;
            }


            if(/  - f\(/)
            {
                push(@freqs,(100*(split)[2]));
            }
            if(/ <-> /)
            {
                push(@rates,(split)[3]);
            }
        }
        close(STAT);

        # write params file
        open(PARAM,">$outfile") || die "# get_PUZZLE_param_from_PHYML3_stat : cannot write to $outfile\n";
        print PARAM "b\n";
        print PARAM "f\n";
        print PARAM "$freqs[0]\n";
        print PARAM "$freqs[1]\n";
        print PARAM "$freqs[2]\n";


       if($model eq 'GTR' or $model eq 'custom')
       {
           print PARAM "m\nm\n";
           print PARAM "1\n$rates[0]\n";
           print PARAM "2\n$rates[1]\n";
           print PARAM "3\n$rates[2]\n";
           print PARAM "4\n$rates[3]\n";
           print PARAM "5\n$rates[4]\n";
           #print PARAM "6\n";
       }
       elsif($model eq 'HKY' or $model eq 'K2P' or $model eq 'F81')
       {
           print PARAM "t\n$ti_tv_ratio\n";
       }
       elsif($model eq 'TN93')
       {
           print PARAM "m\n";
       }
        elsif($model eq 'JC69')
       {
           print PARAM "t\n0.5\n";
       }

       if($alpha_parameter and $p_invariants)
       {
           print PARAM "w\nw\nw\n"; ### i\na\n$alpha_parameter\n
       }
       elsif($alpha_parameter)
       {
           print PARAM "w\n"; ### a\n$alpha_parameter\n
       }
       elsif($p_invariants)
       {
           print PARAM "w\nw\n";  ### i\n$p_invariants\n
       }
       else
       {
           print PARAM "w\n";
       }

        print PARAM "y\n";
        close(PARAM);

        # check the output file is there
        if(!@rates and @freqs and ($model ne 'HKY' && $model ne 'K2P' && $model ne 'F81' && $model ne 'JC69' && $model ne 'TN93')){ return $phyTools::ERROR; }
	
        else{ return $outfile; }
}


sub run_PUZZLE_DIST
{
	# runs puzzle with a phylip-format sequence alignment to calculate a WAG distance
	# matrix with gamma rates
	
	my ($prot_phylip_file, $clean) = @_;
	
	my ($length,$n_of_sequences,$perc_constant_sites,$perc_gaps,$min_dist,$mean_dist,$max_dist,$alpha) = (0,0,0,0,0,0,0,0.0);
	
	my $distfile =   $prot_phylip_file . ".dist";
	my $puzzlefile = $prot_phylip_file . ".puzzle";
	my $tmpfile = (split(/\./,$prot_phylip_file))[0] . "_tmp_log_puzzle";
	
	
	
	# run PUZZLE
   	system("cat $ENV{'PUZZLE_PROT_WAG_DIST_PARAMS_FILE'} | $ENV{'EXE_PUZZLE'} $prot_phylip_file 1> $tmpfile 2>$tmpfile"); # NOTE: do not use &> or perl wont wait for system child proc
	
	# capture alignments stats
	open(PUZZLE,$puzzlefile) || die "# run_PUZZLE_DIST : cannot read $puzzlefile\n";
	while(<PUZZLE>)
	{
		if(/^Input data: (\d+) sequences with (\d+) /)
		{
			($n_of_sequences,$length) = ($1,$2); 
		}      
		elsif(/^Number of constant sites: \d+ \(= (\d+\.\d+)/)
		{
			$perc_constant_sites = $1;
		}
		elsif(/^ Sum\s+\d+\s+\d+\s+\d+\s+(\d+\.\d+)/)
		{
			$perc_gaps = $1;
		}
		elsif(/^Average distance \(over all possible pairs of sequences\):\s+(\d+\.\d+)/)
		{
			$mean_dist = $1;
		}
		elsif(/^\s+minimum\s+: (\d+\.\d+),\s+maximum\s+: (\d+\.\d+)/)
		{
			($min_dist,$max_dist) = ($1,$2);
		}
		elsif(/^Gamma distribution parameter alpha: (\d+\.\d+)/)
		{
			$alpha = $1;
		}
		elsif(/^Gamma distribution parameter alpha \(estimated from data set\): (\d+\.\d+)/) # special case
		{
			$alpha = $1;
		}
	}
	close(PUZZLE);
	
   	# clean files
   	if($clean==1){ system(" rm -f $tmpfile $puzzlefile "); }
	
   	# check the output files are there
   	if(!-s $distfile){ return ($phyTools::ERROR,$phyTools::ERROR); }
   	else{ return ($distfile,$alpha,"# alignment stats: length = $length \%gaps = $perc_gaps \%constant = $perc_constant_sites","# mean distance = $mean_dist"); }
}


sub run_PUZZLE_DIST_HKYG
{
	# runs puzzle with a phylip-formatted sequence alignment to calculate a HKY85+G distance
	# matrix for DNA sequences (with gamma-distributes rates over sites)
	
	my ($dna_phylip_file, $clean) = @_;
	
	my ($length,$n_of_sequences,$perc_constant_sites,$perc_gaps,$min_dist,$mean_dist,$max_dist,$alpha,$titv_ratio) = (0,0,0,0,0,0,0,0,0);
	
	my $distfile =   $dna_phylip_file . ".dist";
	my $puzzlefile = $dna_phylip_file . ".puzzle";
	my $tmpfile = (split(/\./,$dna_phylip_file))[0] . "_tmp_log_puzzle";
	
	
	# run PUZZLE
   	system("cat $ENV{'PUZZLE_DNA_HKY_DIST_PARAMS_FILE'} | $ENV{'EXE_PUZZLE'} $dna_phylip_file 1> $tmpfile 2> $tmpfile");
	
	# capture alignments stats
	open(PUZZLE,$puzzlefile) || die "# run_PUZZLE_DIST : cannot read $puzzlefile\n";
	while(<PUZZLE>)
	{
		if(/^Input data: (\d+) sequences with (\d+) /)
		{
			($n_of_sequences,$length) = ($1,$2); 
		}      
		elsif(/^Number of constant sites: \d+ \(= (\d+\.\d+)/)
		{
			$perc_constant_sites = $1;
		}
		elsif(/^Expected transition\/transversion ratio:\s+(\d+\.\d+)/)
		{
		       $titv_ratio = $1;
		}
		elsif(/^ Sum\s+\d+\s+\d+\s+\d+\s+(\d+\.\d+)/)
		{
			$perc_gaps = $1;
		}
		elsif(/^Average distance \(over all possible pairs of sequences\):\s+(\d+\.\d+)/)
		{
			$mean_dist = $1;
		}
		elsif(/^\s+minimum\s+: (\d+\.\d+),\s+maximum\s+: (\d+\.\d+)/)
		{
			($min_dist,$max_dist) = ($1,$2);
		}
		elsif(/^Gamma distribution parameter alpha: (\d+\.\d+)/)
		{
			$alpha = $1;
		}
		elsif(/^Gamma distribution parameter alpha \(estimated from data set\): (\d+\.\d+)/) # special case
		{
			$alpha = $1;
		}
	}
	close(PUZZLE);
	
   	# clean files
   	if($clean==1){ system(" rm -f $tmpfile $puzzlefile "); }
	
   	# check the output files are there
   	if(!-s $distfile){ return ($phyTools::ERROR,$phyTools::ERROR); }
   	else{ return ($distfile,$alpha,"# alignment stats: length = $length \%gaps = $perc_gaps \%constant = $perc_constant_sites","# mean distance = $mean_dist ti-tv-ratio = $titv_ratio"); }
}

sub run_PUZZLE_LM
{
        my ($DNA_phylip_file,$PUZZLEparams, $clean) = @_;

        my $epsfile =    $DNA_phylip_file . ".eps";
        my $distfile =   $DNA_phylip_file . ".dist";
        my $puzzlefile = $DNA_phylip_file . ".puzzle";
	my $LMpuzzlefile = $DNA_phylip_file . "_LM.puzzle";
        my $tmpfile = (split(/\./,$DNA_phylip_file))[0] . "_tmp_log_puzzle";

	my ($alpha,$n_of_sequences,$length,$perc_constant_sites,$perc_gaps);
	my ($mean_dist,$min_dist,$max_dist,$n_of_site_patterns,$perc_site_patterns);


   	# run PUZZLE
   	print "# run_PUZZLE_LM running: system(cat $PUZZLEparams | $ENV{'EXE_PUZZLE'} $DNA_phylip_file &> $tmpfile)\n";
	system("cat $PUZZLEparams | $ENV{'EXE_PUZZLE'} $DNA_phylip_file 1> $tmpfile 2> $tmpfile");
	print "# run_PUZZLE_LM running finished\n";
         
	unless(-s $puzzlefile)
	{  # WATCH OUT: need to sleep a bit in order for the file to be written by puzzle; this I/O operation is too slow because 64bit linux is faster than 32 bits puzzle!!!
           #		without sleeping $puzzlefile is not found, at least on Tenerife!!!! Cause of a difficult2catch bug
	    sleep(1); 
	}

	# capture alignments stats
	$alpha = 0.0;
	open (PUZZLE, '<', "$puzzlefile") or die "# ERROR: run_PUZZLE_LM : cannot read $puzzlefile : $!\n";
	while(<PUZZLE>)
	{
		if(/^Input data: (\d+) sequences with (\d+) /)
		{
			($n_of_sequences,$length) = ($1,$2); 
		}      
		elsif(/^Number of constant sites: \d+ \(= (\d+\.\d+)/)
		{
			$perc_constant_sites = $1;
		}
		elsif(/^Number of site patterns: (\d+)/)
		{
			$n_of_site_patterns = $1;
			$perc_site_patterns = sprintf("%1.1f",100 * ($n_of_site_patterns / $length));
		}
		elsif(/^ Sum\s+\d+\s+\d+\s+\d+\s+(\d+\.\d)/)
		{
			$perc_gaps = $1;
		}
		elsif(/^Average distance \(over all possible pairs of sequences\):\s+(\d+\.\d+)/)
		{
			$mean_dist = $1;
		}
		elsif(/^\s+minimum\s+: (\d+\.\d+),\s+maximum\s+: (\d+\.\d+)/)
		{
			($min_dist,$max_dist) = ($1,$2);
		}
		elsif(/^Gamma distribution parameter alpha: (\d+\.\d+)/ || /Gamma distribution parameter alpha \(estimated from data set\): (\d+\.\d+)/)
		{
			$alpha = $1;
		}
	}
	close(PUZZLE);


   # clean files
   if($clean==1){ system(" rm -f $tmpfile $epsfile "); }
	
	# rename outfile
	system("mv $puzzlefile $LMpuzzlefile");
	
   # check the output files are there
   if(!(-s $distfile) || !(-s $LMpuzzlefile)){ return ($phyTools::ERROR,$phyTools::ERROR); }
   else
   { 
   	return ($distfile,$LMpuzzlefile,
		$alpha,
		"length = $length \%gaps = $perc_gaps \%constant = $perc_constant_sites \%site_patterns = $perc_site_patterns",
		"$min_dist <= $mean_dist <= $max_dist"); 
   }
}


sub parse_LM_PUZZLE_file
{
	my ($PUZZLE_LM_file) = @_;
   
	my $DB_row = '';
	my %LM_data = ();
	my $Alignment_ID = 0;
	my $ChiSQ_failed_counter = 0;
	
	$LM_data{'ChiSQ_failed'} = 0;
	$LM_data{'Acc_No_ChiSQ_perc'} = '';
	
	open(LM,$PUZZLE_LM_file) || die "# parse_LM_PUZZLE_file : cannot read $PUZZLE_LM_file\n";
	while(<LM>)
	{
		# INPUT FILE
		if(/^Input file name: (\S+\.?\S+)/) # ok
		{
		 	$Alignment_ID++;
			$LM_data{$Alignment_ID} = $1;
		 	#print "this is the Alignemnt_ID $Alignment_ID => $LM_data{$Alignment_ID}\n\n"
		}
		# SEQUENCE ALIGNMENT
		if(/^Input data: (\d+) sequences with (\d+)/) # ok
		{ 
			$LM_data{'n_of_sequences'} = $1; 
			$LM_data{'nt_sites'} = $2; 
		}
		elsif(/^Number of constant sites: \d+ \(= (\d+\.\d+)/) # ok
		{ 
			$LM_data{'percent_constant_sites'} = $1; 
		}
		elsif(/Number of site patterns: (\d+)/) # ok
		{
		    $LM_data{'n_of_site_patterns'} = $1;
		}
		elsif(/Number of constant site patterns: \d+ \(= (\d+.\d+)/) # ok
		{
		    $LM_data{'percent_constant_site_patterns'} = $1;
		   # print "percent_constant_site_patterns: $LM_data{'percent_constant_site_patterns'}\n";
		}
		
		
		# AMBIGUOUS CHARACTERS...
		elsif(/^ Sum\s+\d+\s+\d+\s+\d+\s+(\d+\.\d+)/) # ok
		{
			$LM_data{'percent_gaps'} = $1;
		} 
		
		# SEQUENCE COMPOSITION...
		elsif(/^ (\S+)\s+failed\s+(\d+\.\d+)/)
		{     
		      
			$ChiSQ_failed_counter++;
			$LM_data{'ChiSQ_failed'} = $ChiSQ_failed_counter;
			$LM_data{'Acc_No_ChiSQ_perc'} .= $1."_".$2.":";
			#$LM_data{'Acc_No'} = $1;
			#$LM_data{'P_value'} = $3;
			# print "the chi2_comp value is => $LM_data{'chi2_comp'}\n\n";
			#print "the ChiSQ_failed value is => $LM_data{'ChiSQ_failed'}\n\n";
		}
		
		# LIKELIHOOD MAPPING STATS...
		elsif(/^Number of quartets in region 1: \d+ \(= (\d+\.\d+)/)
		{
			$LM_data{'percent_reg1_quartets'} = $1;
		}
		elsif(/^Number of quartets in region 2: \d+ \(= (\d+\.\d+)/)
		{
		  $LM_data{'percent_reg2_quartets'} = $1;
		}
		elsif(/^Number of quartets in region 3: \d+ \(= (\d+\.\d+)/)
		{
		  $LM_data{'percent_reg3_quartets'} = $1;
		}
		elsif(/^Number of quartets in region 4: \d+ \(= (\d+\.\d+)/)
		{
			$LM_data{'percent_reg4_quartets'} = $1;
		}
		elsif(/^Number of quartets in region 5: \d+ \(= (\d+\.\d+)/)
		{
		  $LM_data{'percent_reg5_quartets'} = $1;
		}
		elsif(/^Number of quartets in region 6: \d+ \(= (\d+\.\d+)/)
		{
		  $LM_data{'percent_reg6_quartets'} = $1;
		}
		elsif(/^Number of quartets in region 7: \d+ \(= (\d+\.\d+)/)
		{
		  $LM_data{'percent_nonresolved_quartets'} = $1;
		}
	}
	close(LM); 
	
	if($LM_data{'percent_reg1_quartets'} && 	$LM_data{'percent_reg2_quartets'} && 
		$LM_data{'percent_reg3_quartets'})
	{
		$LM_data{'percent_fullyresolved_quartets'} = $LM_data{'percent_reg1_quartets'} +
		$LM_data{'percent_reg2_quartets'} + $LM_data{'percent_reg3_quartets'};
	}
	else
	{
		$LM_data{'percent_fullyresolved_quartets'} = $phyTools::NERROR;
	}
	if($LM_data{'percent_reg4_quartets'} && $LM_data{'percent_reg5_quartets'} && 
		$LM_data{'percent_reg6_quartets'})
	{
		$LM_data{'percent_partiallyresolved_quartets'} = $LM_data{'percent_reg4_quartets'} +
		$LM_data{'percent_reg5_quartets'} + $LM_data{'percent_reg6_quartets'};
	}
	else
	{
		$LM_data{'percent_partiallyresolved_quartets'} = $phyTools::NERROR;
	}
	$LM_data{'normalized_FRQ'} = $LM_data{'percent_fullyresolved_quartets'}  * 100 / $LM_data{'nt_sites'} ;
	# print "percent_nonresolved_quartets: $LM_data{'percent_nonresolved_quartets'}\n";
	
     
     $DB_row = $Alignment_ID ."\t" . $LM_data{"$Alignment_ID"} ."\t" . $LM_data{'n_of_sequences'} ."\t" .
	  $LM_data{'nt_sites'} ."\t" . $LM_data{'percent_constant_sites'} ."\t" . $LM_data{'n_of_site_patterns'} ."\t" .
	  $LM_data{'percent_constant_site_patterns'} ."\t" . $LM_data{'percent_gaps'} ."\t" . 
	  $LM_data{'ChiSQ_failed'} ."\t" . "\'$LM_data{'Acc_No_ChiSQ_perc'}\'"."\t" . 
	  $LM_data{'percent_fullyresolved_quartets'} ."\t" .
	  $LM_data{'percent_partiallyresolved_quartets'} ."\t" . $LM_data{'percent_nonresolved_quartets'} ."\t" . 
	  $LM_data{'normalized_FRQ'};
    
    return $DB_row;
}	 



sub run_PROTTEST
{
	my ($AA_phylip_file,$REMOVETMPFILES) = @_;
	
	my ($pwd,$PROTTEST_bestAICmodel,$outfile,$skipped_models);
	
	$outfile = $AA_phylip_file . ".prottest";
	$pwd = `pwd`; chomp $pwd;

	$skipped_models = " -MtArt F -MtREV F -MtMam F -RtREV F -CpREV F -+F F ";

	# 1) copy input file and cd to ProtTest directory
	system("cp $AA_phylip_file $ENV{'PROTTEST_DIR'}");
	chdir($ENV{"PROTTEST_DIR"});
	
	# 2) run Prottest and parse output
	system("$ENV{'EXE_PROTTEST'} -i $AA_phylip_file -o $outfile $skipped_models -verbose F");
	
	open(PROTTEST,$outfile) || die "#run_PROTTEST : cannot read $outfile, perhaps prottest failed\n";
	while(<PROTTEST>)
	{	
		if(/^Best model according to AIC:\s+(\S+)/) 
		{ 
			$PROTTEST_bestAICmodel = $1; 
			$PROTTEST_bestAICmodel =~ s/\+//g;
			last;
		}
	}
	close(PROTTEST);	
	
	# 3) clean files
	if($REMOVETMPFILES){ system("rm -f $AA_phylip_file\_*"); }
	
	# 4) mv results file
	system("mv $AA_phylip_file* $pwd");
		
	# 4) cd to initial directory
	chdir($pwd);

	return ($PROTTEST_bestAICmodel,$outfile);
}
#--------------------------------------------------------------------------------------------#
sub run_PROTTEST2
{
	my ($AA_phylip_file,$REMOVETMPFILES,$PUZZLE_MODELS) = @_;
	
	my ($pwd,$PROTTEST_bestAICmodel,$outfile,$skipped_models);
	
	$outfile = $AA_phylip_file . ".prottest";
	$pwd = `pwd`; chomp $pwd;

	#$skipped_models = " -MtArt F -MtREV F -MtMam F -RtREV F -CpREV F -+F F ";
	
	if($PUZZLE_MODELS) # PUZZLE-compatible models
	{
	    $skipped_models = " -LG F -MtArt F -MtMam F -RtREV F -CpREV F -HIVb F -HIVw F -DCMut F -+F F"; # Pablo, Sept 28th, 2011. Adapted for phyml3.0.1 and ProtTest 2.4
	}
	else
	{
	   $skipped_models = " -MtArt F -MtMam F -RtREV F -Dayhoff F -HIVb F -HIVw F -DCMut F -+I F"; # Pablo, Sept 28th, 2011. Adapted for phyml3.0.1 and ProtTest 2.4
        }
	# 1) copy input file and cd to ProtTest directory
	system("cp $AA_phylip_file $ENV{'PROTTEST2_DIR'}");
	chdir($ENV{"PROTTEST2_DIR"});
	
	# 2) run Prottest and parse output
	system("$ENV{'EXE_PROTTEST2'} -i $AA_phylip_file -o $outfile $skipped_models -verbose F");
	
	open(PROTTEST,$outfile) || die "#run_PROTTEST2 : cannot read $outfile, perhaps prottest failed\n";
	while(<PROTTEST>)
	{	
		if(/^Best model according to AIC:\s+(\S+)/) 
		{ 
			$PROTTEST_bestAICmodel = $1; 
			$PROTTEST_bestAICmodel =~ s/\+//g;
			last;
		}
	}
	close(PROTTEST);	
	
	# 3) clean files
	if($REMOVETMPFILES){ system("rm -f $AA_phylip_file\_*"); }
	
	# 4) mv results file
	system("mv $AA_phylip_file* $pwd");
		
	# 4) cd to initial directory
	chdir($pwd);

	return ($PROTTEST_bestAICmodel,$outfile);
}

sub getPHYMLparameters_prot
{
	my ($model) = @_;
	
	# borrowed from MrAIC
	my (%MODELnr,%MODELcmd,$m);
	
	$MODELnr{0}	=	"WAG";
	$MODELnr{1}	=	"WAGG";
	$MODELnr{2}	=	"WAGIG";
	$MODELnr{3}	=	"WAGI";
	$MODELnr{4}	=	"JTT";
	$MODELnr{5}	=	"JTTG";
	$MODELnr{6}	=	"JTTIG";
	$MODELnr{7}	=	"JTTI";
	$MODELnr{8}	=	"VT";
	$MODELnr{9}	=	"VTG";
	$MODELnr{10}	=	"VTIG";
	$MODELnr{11}	=	"VTI";
	$MODELnr{12}	=	"Blosum62";
	$MODELnr{13}	=	"Blosum62G";
	$MODELnr{14}	=	"Blosum62IG";
	$MODELnr{15}	=	"Blosum62I";
	$MODELnr{16}	=	"Dayhoff";
	$MODELnr{17}	=	"DayhoffG";
	$MODELnr{18}	=	"DayhoffIG";
	$MODELnr{19}	=	"DayhoffI";
	$MODELnr{20}	=	"MtREV";	
	$MODELnr{21}	=	"MtREVG";
	$MODELnr{22}	=	"MtREVIG";
	$MODELnr{23}	=	"MtREVI";
	
	$MODELcmd{0}	=	"D\nM\nM\nY\n";	# WAG
	$MODELcmd{1}	=	"D\nM\nM\nR\nC\n6\nA\nY\nY\n";	# WAGG
	$MODELcmd{2}	=	"D\nM\nM\nV\nY\nR\nC\n6\nA\nY\nY\n";	# WAGGI
	$MODELcmd{3}	=	"D\nM\nM\nV\nY\nY\n";	# WAGI
	$MODELcmd{4}	=	"D\nY\n"; # JTT
	$MODELcmd{5}	=	"D\nR\nC\n6\nA\nY\nY\n"; # JTTG
	$MODELcmd{6}	=	"D\nV\nY\nR\nC\n6\nA\nY\nY\n"; # JTTGI
	$MODELcmd{7}	=	"D\nV\nY\nY\n"; #JTTI
	$MODELcmd{8}	=	"D\nM\nM\nM\nM\nM\nM\nY\n";	# VT
	$MODELcmd{9}	=	"D\nM\nM\nM\nM\nM\nM\nR\nC\n6\nA\nY\nY\n";	# VTG
	$MODELcmd{10}	=	"D\nM\nM\nM\nM\nM\nM\nV\nY\nR\nC\n6\nA\nY\nY\n";	# VTGI
	$MODELcmd{11}	=	"D\nM\nM\nM\nM\nM\nM\nV\nY\nY\n";	# VTI
	$MODELcmd{12}	=	"D\nM\nM\nM\nM\nM\nM\nM\nY\n";	# BLOSUM62
	$MODELcmd{13}	=	"D\nM\nM\nM\nM\nM\nM\nM\nR\nC\n6\nA\nY\nY\n";	
	$MODELcmd{14}	=	"D\nM\nM\nM\nM\nM\nM\nM\nV\nY\nR\nC\n6\nA\nY\nY\n";	 
	$MODELcmd{15}	=	"D\nM\nM\nM\nM\nM\nM\nM\nV\nY\nY\n";	
	$MODELcmd{16}	=	"D\nM\nM\nM\nM\nM\nM\nM\nM\nM\nY\n";	# Dayhoff
	$MODELcmd{17}	=	"D\nM\nM\nM\nM\nM\nM\nM\nM\nM\nR\nC\n6\nA\nY\nY\n";	
	$MODELcmd{18}	=	"D\nM\nM\nM\nM\nM\nM\nM\nM\nM\nV\nY\nR\nC\n6\nA\nY\nY\n";	 
	$MODELcmd{19}	=	"D\nM\nM\nM\nM\nM\nM\nM\nM\nM\nV\nY\nY\n";	
	$MODELcmd{20}	=	"D\nM\nY\n";	# MtREV
	$MODELcmd{21}	=	"D\nM\nR\nC\n6\nA\nY\nY\n";	# MtREVG
	$MODELcmd{22}	=	"D\nM\nV\nY\nR\nC\n6\nA\nY\nY\n";	# MtREVGI
	$MODELcmd{23}	=	"D\nM\nV\nY\nY\n";	# MtREVI
	
	foreach $m (keys(%MODELnr))
	{
		if($MODELnr{$m} eq $model)
		{
			return $MODELcmd{$m};
		}
	}

	return $phyTools::ERROR;
}


sub getPHYML3_parameters_prot
{
	# NOTE 1: written by PV, Sep 28th, 2011, to adapt to phyml3.0.1 and ProtTest2.4
	# note: some models are excluded, as defined in &run_PROTTEST
	my ($model) = @_;
	
	my (%MODELnr,%MODELcmd,$m);
	
	$MODELnr{0}	=	"WAG";
	$MODELnr{1}	=	"WAGG";
	$MODELnr{2}	=	"WAGIG";
	$MODELnr{3}	=	"WAGI";
	$MODELnr{4}	=	"JTT";
	$MODELnr{5}	=	"JTTG";
	$MODELnr{6}	=	"JTTIG";
	$MODELnr{7}	=	"JTTI";
	$MODELnr{8}	=	"VT";
	$MODELnr{9}	=	"VTG";
	$MODELnr{10}	=	"VTIG";
	$MODELnr{11}	=	"VTI";
	$MODELnr{12}	=	"Blosum62";
	$MODELnr{13}	=	"Blosum62G";
	$MODELnr{14}	=	"Blosum62IG";
	$MODELnr{15}	=	"Blosum62I";
	$MODELnr{16}	=	"Dayhoff";
	$MODELnr{17}	=	"DayhoffG";
	$MODELnr{18}	=	"DayhoffIG";
	$MODELnr{19}	=	"DayhoffI";
	$MODELnr{20}	=	"MtREV";	
	$MODELnr{21}	=	"MtREVG";
	$MODELnr{22}	=	"MtREVIG";
	$MODELnr{23}	=	"MtREVI";
	$MODELnr{24}	=	"LG";
	$MODELnr{25}	=	"LGG";
	$MODELnr{26}	=	"LGIG";
	$MODELnr{27}	=	"LGI";
 
        # Models considering empirical AA frequencies
	$MODELnr{28}	=	"WAGF";
	$MODELnr{29}	=	"WAGGF";
	$MODELnr{30}	=	"WAGIGF";
	$MODELnr{31}	=	"WAGIF";
	$MODELnr{32}	=	"JTTF";
	$MODELnr{33}	=	"JTTGF";
	$MODELnr{34}	=	"JTTIGF";
	$MODELnr{35}	=	"JTTIF";
	$MODELnr{36}	=	"VTF";
	$MODELnr{37}	=	"VTGF";
	$MODELnr{38}	=	"VTIGF";
	$MODELnr{39}	=	"VTIF";
	$MODELnr{40}	=	"Blosum62F";
	$MODELnr{41}	=	"Blosum62GF";
	$MODELnr{42}	=	"Blosum62IGF";
	$MODELnr{43}	=	"Blosum62IF";
	$MODELnr{44}	=	"DayhoffF";
	$MODELnr{45}	=	"DayhoffGF";
	$MODELnr{46}	=	"DayhoffIGF";
	$MODELnr{47}	=	"DayhoffIF";
	$MODELnr{48}	=	"MtREVF";	
	$MODELnr{49}	=	"MtREVGF";
	$MODELnr{50}	=	"MtREVIGF";
	$MODELnr{51}	=	"MtREVIF";
	$MODELnr{52}	=	"LGF";
	$MODELnr{53}	=	"LGGF";
	$MODELnr{54}	=	"LGIGF";
	$MODELnr{55}	=	"LGIF";

	$MODELcmd{0}	=	"-d aa -m WAG -c 1 -v 0"; 
	$MODELcmd{1}	=	"-d aa -m WAG -c 4 -a e"; 
	$MODELcmd{2}	=	"-d aa -m WAG -c 4 -a e -v e"; 
	$MODELcmd{3}	=	"-d aa -m WAG -c 1 -v e"; 
	$MODELcmd{4}	=	"-d aa -m JTT -c 1 -v 0"; 
	$MODELcmd{5}	=	"-d aa -m JTT -c 4 -a e -v 0"; 
	$MODELcmd{6}	=	"-d aa -m JTT -c 4 -a e -v e"; 
	$MODELcmd{7}	=	"-d aa -m JTT -c 1 -v e"; 
	$MODELcmd{8}	=	"-d aa -m VT -c 1 -v 0"; 
	$MODELcmd{9}	=	"-d aa -m VT -c 4 -a e -v 0"; 
	$MODELcmd{10}	=	"-d aa -m VT -c 4 -a e -v e"; 
	$MODELcmd{11}	=	"-d aa -m VT -c 1 -v e"; 
	$MODELcmd{12}	=	"-d aa -m BLOSUM62 -c 1 -v 0"; 
	$MODELcmd{13}	=	"-d aa -m BLOSUM62 -c 4 -a e -v 0"; 
	$MODELcmd{14}	=	"-d aa -m BLOSUM62 -c 4 -a e -v e"; 
	$MODELcmd{15}	=	"-d aa -m BLOSUM62 -c 1 -v e"; 
	$MODELcmd{16}	=	"-d aa -m Dayhoff -c 1 -v 0"; 
	$MODELcmd{17}	=	"-d aa -m Dayhoff -c 4 -a e -v 0"; 
	$MODELcmd{18}	=	"-d aa -m Dayhoff -c 4 -a e -v e"; 
	$MODELcmd{19}	=	"-d aa -m Dayhoff -c 1 -v e"; 
	$MODELcmd{20}	=	"-d aa -m MtREV -c 1 -v 0"; 
	$MODELcmd{21}	=	"-d aa -m MtREV -c 4 -a e -v 0"; 
	$MODELcmd{22}	=	"-d aa -m MtREV -c 4 -a e -v e"; 
	$MODELcmd{23}	=	"-d aa -m MtREV -c 1 -v 0"; 
	$MODELcmd{24}	=	"-d aa -m LG -c 1 -v 0"; 
	$MODELcmd{25}	=	"-d aa -m LG -c 4 -a e -v 0"; 
	$MODELcmd{26}	=	"-d aa -m LG -c 4 -a e -v e"; 
	$MODELcmd{27}	=	"-d aa -m LG -c 1 -v 0"; 

        
	$MODELcmd{28}	=	"-d aa -m WAG -c 1 -v 0 -f e"; 
	$MODELcmd{29}	=	"-d aa -m WAG -c 4 -a e -f e"; 
	$MODELcmd{30}	=	"-d aa -m WAG -c 4 -a e -v e -f e"; 
	$MODELcmd{31}	=	"-d aa -m WAG -c 1 -v e -f e"; 
	$MODELcmd{32}	=	"-d aa -m JTT -c 1 -v 0 -f e"; 
	$MODELcmd{33}	=	"-d aa -m JTT -c 4 -a e -v 0 -f e"; 
	$MODELcmd{34}	=	"-d aa -m JTT -c 4 -a e -v e -f e"; 
	$MODELcmd{35}	=	"-d aa -m JTT -c 1 -v e -f e"; 
	$MODELcmd{36}	=	"-d aa -m VT -c 1 -v 0 -f e"; 
	$MODELcmd{37}	=	"-d aa -m VT -c 4 -a e -v 0 -f e"; 
	$MODELcmd{38}	=	"-d aa -m VT -c 4 -a e -v e -f e"; 
	$MODELcmd{39}	=	"-d aa -m VT -c 1 -v e -f e"; 
	$MODELcmd{40}	=	"-d aa -m BLOSUM62 -c 1 -v 0 -f e"; 
	$MODELcmd{41}	=	"-d aa -m BLOSUM62 -c 4 -a e -v 0 -f e"; 
	$MODELcmd{42}	=	"-d aa -m BLOSUM62 -c 4 -a e -v e -f e"; 
	$MODELcmd{43}	=	"-d aa -m BLOSUM62 -c 1 -v e -f e"; 
	$MODELcmd{44}	=	"-d aa -m Dayhoff -c 1 -v 0 -f e"; 
	$MODELcmd{45}	=	"-d aa -m Dayhoff -c 4 -a e -v 0 -f e"; 
	$MODELcmd{46}	=	"-d aa -m Dayhoff -c 4 -a e -v e -f e"; 
	$MODELcmd{47}	=	"-d aa -m Dayhoff -c 1 -v e -f e"; 
	$MODELcmd{48}	=	"-d aa -m MtREV -c 1 -v 0 -f e"; 
	$MODELcmd{49}	=	"-d aa -m MtREV -c 4 -a e -v 0 -f e"; 
	$MODELcmd{50}	=	"-d aa -m MtREV -c 4 -a e -v e -f e"; 
	$MODELcmd{51}	=	"-d aa -m MtREV -c 1 -v 0 -f e"; 
	$MODELcmd{52}	=	"-d aa -m LG -c 1 -v 0 -f e"; 
	$MODELcmd{53}	=	"-d aa -m LG -c 4 -a e -v 0 -f e"; 
	$MODELcmd{54}	=	"-d aa -m LG -c 4 -a e -v e -f e"; 
	$MODELcmd{55}	=	"-d aa -m LG -c 1 -v 0 -f e"; 
	
	foreach $m (keys(%MODELnr))
	{
		if($MODELnr{$m} eq $model)
		{
			return $MODELcmd{$m};
		}
	}

	return $phyTools::ERROR;
}


sub getPHYMLparameters_prot_100boot
{
	my ($model) = @_;
	
	# borrowed from MrAIC
	my (%MODELnr,%MODELcmd,$m);
	
	$MODELnr{0}	=	"WAG";
	$MODELnr{1}	=	"WAGG";
	$MODELnr{2}	=	"WAGIG";
	$MODELnr{3}	=	"WAGI";
	$MODELnr{4}	=	"JTT";
	$MODELnr{5}	=	"JTTG";
	$MODELnr{6}	=	"JTTIG";
	$MODELnr{7}	=	"JTTI";
	$MODELnr{8}	=	"VT";
	$MODELnr{9}	=	"VTG";
	$MODELnr{10}	=	"VTIG";
	$MODELnr{11}	=	"VTI";
	$MODELnr{12}	=	"Blosum62";
	$MODELnr{13}	=	"Blosum62G";
	$MODELnr{14}	=	"Blosum62IG";
	$MODELnr{15}	=	"Blosum62I";
	$MODELnr{16}	=	"Dayhoff";
	$MODELnr{17}	=	"DayhoffG";
	$MODELnr{18}	=	"DayhoffIG";
	$MODELnr{19}	=	"DayhoffI";
	$MODELnr{20}	=	"MtREV";	
	$MODELnr{21}	=	"MtREVG";
	$MODELnr{22}	=	"MtREVIG";
	$MODELnr{23}	=	"MtREVI";
	
	$MODELcmd{0}	=	"D\nB\n100\nN\nM\nM\nY\n";	# WAG
	$MODELcmd{1}	=	"D\nB\n100\nN\nM\nM\nR\nC\n6\nA\nY\nY\n";	# WAGG
	$MODELcmd{2}	=	"D\nB\n100\nN\nM\nM\nV\nY\nR\nC\n6\nA\nY\nY\n";	# WAGGI
	$MODELcmd{3}	=	"D\nB\n100\nN\nM\nM\nV\nY\nY\n";	# WAGI
	$MODELcmd{4}	=	"D\nB\n100\nN\nY\n"; # JTT
	$MODELcmd{5}	=	"D\nB\n100\nN\nR\nC\n6\nA\nY\nY\n"; # JTTG
	$MODELcmd{6}	=	"D\nB\n100\nN\nV\nY\nR\nC\n6\nA\nY\nY\n"; # JTTGI
	$MODELcmd{7}	=	"D\nB\n100\nN\nV\nY\nY\n"; #JTTI
	$MODELcmd{8}	=	"D\nB\n100\nN\nM\nM\nM\nM\nM\nM\nY\n";	# VT
	$MODELcmd{9}	=	"D\nB\n100\nN\nM\nM\nM\nM\nM\nM\nR\nC\n6\nA\nY\nY\n";	# VTG
	$MODELcmd{10}	=	"D\nB\n100\nN\nM\nM\nM\nM\nM\nM\nV\nY\nR\nC\n6\nA\nY\nY\n";	# VTGI
	$MODELcmd{11}	=	"D\nB\n100\nN\nM\nM\nM\nM\nM\nM\nV\nY\nY\n";	# VTI
	$MODELcmd{12}	=	"D\nB\n100\nN\nM\nM\nM\nM\nM\nM\nM\nY\n";	# BLOSUM62
	$MODELcmd{13}	=	"D\nB\n100\nN\nM\nM\nM\nM\nM\nM\nM\nR\nC\n6\nA\nY\nY\n";	
	$MODELcmd{14}	=	"D\nB\n100\nN\nM\nM\nM\nM\nM\nM\nM\nV\nY\nR\nC\n6\nA\nY\nY\n";	 
	$MODELcmd{15}	=	"D\nB\n100\nN\nM\nM\nM\nM\nM\nM\nM\nV\nY\nY\n";	
	$MODELcmd{16}	=	"D\nB\n100\nN\nM\nM\nM\nM\nM\nM\nM\nM\nM\nY\n";	# Dayhoff
	$MODELcmd{17}	=	"D\nB\n100\nN\nM\nM\nM\nM\nM\nM\nM\nM\nM\nR\nC\n6\nA\nY\nY\n";	
	$MODELcmd{18}	=	"D\nB\n100\nN\nM\nM\nM\nM\nM\nM\nM\nM\nM\nV\nY\nR\nC\n6\nA\nY\nY\n";	 
	$MODELcmd{19}	=	"D\nB\n100\nN\nM\nM\nM\nM\nM\nM\nM\nM\nM\nV\nY\nY\n";	
	$MODELcmd{20}	=	"D\nB\n100\nN\nM\nY\n";	# MtREV
	$MODELcmd{21}	=	"D\nB\n100\nN\nM\nR\nC\n6\nA\nY\nY\n";	# MtREVG
	$MODELcmd{22}	=	"D\nB\n100\nN\nM\nV\nY\nR\nC\n6\nA\nY\nY\n";	# MtREVGI
	$MODELcmd{23}	=	"D\nB\n100\nN\nM\nV\nY\nY\n";	# MtREVI
	
	foreach $m (keys(%MODELnr))
	{
		if($MODELnr{$m} eq $model)
		{
			return $MODELcmd{$m};
		}
	}

	return $phyTools::ERROR;
}



sub create_PUZZLE_LM_prot_parameters_file # Updated by PV on Oct 1st 2011, to read the outuput produced by phyml3.0.1
{
	my ($AA_phylip_file,$PHYML_stat,$PROTTEST_bestAICmodel) = @_;
	
	my ($myAlpha, $myprop_invariants,$GIparams) = ('','','');
	
	my $outfile = $AA_phylip_file . "_puzzle_LM_prot_params";
	
	open(PARAMS,">$outfile") || die "# create_PUZZLE_LM_prot_parameters_file : cannot write to $outfile\n";
	
	# 1) handle simple models 
	if($PROTTEST_bestAICmodel eq 'Dayhoff')
	{
		print PARAMS "b\nm\ny\n";
	}
	elsif($PROTTEST_bestAICmodel eq 'JTT')
	{
		print PARAMS "b\nm\nm\ny\n";
	}
	elsif($PROTTEST_bestAICmodel eq 'MtREV')
	{
		print PARAMS "b\nm\nm\nm\ny\n";
	}
	elsif($PROTTEST_bestAICmodel eq 'Blosum62')
	{
		print PARAMS "b\nm\nm\nm\nm\ny\n";
	}
	elsif($PROTTEST_bestAICmodel eq 'VT')
	{
		print PARAMS "b\nm\nm\nm\nm\nm\ny\n";
	}
	elsif($PROTTEST_bestAICmodel eq 'WAG')
	{
		print PARAMS "b\nm\nm\nm\nm\nm\nm\ny\n";
	}
	else # 2) in the following cases we must read gamma shape param and prop of invariant sites
	{
	
		# 2.1) read PHYL params
		open(PHYMLSTAT,$PHYML_stat) || die "# create_PUZZLE_LM_prot_parameters_file : cannot read $PHYML_stat\n";
		while(<PHYMLSTAT>)
		{
			#print "# using PHYML_stats file: $PHYML_stat\n";
			if(/  - Gamma shape parameter: 		(\d+.\d+)/){ $myAlpha = $1; }
			elsif(/ Proportion of invariant: 		(\d+.\d+)/){ $myprop_invariants = $1; }
		}
		close(PHYMLSTAT);
	
		# 2.2) choose proper PUZZLE command sequence for among site rate var (IG)
		if($PROTTEST_bestAICmodel =~ 'IG')
		{
			$GIparams = "w\nw\nw\ni\n$myprop_invariants\na\n$myAlpha\ny\n";
		}
		elsif($PROTTEST_bestAICmodel =~ 'I')
		{
			$GIparams = "w\nw\ni\n$myprop_invariants\ny\n";
		}
		elsif($PROTTEST_bestAICmodel =~ 'G')
		{
			$GIparams = "w\na\n$myAlpha\ny\n";
		}
	
		# 2.3) write complete PUZZLE command sequence
		if($PROTTEST_bestAICmodel =~ 'Dayhoff')
		{
			print PARAMS "b\nm\n$GIparams";
		}
		elsif($PROTTEST_bestAICmodel =~ 'JTT')
		{
			print PARAMS "b\nm\nm\n$GIparams";
		}
		elsif($PROTTEST_bestAICmodel =~ 'MtREV')
		{
			print PARAMS "b\nm\nm\nm\n$GIparams";
		}
		elsif($PROTTEST_bestAICmodel =~ 'Blosum62')
		{
			print PARAMS "b\nm\nm\nm\nm\n$GIparams";
		}		
		elsif($PROTTEST_bestAICmodel =~ 'VT')
		{
			print PARAMS "b\nm\nm\nm\nm\nm\n$GIparams";
		}
		elsif($PROTTEST_bestAICmodel =~ 'WAG')
		{
			print PARAMS "b\nm\nm\nm\nm\nm\nm\n$GIparams";
		}	
	}
	close(PARAMS);
	
	return $outfile;
}


sub put_fullname_in_topology
{
	my ($tree_topology , $number2name_equiv_file ) = @_;

	my ($fullname_topol,$dum,$num,$fullname,%EQUIV);
	
	# 1) read equiv file
	open(EQUIV,$number2name_equiv_file)|| die "# put_fullname_in_topology : cannot read $number2name_equiv_file\n";
	while(<EQUIV>)
	{
		next if not(/=>/);
		($dum,$num,$dum,$fullname) = split(/ /,$_,4);
		
		if($fullname =~ /\[(\S+)\]/){	$fullname = $1; }
		#elsif($fullname =~ /\[(\S+)\]/){	$fullname = $1; }
		
		$fullname =~ s/[,|\(|\)]//g;
		chomp $fullname;
		
		$EQUIV{$num} = $fullname;
		
		#print "$num $fullname\n";
	}
	close(EQUIV);
	
	# substitute numbers for fullnames in topologies
	$fullname_topol = $tree_topology;
	
	$fullname_topol =~ s/(\d+):/$EQUIV{$1}:/g;
	
	return $fullname_topol;
}


#sub oligo_megablast
#{
#	my ($blastDB,$blastinfile) = @_;
#
#	my (@blasthits,$hit,$refleft,$refright,$righthit,$lefthit);
#	my ($no_overlaps) = 0;
#
#	# 1) run megablast	
#	open(BLAST,"$ENV{'EXE_MEGABLAST'} -d $blastDB -i $blastinfile -W 15 -e $ENV{'MEGABLAST_EVAL'} -m 8 -D 2 -F F |") ||
#             die "# oligo_megablast : cannot run $ENV{'EXE_MEGABLAST'} -d $blastDB -i $blastinfile\n";
#        while(<BLAST>)
#        {
#        	push(@blasthits,[ split ]); #print;
#        }
#        close(BLAST);
#
#	# 1.1) results are not sorted by evalue!! 
#	@blasthits = sort {$a->[10]<=>$b->[10]} (@blasthits);
#
#        # 2) comprueba si hay hits solapantes
#        if(@blasthits)
#        {
#        	if($blasthits[0][8] > $blasthits[0][9])
#                {
#        	        $refright = $blasthits[0][8];
#                        $refleft  = $blasthits[0][9];
#                }
#                else
#                {
#                        $refright = $blasthits[0][9];
#                        $refleft  = $blasthits[0][8];
#                }
#
#                foreach $hit (@blasthits)
#                {
#                        $lefthit  = $hit->[9];
#                        $righthit = $hit->[8];
#                        if($righthit < $lefthit)
#                        {
#                	        $lefthit  = $hit->[8];
#                                $righthit = $hit->[9];
#                        }
#
#                        if($lefthit > $refright || $righthit < $refleft)
#                        {
#                        	$no_overlaps++;
#                        }
#		}
#	}
#        else{ $no_overlaps = $phyTools::NERROR; }
#
#	return ($no_overlaps,$blasthits[0][1],$refleft,$refright,$blasthits[0][10],$blasthits[0][2]);
#}


sub filter_aln_rate4site
{
	# takes a FASTA hash made by read_FASTA_sequence ($ref_FASTAaln)

	my ( $filename, $evol_model, $ref_FASTAaln , $evolrate_cutoff ) = @_;

	my ( %FASTAaln, $tmpFASTAfile, $tmpTREEfile, $treefile, $header, $origtree );
	my ( $tmpR4S, $tmpR4Sorig, $tmpR4Stree, %site_rates, $aln_mean, $aln_stddev );
	my ( $sequence, $site, $n_of_sites, @data, $rate4site_info, $removed_sites ); 

	# 1) print FASTA file padding names to have at least 6 chars '______'
	$tmpFASTAfile = $filename . "_r4s.faa";
	open(FASTA,">$tmpFASTAfile") || die "# filter_aln_rate4site : cannot write to $tmpFASTAfile\n";
	
	foreach $header (sort {$a<=>$b} (keys (%{$ref_FASTAaln}))) 
	{
		print FASTA ">$header\______\n$ref_FASTAaln->{$header}{'SEQ'}\n";
	}
	close(FASTA);

	# 2) read original tree and print a temporary file for rate4site with padded sequence names
	$treefile = $filename;
	#$treefile =~ s/\.faa/\.phy_phyml_tree\.txt/;
	$treefile =~ s/\.aa_amp/\_aln.phy_phyml_tree\.txt/;
	
	open(ORIGTREE,$treefile) || die "# $0 : cannot open to $treefile $!\n";
	while(<ORIGTREE>)
	{
		$origtree = $_;	
	}
	close(ORIGTREE);	

	while( $origtree =~ s/\((\d+):/\($1\______:/g ){ }
	while( $origtree =~ s/\,(\d+):/\,$1\______:/g ){ }	

	$tmpTREEfile = $treefile . '_r4s_tree.tmp';
        open(TREE,">$tmpTREEfile") || die "# filter_aln_rate4site : can't write to tmpTREEfile $!\n";
	print TREE $origtree;
	close(TREE);	

	# 3) run and  parse rate4site output and select fast evolving columns
	if($evol_model =~ /WAG/){ $evol_model = '-Mw'; }
	elsif($evol_model =~ /Dayhoff/){ $evol_model = '-Md'; }
        else{ $evol_model = '-Mj'; }

	$tmpR4Sorig = "r4sOrig.res";
	$tmpR4S     = "r4s.res";
	$tmpR4Stree = $filename . "_r4s_tree.out";

	system("$ENV{'EXE_RATE4SITE'} -s $tmpFASTAfile -t $tmpTREEfile -bn $evol_model -x $tmpR4Stree 1> /dev/null 2> /dev/null");
	
	open(R4SOUT,$tmpR4Sorig) || die "# filter_aln_rate4site : cannot read $tmpR4Sorig\n";
	while(<R4SOUT>)
	{
		if(/^#/)
		{
			if(/Average/){ $aln_mean = (split)[2] }
			if(/Standard Deviation/){ $aln_stddev = (split)[3] }
		}
		elsif(/^\s+\d+/)
		{
			@data = split;
			$site_rates{$data[0]} = $data[2];		
		}
	}
	close(R4SOUT);

	$rate4site_info = '';
	$removed_sites=0;
	foreach $header (sort {$a<=>$b} (keys (%{$ref_FASTAaln})))
        {
		$FASTAaln{$header}{'NAME'} = $ref_FASTAaln->{$header}{'NAME'};
		if($header == 1) 
		{
			$rate4site_info .= "R4S siterate= $aln_mean += $aln_stddev | removed sites: ";
		}

		$sequence = $ref_FASTAaln->{$header}{'SEQ'};
		$n_of_sites = length($sequence);
		
		for($site=0;$site<$n_of_sites;$site++)
		{
			if($site_rates{$site+1} <= $evolrate_cutoff)
			{
				$FASTAaln{$header}{'SEQ'} .= substr($sequence,$site,1);
			} 
			else
			{
				if($header == 1)
				{ 
					$rate4site_info .= ($site+1).","; 
					$removed_sites++;	
				}
			}
		}
        }
	$rate4site_info = "# $rate4site_info total_removed_sites: $removed_sites\n";
	
	# 4) clean temporary files
	system("rm -f $tmpFASTAfile $tmpTREEfile $tmpR4Sorig $tmpR4S $tmpR4Stree");	

	# 5) return hash with filtered sequences, removed columns recorded in FASTA headers 
	return ($rate4site_info,%FASTAaln);
}


sub get_GO_terms4gi 
{
	# returns a string with GO annotation for a GI and a species name
	my ($gi,$species) = @_;
	
	my (@GOfiles, @GOterms);
	my ($fileOK,$giOK) = (0,0);

	# 1) try to find GOfile for this species
	opendir(GODIR,$ENV{"GO_DIR"}) || die "# get_GO_annotation4gi : cannot list $ENV{'GO_DIR'}\n";
	@GOfiles = grep {/$species/} readdir(GODIR); 
	closedir(GODIR);

	if(!$GOfiles[0]){ print "# get_GO_annotation4gi : cannot find GO file for species $species\n";	}
	else{ $fileOK = 1 }

	# 2) find right gi in GOfile
	if($fileOK)
	{
		open(GOFILE,$ENV{"GO_DIR"}.$GOfiles[0])|| die "# get_GO_annotation4gi : cannot read GO file $GOfiles[0]\n";
		while(<GOFILE>)
		{
			if(/^$gi/)
			{
				while($_ =~ /(GO:\d+)/g){ push(@GOterms,$1); }
				$giOK = 1;
				last;
			}
		}
		close(GOFILE);
	}

	return @GOterms; 
}


sub translate_GO_terms
{
	my(@GOterms) = @_;

	my (@sentences,$term,$GOid);
	my $n_of_terms = scalar(@GOterms);
	my $terms_found = 0;

	open(GODICT,$ENV{"GO_DICTIONARY"})|| die "# translate_GO_terms : cannot read $ENV{'GO_DICTIONARY'}\n";
	while(<GODICT>)
        {
		next if(/^!/);
		$GOid = (split)[0];

		if(grep(/$GOid/,@GOterms))
		{
			push(@sentences,(split)[1]);
			++$terms_found;
			last if($terms_found == $n_of_terms);
		}	
	}
	close(GODICT);

	return @sentences;
}


sub calc_correlation_coef
{
	# calcula el coeficiente de correlacion de dos vectores X e Y
        # pasados por referencia

        my ($x,$y) = @_;

       	my ($sumx, $sumxs, $meanx, $meany, $sumxy, $sumy, $sumys);
        my ($i,$numfeat);

        $numfeat = scalar(@{$y});
        $sumy=$sumys=$sumx=$sumxs=$meanx=$meany=$sumxy=0;

        for ($i=0; $i < $numfeat; $i++)
        {
                $sumy += $y->[$i];
                $sumys += $$y[$i] * $$y[$i];
                $sumx += $$x[$i];
                $sumxs += $$x[$i] * $$x[$i];
                $meanx += $$x[$i];
                $meany += $$y[$i];
                $sumxy += $$y[$i] * $$x[$i];
        }
        $meanx /= $numfeat;
        $meany /= $numfeat;
                                                                                                             
        if ($sumx*$sumx/$numfeat == $sumxs  || $sumy*$sumy/$numfeat == $sumys)
        {
                return 0;
        }

        return ($sumxy - ($sumx*$sumy)/$numfeat)/sqrt( ($sumxs - ($sumx * $sumx)/$numfeat) * ($sumys - ($sumy* $sumy) /$numfeat));
}

sub find_taxa_FASTA_headers
{
        # takes a read_FASTA_sequence-returned hash ref and finds all present taxa, parsed as [strings]
        # expects [taxon name] to the righttmost side, and applies some naive filters
        # returns a hash with taxa names as keys and occurences as values

        my($ref_FASTA) = @_;

        my ($seq,$taxon,$header,%TAXA);

        foreach $seq (keys(%{$ref_FASTA}))
        {
                # assume last [] is the taxon name, to avoid problems with headers such as
                # [acyl-carrier-protein] synthase I [Buchnera aphidicola str. Bp (Baizongia pistaciae)]
                $header = reverse($ref_FASTA->{$seq}{'NAME'});
                while($header =~ /(\].+?\[)/g)
                {
                        $taxon = reverse($1);
			# adhoc rules to filter out common annotations such as:
                        # [5'-phosphate] , [3-hydroxymyristoyl] , 
                        # [acyl-carrier-protein] , [NAD(P)H] , [Cu-Zn]
                        next if(substr($taxon,1,1) !~ /[A-Z]/ || $taxon !~ /\w+\s+\w+/);

                        $TAXA{$taxon}{'SIZE'}++;
                        push(@{$TAXA{$taxon}{'MEMBERS'}},$seq);
			last;
                }
        }

        return %TAXA;
}

sub find_taxa_FASTA_array_headers
{
	# takes a read_FASTA_file_array-returned array ref and finds all present taxa, parsed as [strings]
	# expects [taxon name] to the righttmost side, and applies some naive filters
	# returns a hash with taxa names as keys and occurences as values
	# if $single_word_taxons_are_OK is passed, [labels] are parsed as valid taxon names, as opposed to
	# default [Genus species ...]
	my($ref_FASTA,$single_word_taxons_are_OK) = @_;
	
	my ($seq,$taxon,$header,%TAXA);
	
	foreach $seq (0 .. $#{ $ref_FASTA } )
	{
		# assume last [] is the taxon name, to avoid problems with headers such as
		# [acyl-carrier-protein] synthase I [Buchnera aphidicola str. Bp (Baizongia pistaciae)]
		$header = reverse($ref_FASTA->[$seq][NAME]); 
		while($header =~ /(\].+?\[)/g)
		{
			$taxon = reverse($1);
			# adhoc rules to filter out common annotations such as:
			# [5'-phosphate] , [3-hydroxymyristoyl] , 
			# [acyl-carrier-protein] , [NAD(P)H] , [Cu-Zn]
			
			next if(!$single_word_taxons_are_OK && 
				(substr($taxon,1,1) !~ /[A-Z]/ || $taxon !~ /\w+\s+\w+/));

			$TAXA{$taxon}{'SIZE'}++;
			push(@{$TAXA{$taxon}{'MEMBERS'}},$seq); 
			last;
		}
	}
	
	return %TAXA;
}


sub find_exon_boundaries_FASTA_headers
{
        # takes a read_FASTA_sequence-returned hash and finds all present exon boundadies, parsed as <<bound..ries,>>
        # returns a hash with fasta keys and boundaries as values in a array

        my($ref_FASTA) = @_;

        my ($seq,$boundaries,$right,$left,$ieb,%exon_boundaries);

        foreach $seq (sort {$a<=>$b}(keys(%{$ref_FASTA})))
        {
                if($ref_FASTA->{$seq}{'NAME'} =~ /<<(.+)>>/)
                {
                        $boundaries = $1;
			foreach $ieb (split(/,/,$boundaries))
			{
				if($ieb =~ /(\d+)\.\.(\d+)/)
				{ 
					($left,$right) = ($1,$2); 
					push(@{$exon_boundaries{$seq}},[$left,$right]);
				}
			}
                }
        }

        return %exon_boundaries;
}

sub find_introns_FASTA_headers
{
	my($ref_FASTA,$use_AA_coords) = @_;

	my ($seq,$intron,$ilength,$icoord,%introns);

	foreach $seq (sort {$a<=>$b} (keys(%{$ref_FASTA})))
        {
                if($ref_FASTA->{$seq}{'NAME'} =~ /==(.+)==/)
                {
			foreach $intron (split(/,/,$1))
			{
				if($intron =~ /(\d+)\((\d+)\)/)
				{
					($icoord,$ilength) = ($1,$2);
					if($use_AA_coords)
					{ 
						$icoord = int(($icoord/3)+0.5); 
					}
					$introns{$seq}{$icoord} = $ilength; 
				}
			}
		}
	}

	return %introns;
}


sub adjust_introns_alignment
{
	# adjust intron positions to alignment columns and returns a hash of alginment columns that span introns
	#if 1+ introns happen at the same colum, reports the longest intron size
	my($ref_FASTA,$ref_introns) = @_;

	my ($seq,$intron,$ilength,$icoord,%alnintrons);
	my ($pos,$n_of_cols);
	
	foreach $seq (sort {$a<=>$b} (keys(%{$ref_FASTA})))
        {
		$pos=$n_of_cols=0;
		foreach $intron (sort {$a<=>$b} (keys(%{$ref_introns->{$seq}})))
		{
			$ilength = $ref_introns->{$seq}{$intron};
			while($n_of_cols<$intron)
			{
				if(substr($ref_FASTA->{$seq}{'SEQ'},$pos,1) ne '-'){ $n_of_cols++ }
				$pos++;
			} #print "#$seq:$intron->$pos\n";

			if($alnintrons{$pos})
			{
				if($ilength > $alnintrons{$pos}){ $alnintrons{$pos} = $ilength }
			}
			else{ $alnintrons{$pos} = $ilength }	
		}
	}

	return %alnintrons;
}

sub translate_dna2prot
{
	# takes a dna sequence and translates it to protein using a passed genetic code
	# from the list of 16 available in Bioperl 
	# http://www.bioperl.org/Core/Latest/bptutorial.html
	# http://doc.bioperl.org/releases/bioperl-1.4/Bio/Tools/CodonTable.html
	# NO PARECE ACEPTAR TABLAS DE CODONES COMO ARGUMENTO con bioperl1.4 y algunas secuencias en particular (grandes.fna)
	
	my($dnaseq, $gencode) = @_;
	
	my $prot;

	my $seq = Bio::Seq->new( -display_id => '1', -seq => $dnaseq );
	if(!$gencode)
	{
		$prot = $seq->translate();
		#print "#prot1 is $prot ".$prot->seq()."\n";
	}
	else
	{
		$prot = $seq->translate(-codontable_id => $gencode); #print "#prot2 is ".$prot->seq()."\n";
		#$prot = $seq->translate(); print "#prot2 is ".$prot->seq()."\n";
	}
	
	return $prot->seq();
}


sub convert_FASTA_dna2prot
{
	# gets a DNA read_FASTA_sequence-returned hash and returns a translated hash of protein sequences
	# can take a non-standard code table/genetic code, NO PARECE ACEPTAR TABLAS DE CODONES COMO ARGUMENTO
	
	my ($ref_FASTA, $gencode) = @_;

	my ($seq,%FASTAprot);
	
	foreach $seq (keys(%{$ref_FASTA}))
	{
		$FASTAprot{$seq}{'NAME'} = $ref_FASTA->{$seq}{'NAME'};
		if($gencode eq '' || $gencode eq 'universal')
		{
			$gencode = 'universal';
			$FASTAprot{$seq}{'SEQ'} = translate_dna2prot($ref_FASTA->{$seq}{'SEQ'});
			#$FASTAprot{$seq}{'SEQ'} = (split(/\*/,$FASTAprot{$seq}{'SEQ'}))[0];
		}
		else
		{ 
			$FASTAprot{$seq}{'SEQ'} = translate_dna2prot($ref_FASTA->{$seq}{'SEQ'},$gencode);
			#$FASTAprot{$seq}{'SEQ'} = (split(/\*/,$FASTAprot{$seq}{'SEQ'}))[0];
		}
	}

	# debugging
	#foreach $seq (keys(%FASTAprot))
	#{
	#	print "$seq $FASTAprot{$seq}{'NAME'}$FASTAprot{$seq}{'SEQ'}\n";
	#}	
	
	return %FASTAprot;
}


sub cluster_sequences
{
	# clusters sequences based on their distance matrix
	# returns a printable matrix, a labelled newick tree to be displayed, and a hash of sequence clusters of size > 1
	# also returns a labelled graph tree and a scale string statin max,min distances 
	# cluster -1 includes all orphan sequences 

	my ($dist_file,$maxdist,$filter_redundant,$treefile,$graphtreefile,$ref_FASTAhash,%user_cluster) = @_;
	
	my (%clusters,$id,$n_of_sequences,$n_of_clusters,@distances);
	my (%IDs,$s1,$s2,$n_of_close_sequences,%clustered,$cluster);
	my ($n_of_c_sequences,%redundant,@tmpdistances,$dist_matrix);
	my ($fully_labelled_tree_file,$fully_labelled_tree,$label,$oldlabel);
	my ($fully_labelled_graph_tree_file,$fully_labelled_graph_tree);
	my (%scale,$scale_string);
	my ($tree_string,$member,$user_clusterOK) = ('','',0);
	
	# 0) check whether a user-defined cluster is to be used 
	if(scalar(values(%user_cluster)) > 0){ $user_clusterOK = 1 }
	
	# 1) read distance file
	$id = '';
	$n_of_sequences = -1;
	open(DIST,$dist_file) || die "# cluster_sequences: cannot read $dist_file\n";
	while(<DIST>)
	{
		next if(/^$/);
		if(/(^\S+)\s+/)
		{
			$id = $1;
			$n_of_sequences++;
			
			$IDs{$n_of_sequences} = $id;
			
			@tmpdistances = split;
			push(@{$distances[$n_of_sequences]}, @tmpdistances );
		}
		elsif($id)
		{
			@tmpdistances = split;
			push(@{$distances[$n_of_sequences]}, @tmpdistances);
		}
	}
	close(DIST);

	#print "#$distances[0][0]#\n";print "#$distances[1][3]#\n"; # primera etiqueta, tercer valor de segunda fila
	
	# 2) filter redundant sequences
	if(!$user_clusterOK && $filter_redundant)
	{
		for($s1=0;$s1<=$n_of_sequences;$s1++)
		{
			next if($redundant{$s1});	
			for($s2=1;$s2<=$n_of_sequences+1;$s2++)
			{
				next if($s1 == $s2-1);
					
				if($distances[$s1][$s2] <= $ENV{"MIN_PAIR_DIST_REDUND"})
				{
					#print "|$s1|$s2|$IDs{$s1}|",$IDs{$s2-1},"|$distances[$s1][$s2]|\n";
					$redundant{$s2-1} = 1;
				}
			}
		}		
	}
	
	
	# 3) compile printable distance matrix
	$dist_matrix = '######## ';
	for($s1=0;$s1<=$n_of_sequences;$s1++){ $dist_matrix .= sprintf("%04s ",$s1); }
	$dist_matrix .= "\n";
	for($s1=0;$s1<=$n_of_sequences;$s1++)
	{	
		$dist_matrix .= "$IDs{$s1} ";
		for($s2=1;$s2<=$n_of_sequences+1;$s2++)
		{
			#$dist_matrix .= sprintf("%1.1f ",$distances[$s1][$s2]);
			if($distances[$s1][$s2] <= $maxdist){ $dist_matrix .= sprintf("%1.2f ",$distances[$s1][$s2]); }
			else { $dist_matrix .= "     "; } 
			
			# 3.1) keep track of min,max matrix values
			if(!$scale{'max'} || $distances[$s1][$s2] > $scale{'max'}{'dist'})
			{
				$scale{'max'}{'dist'} = $distances[$s1][$s2];
				$scale{'max'}{'label'} = "$IDs{$s1} <=> $IDs{$s2-1}";
			}
			if($IDs{$s1} ne $IDs{$s2-1} && (!$scale{'min'} || $distances[$s1][$s2] < $scale{'min'}{'dist'}))
			{
				$scale{'min'}{'dist'} = $distances[$s1][$s2];
				$scale{'min'}{'label'} = "$IDs{$s1} <=> $IDs{$s2-1}";
			} #print "$distances[$s1][$s2] $IDs{$s1},$IDs{$s2-1}\n";
		}
		$dist_matrix .= "\n";
	}
	$dist_matrix .= "\n";
	
		
	# 4) check sequence distance and make clusters
	$n_of_clusters = 0;
	
	if(!$user_clusterOK)
	{
		for($s1=0;$s1<=$n_of_sequences;$s1++)
		{
			next if($clustered{$s1} || $redundant{$s1});
		
			$n_of_close_sequences = 0;
			
			for($s2=1;$s2<=$n_of_sequences+1;$s2++)
			{
				next if($s1 == $s2-1);
				next if($filter_redundant && $redundant{$s2-1});
				if(!$clustered{$s2-1} && $distances[$s1][$s2] <= $maxdist)
				{		
					push(@{$clusters{$n_of_clusters}}, $IDs{$s2-1});  
					$n_of_close_sequences++;
					
					$clustered{$s2-1} = 1;
				}
			}
	
			if($n_of_close_sequences)
			{ 
				push(@{$clusters{$n_of_clusters}}, $IDs{$s1});
				$clustered{$s1} = 1;
				$n_of_clusters++; 
			}
			else
			{
				push(@{$clusters{'-1'}}, $IDs{$s1});
			}
		}
	}
	else
	{
		$n_of_clusters++;
		%clusters = %user_cluster;
	}

	# 5) label tree leafs with cluster numbers
	$tree_string = '';
	open(TREE,$treefile) || die "# cluster_sequences: cannot read $treefile\n";
	while(<TREE>){ $tree_string .= $_; }
	close(TREE);	

	$tree_string =~ s/__\d+:/:/g; # remove previous cluster assignments from labels	#print "#$tree_string#\n"; exit;
	
	foreach $n_of_clusters (keys(%clusters))
	{	
		next if($n_of_clusters eq '-1'); # orphan sequences with no associated cluster 
		foreach $member (@{$clusters{$n_of_clusters}})
		{
			$tree_string =~ s/$member:/$member\__$n_of_clusters:/; 
		}
	}

	open(TREE,">$treefile") || die "# cluster_sequences: cannot rewrite $treefile\n";
	print TREE $tree_string;
	close(TREE);
	
	# 6) generate fully labelled tree
	$fully_labelled_tree = add_labels2newick_tree($tree_string,$ref_FASTAhash);	
	
	$fully_labelled_tree_file = $treefile;
	$fully_labelled_tree_file =~ s/nj\.ph/nj_labelled\.ph/;
	
	open(LABELTREE,">$fully_labelled_tree_file") || die "# cluster_sequences: cannot create $fully_labelled_tree_file\n";
	print LABELTREE $fully_labelled_tree;
	close(LABELTREE);
	
	# 7) generate fully labelled graph tree
	$tree_string = '';
	open(TREE,$graphtreefile) || die "# cluster_sequences: cannot read $graphtreefile\n";
	while(<TREE>){ $tree_string .= $_; }
	close(TREE);
	
	$tree_string =~ s/__\d+:/:/g; # remove previous cluster assignments from labels	#print "#$tree_string#\n"; exit;
	
	foreach $n_of_clusters (keys(%clusters))
	{	
		next if($n_of_clusters eq '-1'); # orphan sequences with no associated cluster 
		foreach $member (@{$clusters{$n_of_clusters}})
		{
			$tree_string =~ s/$member/$member\__$n_of_clusters/; 
		}
	}

	$fully_labelled_graph_tree = add_labels2newick_tree($tree_string,$ref_FASTAhash);
	
	$fully_labelled_graph_tree_file = $graphtreefile;
	$fully_labelled_graph_tree_file =~ s/nj_txt\.graph/nj_txt_labelled\.graph/;
	open(TREE,">$fully_labelled_graph_tree_file") || die "# cluster_sequences: cannot create $fully_labelled_graph_tree_file\n";
	print TREE $fully_labelled_graph_tree;
	close(TREE);
	
	# 8) write scale_string
	$scale_string = "# max distance = $scale{'max'}{'dist'} ( $scale{'max'}{'label'} )\n";
	$scale_string .= "# min distance = $scale{'min'}{'dist'} ( $scale{'min'}{'label'} )\n";
	
	return ($dist_matrix,$fully_labelled_tree_file,$fully_labelled_graph_tree_file,$scale_string,%clusters);	
}

sub add_labels2newick_tree
{
	# assumes a hash f hashes with {001}{NAME} = "ftz Ecoli"
	my ($fully_labelled_tree,$ref_label_hash) = @_;

	my ($label,$oldlabel,$taxon);

	foreach my $seq (keys(%{$ref_label_hash}))
        {
		$taxon = '';
		$label = $ref_label_hash->{$seq}{'NAME'}; 
                $label =~ s/^>//g;
                $label = (split(/#\d+\|/,$label))[0]; # ignora primers  #164|352|0|1|1|0|55.4|57.8|AGTCCAAG...
		$label =~ s/<<.*?>>//g;
		$label =~ s/\s+|;|\-|,|:|\*|\.|\(|\)|#/_/g;
		$label =~ s/_{1,}/_/g;
		$label =~ s/==\S*?==//g;
		$label =~ s/$seq\_//g; #print "mira:$label\n";
                if($label =~ /(\[.*?\])/){ $taxon = $1 }
		if($taxon && substr($label,0,100) !~ /$taxon/){ $label = substr($label,0,15) . '...' . $taxon; } #5May09 PV ~ if($taxon && substr($label,0,25) !~ /$taxon/){ $label = substr($label,0,15) . '...' . $taxon; }
                else{ $label = substr($label,0,100) . '...'; } #print  "mira:$label\n";

		if($fully_labelled_tree =~ /($seq\__\d+)/)
                {
                        $oldlabel = $1;
                        $fully_labelled_tree =~ s/$oldlabel/$oldlabel $label/;
                }
                elsif($fully_labelled_tree =~ /($seq):/)
                {
                        $oldlabel = $1;
                        $fully_labelled_tree =~ s/$oldlabel:/$oldlabel $label:/;
                }
	}

	#die $fully_labelled_tree;

	return $fully_labelled_tree;
}


sub select_codon_tables
{
	my ($ref_Fasta,$tabletype) = @_;
	
	my ($seq,$genus,$sp,$taxon,%species,@tmp,@codon_usage_files,$codon_list,$taxon_info,$tgenus,$tsp);
	my (%codondata,%genus_codondata,%codontables,$n_of_data,%tmpcodontables,%genus_tmpcodontables,$file);
	
	if(!$tabletype){ $tabletype = '' }
	
	# 1) get codon files
	opendir(DIR, $ENV{'GB_SPSUM_DIR'}) || die "# select_codon_tables: cannot list $ENV{'GB_SPSUM_DIR'}\n";
	@tmp = sort(grep { /\.spsum/ } readdir(DIR));
	closedir DIR;
	
	# 1.1) sort them and put phages and virues towards the end as they contain species names and give uncorrect tables
	while(scalar(@tmp)>2)
	{ 
		$file = shift(@tmp); 
		if($file =~ /phg/ || $file =~ /vrl/){ push(@tmp,$file) }
		else{ push(@codon_usage_files,$file) }
	}
	foreach $file (sort(@tmp)){ push(@codon_usage_files,$file) }#print "@codon_usage_files\n";

	
	# 2) identify suitable codon tables for each taxon in $ref_Fasta
	foreach $seq (sort {$a<=>$b} (keys(%{$ref_Fasta})))
	{
		$taxon_info = $ref_Fasta->{$seq}{'NAME'};
		$taxon_info =~ s/_/ /g;
		
		if($taxon_info =~ /\[\s?(\w+)\s+(\w+)/)
		{
			$genus = lc($1);
			$genus = uc(substr($genus,0,1)).substr($genus,1);
			$sp = lc($2);
			if($sp =~ 'spp'){ $sp = 'sp.' }
			
			$taxon = "$genus $sp";   #print "#$genus#$sp#\n";
			
			next if($species{$taxon});
			
			foreach $codon_list (@codon_usage_files)
			{
				#print "|$codon_list|$taxon|\n";
			
				open(LIST,$ENV{'GB_SPSUM_DIR'}.$codon_list) || die "# select_codon_tables: cannot read $codon_list\n";
				while(<LIST>)
				{
					next if not (/\:/);
					if($tabletype eq 'mt')
					{
						next unless(/mitochondrion/);
					}
					elsif($tabletype eq 'cp')
					{
						next unless(/chloroplast/);
					}
					elsif($tabletype eq 'pl')
					{
						next unless(/plastid/);
					}
					
					if(length($genus) > 1)
					{
						if(/$taxon.*\: (\d+)/)
						{
							$n_of_data = $1;
							if(!$codondata{$taxon} || $n_of_data > $codondata{$taxon})
							{
								#print "# $_";
								$codondata{$taxon} = $n_of_data;
								$tmpcodontables{$taxon} = (split(/\: /,$_))[0];
							}	
							$species{$taxon} = 1;
						}	
						elsif(length($genus) > 1 && /$genus.*\: (\d+)/)
						{
							#print "#$genus#$_";
							$n_of_data = $1;
							if(!$genus_codondata{$genus} || $n_of_data > $genus_codondata{$genus})
							{
								$genus_codondata{$genus} = $n_of_data;
								$genus_tmpcodontables{$genus} = (split(/\: /,$_))[0];
							}
							$species{$taxon} = 1;	
						}
					}			
					elsif(/(\w+)\s+(\w+).*\: (\d+)/)
					{
					       $tgenus = uc(substr($1,0,1));
					       $tsp = $2;
					       $n_of_data = $3; # print; print "#$tgenus#$tsp#\n";
					       
					       if($tsp ne 'sp' && $tsp eq $sp && $genus eq $tgenus) 
					       {		       
						       #print "$genus $sp <=> $tgenus $tsp\n";
						       if(!$codondata{$taxon} || $n_of_data > $codondata{$taxon})
						       {
							       $codondata{$taxon} = $n_of_data;
							       $tmpcodontables{$taxon} = (split(/\: /,$_))[0];
						       }
						       $species{$taxon} = 1;
					       }       
					}
				}
				close(LIST);
			
				next if(!$species{$taxon});
			
				#print "|$taxon|$tmpcodontables{$taxon}|$codondata{$taxon}|\n";
				#print "|$genus|$genus_tmpcodontables{$genus}|$genus_codondata{$genus}|\n";
			
				if($sp eq 'sp.')
				{
					if($tmpcodontables{$taxon} && $genus_tmpcodontables{$genus})
					{
						if($codondata{$taxon} > $genus_codondata{$genus})
						{
							$codontables{$tmpcodontables{$taxon}} = $codon_list;
						}
						else{ $codontables{$genus_tmpcodontables{$genus}} = $codon_list; }
					}
					elsif($tmpcodontables{$taxon})
					{
						$codontables{$tmpcodontables{$taxon}} = $codon_list;
					}
					elsif($genus_tmpcodontables{$genus})
					{
						$codontables{$genus_tmpcodontables{$genus}} = $codon_list;
					}
					last;
				}
				else
				{
					if($tmpcodontables{$taxon})
					{
						$codontables{$tmpcodontables{$taxon}} = $codon_list;
					}
					elsif($genus_tmpcodontables{$genus})
					{
						$codontables{$genus_tmpcodontables{$genus}} = $codon_list;
					}	
					last;
				}
				#print "@codontables\n";		
			}
		}
	}
	
	#print "#\n",keys(%codontables),"\n";exit;
	
	return %codontables;
}	

sub extract_intergenic_from_genbank
{
	# takes a genbank input file and creates a FNA file containing all intergenic sequences found
	# inspired by http://bioperl.org/pipermail/bioperl-l/2006-March/021065.html
	# use $min_intergenic_size = 0 to skip minimal length tests
	# use $max_intergenic_size = 0 to skip maximum length test
	# use $length_flanking_ORFs > 0 if you want to cut oligonucleotides from both flanking ORFs to be used as PCR anchors

	my ($infile,$out_intergenic_file,$min_intergenic_size,$max_intergenic_size,$length_flanking_ORFs) = @_;	

	if(!defined($length_flanking_ORFs) || $length_flanking_ORFs < 0){ $length_flanking_ORFs = 0 }

	my ($n_of_intergenic,$gi,$start,$end,$length,$strand,$genename,$taxon) = (0);
	my $in = new Bio::SeqIO(-file => $infile, -format => 'genbank' );

	# first create tmp outfile, in case sub does not finish cleanly
	open(FNA,">$out_intergenic_file.tmp") || die "# extract_intergenic_from_genbank : cannot create $out_intergenic_file.tmp\n";

	while( my $seq = $in->next_seq) # scan all sequences found in $input  
        {
		my ($gbaccession,$sequence,$gen,@genes) = ( $seq->accession() );
		$sequence = $seq->primary_seq()->seq() || 'empty, need a full genbank entry!';
		$taxon = '';

		for my $f ($seq->get_SeqFeatures) 
		{
			if($f->primary_tag =~ /CDS|rRNA|tRNA/) # campos de 'genes'
			{
				$gi = $genename = ''; # compatible con subrutina extract_CDS_from_genbank
				if($f->has_tag('db_xref'))
				{
					my $crossrefs = join(',',sort $f->each_tag_value('db_xref'));
					if($crossrefs =~ /(GI\:\d+)/){ $gi = $1 }
				}
				elsif($f->has_tag('locus_tag'))
				{
                               		if($gi eq '' && $f->has_tag('locus_tag')){ $gi = "ID:".join(',',sort $f->each_tag_value('locus_tag')); }
				}
				if($f->has_tag('gene'))
				{
					$genename = join(',',sort $f->each_tag_value('gene'));
				}				
				else
				{
					if($f->has_tag('product'))
                                        {
                                 		$genename = join(',',sort $f->each_tag_value('product'));
                                        }
				}

				next if($gi eq '');
 
				$start = $f->start(); # 1..n naturals
				$end   = $f->end();
				$strand = $f->location()->strand(); 
				push(@genes,[$start,$end,$gi,$strand,$genename]);
			}
			elsif($f->primary_tag() =~ /source/)
                        {
                                if($f->has_tag('organism'))
                                {
                                        foreach my $element ($f->each_tag_value('organism')){ $taxon .= "[$element],"; }
                                        chop $taxon;
                                }
			}
		}		

		# calcular ORFs vecinos y corta las secuencias intergenicas
		for($gen=1;$gen<scalar(@genes);$gen++)
		{
			next if($genes[$gen-1][1] >= $genes[$gen][0]); # evita solapar, asume genes estan en orden en GenBank
			next if( defined($min_intergenic_size) && $min_intergenic_size > 0 && 
				$genes[$gen][0]-$genes[$gen-1][1]+1 < $min_intergenic_size); # skip short intergenes
			next if( defined($max_intergenic_size) && $max_intergenic_size > 0 && 
				$genes[$gen][0]-$genes[$gen-1][1]+1 > $max_intergenic_size); # skip long intergenes
			next if($genes[$gen-1][1]-$genes[$gen-1][0]+1<$length_flanking_ORFs); # skip if neighbor ORFs are too short
			next if($genes[$gen][1]-$genes[$gen][0]+1<$length_flanking_ORFs);
	
			$n_of_intergenic++;

			# calculate intergene boundaries for header (in 1 .. n naturals)
			$start  = $genes[$gen-1][1] + 1 - $length_flanking_ORFs; # first intergene nt - flank
			$end    = $genes[$gen][0] - 1 + $length_flanking_ORFs; # last intergene nt + flank 
			$length = $end - $start + 1;
			
			print FNA ">intergenic$n_of_intergenic|$gbaccession|coords:$start..$end|length:$length|".
				"neighbours:$genes[$gen-1][2]($genes[$gen-1][3]),$genes[$gen][2]($genes[$gen][3])|".
				"neighbour_genes:$genes[$gen-1][4],$genes[$gen][4]|$taxon\n";

			# calculate sequence boundaries for substr, which starts in 0
			#print FNA lc(substr($sequence,$start-1,$length))."\n"; # in one go!
			# left ORF upper case
			print FNA substr($sequence,$genes[$gen-1][1]-$length_flanking_ORFs,$length_flanking_ORFs);
			# intergene lower case
			print FNA lc(substr($sequence,$genes[$gen-1][1],$genes[$gen][0]-$genes[$gen-1][1]-1));
			# right ORF upper case
			print FNA substr($sequence,$genes[$gen][0]-1,$length_flanking_ORFs)."\n";
		}
	}

	close(FNA);

	# make final outfiles only when process finishes cleanly
	rename( "$out_intergenic_file.tmp" , $out_intergenic_file );

	return $n_of_intergenic;
}

sub extract_features_from_genbank
{
#     gene            230511..230630
#                     /locus_tag="ECSE_5SrRNA01"
#     rRNA            230511..230630
#                     /locus_tag="ECSE_5SrRNA01"
#                     /product="5S ribosomal RNA"
#     gene            230683..230759
#                     /locus_tag="ECSE_tRNA03"
#     tRNA            230683..230759
#                     /locus_tag="ECSE_tRNA03"
#                     /product="tRNA-Asp"
#     gene            230922..231725
#                     /locus_tag="ECSE_0203"
#     CDS             230922..231725
#                     /locus_tag="ECSE_0203"
#                     /product="2,5-diketo-D-gluconate reductase"

	# takes a genbank input file and creates a single FASTA nucleotide file containing all features 
        # defined in argument $features_to_parse or $ENV{'GBKFEATURES'}
	# if $out_feature_file == 0 does not create any outfile
        # returns reference to features parsed

	my ($infile,$out_feature_file,$features_to_parse) = @_;
    
	if(!$features_to_parse){ $features_to_parse = $ENV{'GBKFEATURES'} }

	my ($gi,$genename,$coords,$gbaccession,$header,$genelength,%already_seen);
	my ($start,$end,$strand,$source);
	my ($taxon,$strain)=('','');
	my $in = new Bio::SeqIO(-file => $infile, -format => 'genbank');

    	# first create tmp outfile, in case sub does not finish cleanly
	if($out_feature_file)
	{
    		open(FNA,">$out_feature_file.tmp") || 
			die "# extract_features_from_genbank : cannot create $out_feature_file.tmp\n";
	}
	while( my $seq = $in->next_seq) # scan all sequences found in $input
    	{
        	$gbaccession = $seq->accession();
        	$taxon = $coords = $genelength = $source = '';
	
        	for my $f ($seq->get_SeqFeatures)
        	{
			if($f->primary_tag() =~ /source/)
                        {
				$source = $f->end(); 
                                if($f->has_tag('organism'))
                                {
                                    foreach my $element ($f->each_tag_value('organism')){ $taxon .= "[$element],"; }
                                    chop $taxon;
                                }
                        }
		 	elsif($f->has_tag('strain') && $strain eq '')
                        {
                        	$strain = join(',',sort $f->each_tag_value('strain'));
                        }
            		elsif($features_to_parse =~ $f->primary_tag())
            		{
		                $gi = $genename = ''; # compatible con subrutina extract_CDS_from_genbank

 				$coords = $f->spliced_seq();

                		if($f->has_tag('locus_tag') && $gi eq '')
                		{
		                    $gi = "ID:".join(',',sort $f->each_tag_value('locus_tag')); 
                		}
				elsif($f->has_tag('db_xref') && $gi eq '')
                                {
                                    $gi = "ID:".join(',',sort $f->each_tag_value('locus_tag')); 
                                }

		                if($f->has_tag('gene'))
                		{
		                    $genename = join(',',sort $f->each_tag_value('gene'));
				    $genename =~ s/\s+/_/g;
                		}
                		elsif($f->has_tag('product'))
                    		{ 
                        		$genename = join(',',sort $f->each_tag_value('product'));
					$genename =~ s/\s+/_/g;
                		}

		   		next if(!$coords->{'seq'});
	
				$start = $f->start(); # 1..n naturals
                                $end   = $f->end();
                                $strand = $f->location()->strand();

				if(!$gi) # in cases such as nonCDS items in Smelilloti Jul2012
				{
					$gi = "ID:$gbaccession:$start-$end";
				}
		
				if(!$already_seen{$source}{"$gi,$start,$end,$strand"})
                		{
					if($out_feature_file)
					{
						$genelength=length($coords->{'seq'});
                                		$header = $gi." |".$taxon."|".$strain."|".$genename."|".$genelength."|".
		                                        "$gbaccession($source):$start-$end:$strand ".$coords->desc();	
						print FNA ">$header\n$coords->{'seq'}\n" if($out_feature_file);
					}
					$already_seen{$source}{"$gi,$start,$end,$strand"} = 1;
                		}
			}
        	}
	}

	if($out_feature_file)
	{
		close(FNA);

	    	# make final outfiles only when process finishes cleanly
	    	rename( "$out_feature_file.tmp" , $out_feature_file );
	}

	return \%already_seen;
}

sub extract_CDSs_from_genbank
{
	# takes a genbank input file and creates two FASTA files containing all CDSs in
	# aminoacid and dna sequences, respectively
	# returns number of CDSs found

	my ($infile,$out_prot_file,$out_dna_file,$verbose) = @_;
	
	my ($gene,$gi,$crossrefs,$gbaccession,$header,$genelength,$CDScoords);
	my ($start,$end,$strand,$source,$protsequence,%already_seen,%genes);
	my $n_of_CDS = 0;
	my ($taxon,$strain)=('','');
	if($verbose){ print "# extract_CDSs_from_genbank : extracting file $infile...\n" }

	my $in = new Bio::SeqIO(-file => $infile, -format => 'genbank' );

	# parse all CDS sequences and store their neighbors
	while( my $seq = $in->next_seq) # scan all sequences found in $input  
	{    
		$gbaccession = $seq->accession();   
		$source = '';

	       	foreach my $f ($seq->get_SeqFeatures) 
	       	{
			if($f->primary_tag() =~ /source/)
			{
				$source = $f->end();			
				if($f->has_tag('organism') && $taxon eq '')
				{
					foreach my $element ($f->each_tag_value('organism')){ $taxon .= "[$element],"; }
					chop $taxon; 
				}	
				if($f->has_tag('strain') && $strain eq '')
				{
					$strain = join(',',sort $f->each_tag_value('strain'));
				}
			}
			elsif($f->primary_tag() =~ /CDS/)
			{
				$gene=$gi=$crossrefs=$genelength=$protsequence='';
			
				#will fail obtaining the nucleotide sequence with frame shifts	
				#LOCUS       NC_002505            2961149 bp    DNA     circular BCT 24-MAY-2010
				#DEFINITION  Vibrio cholerae O1 biovar El Tor str. N16961 chromosome I, complete genome.
				#ACCESSION   NC_002505
				#VERSION     NC_002505.1  GI:1564003
				#dnaG
				# complement(join(552383..553567,553567..554145))
				$CDScoords = $f->spliced_seq(); 
			
				if($f->has_tag('db_xref'))
				{
					$crossrefs = join(',',sort $f->each_tag_value('db_xref'));
					if($crossrefs =~ /(GI\:\d+)/){ $gi = $1; $crossrefs =~ s/$gi// } 
					next if($crossrefs =~ /PSEUDO:/); # no sabemos si es universal, funciona con Bradyrizobium_ORS278.gb
				}
				
				if($f->has_tag('translation'))
				{
					$protsequence = join('',sort $f->each_tag_value('translation'));
				}
				
				if($gi eq '' && $f->has_tag('locus_tag'))
				{
					$gi = "ID:".join(',',sort $f->each_tag_value('locus_tag')); #print "$gi\n";
				}
				
			
				if($f->has_tag('gene'))
				{
					$gene = join(',',sort $f->each_tag_value('gene'));
				}
				else
				{
					if($f->has_tag('product'))
					{ 
						$gene = join(',',sort $f->each_tag_value('product')); 
					}
				}			
	
				# often a protein sequence is reported in absence of a validated DNA region
				next if($protsequence eq '');

				# NA is assigned when a valid DNA region was not reported
				$genelength= length($CDScoords->{'seq'}) || 'NA';
		
				$start = $f->start(); # 1..n naturals
                                $end   = $f->end();
                                $strand = $f->location()->strand();

				$header = $gi." |".$taxon."|".$strain."|".$gene."|".$genelength."|".
					"$gbaccession($source):$start-$end:$strand" .
					" ^$crossrefs^ ".$CDScoords->desc();     
				
				if(!$already_seen{$header})
				{
					push(@{$genes{$gbaccession}},[$gi,$gene,$strand,$header,$CDScoords->{'seq'},$protsequence]);
					$already_seen{$header} = 1;
					$n_of_CDS++;	
				}
			}
	    	}
	}


	# first create tmp files, in case sub does not finish cleanly
        open(DNA,">$out_dna_file.tmp") || die "# extract_CDSs_from_genbank : cannot create $out_dna_file.tmp\n";
        open(PROT,">$out_prot_file.tmp") || die "# extract_CDSs_from_genbank : cannot create $out_prot_file.tmp\n";	

	foreach $gbaccession (keys(%genes))
	{
		my ($contig_n_of_genes,$contig) = (scalar(@{$genes{$gbaccession}}),$genes{$gbaccession});
		for(my $g=0;$g<scalar(@$contig);$g++)
		{
			$header = "$contig->[$g][3]|";
			if($g==0){ $header .= "neighbours:start(),"; }
			else{ $header .= "neighbours:$contig->[$g-1][0]($contig->[$g-1][2]),"; }
			if($g<$contig_n_of_genes-1)
			{ 
				$header .= "$contig->[$g+1][0]($contig->[$g+1][2])|"; 
			}
			else{ $header .= "end()|" }
			
			if($g==0){ $header .= "neighbour_genes:start(),"; }
                        else{ $header .= "neighbour_genes:$contig->[$g-1][1],"; }
                        if($g<$contig_n_of_genes-1)
                        {
                                $header .= "$contig->[$g+1][1]|"; 
                        }
                        else{ $header .= "end()|" }

			print DNA ">$header\n$contig->[$g][4]\n";
			print PROT ">$header\n$contig->[$g][5]\n";
		}
	}
	close(DNA);
	close(PROT);
      
	# make final outfiles only when process finishes cleanly
        rename( "$out_prot_file.tmp" , $out_prot_file );
	rename( "$out_dna_file.tmp" , $out_dna_file );
	
       	return $n_of_CDS;
}	

# Nov2012 to be used when extract_CDSs fails, in cases such as FN869568.gbk
sub extract_genes_from_genbank
{
	# takes a genbank input file and creates two FASTA files containing all genes in
        # aminoacid and dna sequences, respectively
        # returns number of genes found

        my ($infile,$out_prot_file,$out_dna_file,$verbose) = @_;

        my ($gene,$gi,$gbaccession,$header,$genelength,$CDScoords);
        my ($start,$end,$strand,$source,$protsequence,%already_seen,%genes);
        my $n_of_genes = 0;
        my ($taxon,$strain)=('','');
        if($verbose){ print "# extract_genes_from_genbank : extracting file $infile...\n" }

        my $in = new Bio::SeqIO(-file => $infile, -format => 'genbank' );
	my $seqobj = Bio::Seq->new( -display_id => 'tmp' );

	while( my $seq = $in->next_seq) # scan all sequences found in $input
        {
                $gbaccession = $seq->accession();
                $source = '';
                foreach my $f ($seq->get_SeqFeatures)
                {
                        if($f->primary_tag() =~ /source/)
                        {
                                $source = $f->end();
                                if($f->has_tag('organism') && $taxon eq '')
                                {
                                        foreach my $element ($f->each_tag_value('organism')){ $taxon .= "[$element],"; }
                                        chop $taxon;
                                }
                                if($f->has_tag('strain') && $strain eq '')
                                {
                                        $strain = join(',',sort $f->each_tag_value('strain'));
                                }
                        }
                        elsif($f->primary_tag() =~ /gene/)
                        {
                                $gene=$gi=$genelength=$protsequence='';

                                $CDScoords = $f->spliced_seq();
				$seqobj->seq($CDScoords->{'seq'});
				$protsequence = $seqobj->translate()->seq();chop($protsequence); # faster 
				#$protsequence = translate_dna2prot($CDScoords->{'seq'}); chop($protsequence);
				#print "$CDScoords->{'seq'}\n$protsequence\n";

				# check internal stop codons
                                next if($protsequence =~ /\*/);

                                if($gi eq '' && $f->has_tag('locus_tag'))
                                {
                                        $gi = "ID:".join(',',sort $f->each_tag_value('locus_tag')); #print "$gi\n";
                                }

                                if($f->has_tag('gene'))
                                {
                                        $gene = join(',',sort $f->each_tag_value('gene'));
                                }

                                # NA is assigned when a valid DNA region was not reported
                                $genelength= length($CDScoords->{'seq'}) || 'NA';

                                $start = $f->start(); # 1..n naturals	
				$end   = $f->end();
                                $strand = $f->location()->strand();

                                $header = $gi." |".$taxon."|".$strain."|".$gene."|".$genelength."|".
                                        "$gbaccession($source):$start-$end:$strand" .
                                        " ^^ ".$CDScoords->desc();
                                
                                if(!$already_seen{$header})
                                {
					push(@{$genes{$gbaccession}},[$gi,$gene,$strand,$header,$CDScoords->{'seq'},$protsequence]);
                                        $already_seen{$header} = 1;
                                        $n_of_genes++;    
                                }
			}
		}
	}

	# first create tmp files, in case sub does not finish cleanly
        open(DNA,">$out_dna_file.tmp") || die "# extract_genes_from_genbank : cannot create $out_dna_file.tmp\n";
        open(PROT,">$out_prot_file.tmp") || die "# extract_genes_from_genbank : cannot create $out_prot_file.tmp\n";
	
	foreach $gbaccession (keys(%genes))
        {
                my ($contig_n_of_genes,$contig) = (scalar(@{$genes{$gbaccession}}),$genes{$gbaccession});
                for(my $g=0;$g<scalar(@$contig);$g++)
                {
                        $header = "$contig->[$g][3]|";
                        if($g==0){ $header .= "neighbours:start(),"; }
                        else{ $header .= "neighbours:$contig->[$g-1][0]($contig->[$g-1][2]),"; }
                        if($g<$contig_n_of_genes-1)
                        {
                                $header .= "$contig->[$g+1][0]($contig->[$g+1][2])|";
                        }
                        else{ $header .= "end()|" }

			if($contig->[$g-1][1] || $contig->[$g+1][1])
			{
	                        if($g==0){ $header .= "neighbour_genes:start(),"; }
	                        else{ $header .= "neighbour_genes:$contig->[$g-1][1],"; }
	                        if($g<$contig_n_of_genes-1)
	                        {
	                                $header .= "$contig->[$g+1][1]|";
	                        }
	                        else{ $header .= "end()|" }
			}

                        print DNA ">$header\n$contig->[$g][4]\n";
                        print PROT ">$header\n$contig->[$g][5]\n";
                }
        }
	close(DNA);
        close(PROT);

        # make final outfiles only when process finishes cleanly
        rename( "$out_prot_file.tmp" , $out_prot_file );
        rename( "$out_dna_file.tmp" , $out_dna_file );

        return $n_of_genes;
}

sub extract_gene_name
{
	# returns gene name from a FASTA header generated by extract_CDSs_from_genbank
	return (split(/\|/,$_[0]))[3];
}

# sub DB_row_Puzzle_LM_DNA
# {
# 
#    my (%LM_data) = @_;
#    my  $DB_row = $LM_data{$Alignment_ID} ."\t" . $LM_data{'n_of_sequences'} ."\t" . $LM_data{'nt_sites'} ."\t" .
#           $LM_data{'percent_constant_sites'} ."\t" . $LM_data{'n_of_site_patterns'} ."\t" .
# 	  $LM_data{'percent_constant_site_patterns'} ."\t" . $LM_data{'percent_gaps'} ."\t" . 
# 	  $LM_data{'ChiSQ_failed'} ."\t" . $LM_data{'Acc_No_ChiSQ_perc'} ."\t" . $LM_data{'percent_fullyresolved_quartets'} ."\t" .
# 	  $LM_data{'percent_partiallyresolved_quartets'} ."\t" . $LM_data{'percent_nonresolved_quartets'} . "\n"; 
#    
#      
#    return $DB_row;
# 
# }





# sub upload_LM_data2DB
# {
# 	my (%LM_puzzle_DNA_data) = @_;
# 
# 	
# 	
# 	
# 	
# 	my $key;
# 	my $DB_row;
# 
# 	foreach $key (sort keys %LM_puzzle_DNA_data)
# 	{
# 	  $DB_row = $LM_data{$Alignment_ID} ."\t" . $LM_data{'n_of_sequences'} ."\t" . $LM_data{'nt_sites'} ."\t" .
#           $LM_data{'percent_constant_sites'} ."\t" . $LM_data{'n_of_site_patterns'} ."\t" .
# 	  $LM_data{'percent_constant_site_patterns'} ."\t" . $LM_data{'percent_gaps'} ."\t" . 
# 	  $LM_data{'ChiSQ_failed'} ."\t" . $LM_data{'Acc_No_ChiSQ_perc'} ."\t" . $LM_data{'percent_fullyresolved_quartets'} ."\t" .
# 	  $LM_data{'percent_partiallyresolved_quartets'} ."\t" . $LM_data{'percent_nonresolved_quartets'} ."\n";
# 	}
# 
# 	# 1) concatena separando por comas los valroes del hash en el orden de la tabla que llenes,
# 	# las cadenas entre comillas simples
# 	
# 	$DB_row = $LM_data{$Alignment_ID} ."\t" . $LM_data{'n_of_sequences'} ."\t" . $LM_data{'nt_sites'} ."\t" .
#           $LM_data{'percent_constant_sites'} ."\t" . $LM_data{'n_of_site_patterns'} ."\t" .
# 	  $LM_data{'percent_constant_site_patterns'} ."\t" . $LM_data{'percent_gaps'} ."\t" . 
# 	  $LM_data{'ChiSQ_failed'} ."\t" . $LM_data{'Acc_No_ChiSQ_perc'} ."\t" . $LM_data{'percent_fullyresolved_quartets'} ."\t" .
# 	  $LM_data{'percent_partiallyresolved_quartets'} ."\t" . $LM_data{'percent_nonresolved_quartets'};
# 	
# 	
# 	
	
	# 2) inserta la concatenaci? resultante en VALUES($concatenaci?)


#my $conexion = DBI->connect("dbi:mysql:PhyMarkers:servidor.ccg.unam.mx","usuario","clave",{PrintError,0,RaiseError,1,AutoCommit,0});
#
## filas contiene el nmero de filas afectadas por la sentencia,
## en caso de error devuelve falso
#my $filas = $conexion->do("INSERT INTO variables VALUES ('grafo', 0)");
#if(!$filas)
#{
#	print "no pude ejecutar la expresi? SQL\n";
#}
#else
#{
#	# si quieres fijar estos cambios:
#	$conexion->commit(); 
#	
#	# si quisieras descartalos:
#	$conexion->rollback();
#}
#
#$dbh->rollback();
#
#$conexion->disconnect();
# 
# 	return 1;
# 	
# 	#if error return  0;
# 
# }

#sub chisquaretable = (
#50,{0.05,67.51,0.01,76.15,0.001,86.66},
#51,{0.05,68.67,0.01,77.39,0.001,87.97},
#52,{0.05,69.83,0.01,78.62,0.001,89.27},
#53,{0.05,70.99,0.01,79.84,0.001,90.57},
#54,{0.05,72.15,0.01,81.07,0.001,91.88},
#55,{0.05,73.31,0.01,82.29,0.001,93.17},
#56,{0.05,74.47,0.01,83.52,0.001,94.47},
#57,{0.05,75.62,0.01,84.73,0.001,95.75},
#58,{0.05,76.78,0.01,85.95,0.001,97.03},
#59,{0.05,77.93,0.01,87.17,0.001,98.34},
#60,{0.05,79.08,0.01,88.38,0.001,99.62},
#61,{0.05,80.23,0.01,89.59,0.001,100.88},
#62,{0.05,81.38,0.01,90.80,0.001,102.15},
#63,{0.05,82.53,0.01,92.01,0.001,103.46});
#}


1;



