run_phylip.sh

#!/usr/bin/env bash

#-------------------------------------------------------------------------------------------------------
#: PROGRAM: run_phylip.sh
#: AUTHOR: Pablo Vinuesa, Center for Genomic Sciences, UNAM, Mexico
#:         https://www.ccg.unam.mx/~vinuesa/ twitter: @pvinmex
#
#: PROJECT START: October 16th, 2013
#    This program has been developed mainly for teaching purposes, 
#    with improvements/new features added as the script was used in 
#    diverse courses taught to undergrads at https://www.lcg.unam.mx
#    (LCG-UNAM) and the International Workshops on Bioinformatics (TIB)

#: AIM: run PHYLIP's distance methods [NJ|UPGMA] for DNA and proteins (dnadist|protdist) with optional bootstrapping
#       This script was written to teach intermediate Bash scripting to my students at the 
#       Bachelor's Program in Genome Sciences at the Center for Genome Sciences, UNAM, Mexico
#       https://www.lcg.unam.mx
#
#: INPUT: multiple sequence alignments (PROT|DNA) with at least 4 distinct sequences in phylip format
#
#: OUTPUT: [NJ|UPGMA] phylogenies and, if requested, bootstrap consensus trees

#: SOURCE: Freely available on GitHub @ https://github.com/vinuesa/intro2linux
#          Released under the GPLv3 License. 
#          http://www.gnu.org/copyleft/gpl.html

#### DEPENDENCIES 
# Assumes that the following binaries and scripts are all in $PATH, checked by check_dependencies()
#
#   1) Binaries from the PHYLIP package: 
#	seqboot dnadist protdist neighbor consense
#   NOTE: Linux-compatible PHYLIP binaries are supplied in the distro\'s bin/ directory
#           https://github.com/vinuesa/intro2linux/tree/master/bin  

#: TODO:
#    implement also parsimony and ML analyses and parallelize bootstrapping

#: KNOWN BUGS: None.
#    Please report any errors you may encounter through the GitHub issue pages
#-------------------------------------------------------------------------------------------------------

# make sure the user has at least bash version 4, since the script uses standard arrays (introduced in version 4),
#  but future development may require hashes, introduced in version 4
[ "${BASH_VERSION%%.*}" -lt 4 ] && echo "$HOSTNAME is running an ancient bash: ${BASH_VERSION}; ${0##*/} requires bash version 4 or higher" && exit 1

# 0. Define strict bash settings
set -e              # exit on non-zero exit status
set -u              # exit if unset variables are encountered
set -o pipefail     # exit after unsuccessful UNIX pipe command

# set and export LC_NUMERIC=en_US.UTF-8, to avoid problems with locales tha use 1,32
LC_NUMERIC=en_US.UTF-8
export LC_NUMERIC

args="$*"
progname=${0##*/} # run_phylip.sh
VERSION=2.0 

# GLOBALS
#DATEFORMAT_SHORT="%d%b%y" # 16Oct13
#TIMESTAMP_SHORT=$(date +${DATEFORMAT_SHORT})

date_F=$(date +%F |sed 's/-/_/g')-   # 2013_10_20
date_T=$(date +%T |sed 's/:/./g')    # 23.28.22 (hr.min.secs)
start_time="$date_F$date_T"
#current_year=$(date +%Y)

wkdir=$(pwd)

# initialize variables
def_DNA_model=F84
def_prot_model=JTT
input_phylip=
runmode=
boot=100
CV=1
model=
DEBUG=0
sequential=0
TiTv=2
upgma=0
outgroup=0
gamma="0.0"
outgroup=1

nw_utils_ok=

declare -a outfiles # to collect the output files written to disk

#---------------------------------------------------------------------------------#
#>>>>>>>>>>>>>>>>>>>>>>>>>>>> FUNCTION DEFINITIONS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<#
#---------------------------------------------------------------------------------#
function print_dev_history()
{
    cat <<EOF_
    
    $progname $VERSION has been developed mainly for teaching purposes, 
      with improvements/new features added as the script was used in 
      diverse courses taught to undergrads at https://www.lcg.unam.mx
      and the International Workshops on Bioinformatics (TIB)
      
    TODO: * implement also parsimony and ML analyses
          * parallelize bootstrapping  
    
    # v2.0 2020-11-19; * fixed bug in write_protdist_params in the if [ $model == JTT ] || [ ...] conditional
   		       * fixed an error in the model name input checks, changing PMG for PMB
   		       * fixed grep in check_nw_utils_ok and made more robust conditional checking to call the function
   		       * use variable expansion only once in naming outfiles outtress, saving it to a variable
   			 to make the code more streamlined and possibly a tiny bit faster
   		       * removed useless [ -s outtre|outfile ] tests for the orignal tress|outfiles in the bootstrapping block, 
   			   as these are already checked in the first code block computing tress from the original alignment
    
    # v1.6 2020-11-19; * streamlined write_PHYLIP_param functions by grouping echo calls in { } > params; avoiding concatenation 
    #                  * improved description of write_PHYLIP_param functionality

    # v1.5 2020-11-15; * added check_phylip_ok to validate input phylip file with -ge 4 sequences and -ge 4 characters
    		       * added remove_phylip_param_files
    		       * added option -I to call print_install_notes
    		       * added [ "${BASH_VERSION%%.*}" -lt 4 ] && die
    		       * makes more extensive use of internal Bash variables and variable expansions, including \${var/pattern/string}

    # v1.4 2020-11-14; * added function check_nw_utils_ok, to check that nw_display and nw_support can be executed
    			  as they may be in path but cannot find /lib64/libc.so.6: version GLIBC_2.14
    			  when the binaries are compiled with dynamic linking to the libs
    
    # v1.3 2020-11-14; * improved layout of output messages; 
    		       * improved regex in extract_tree_from_outfile (now also for NJ tree)
    		       * if nw_support and nw_display are not available, 
    			 prints both the NJ and boot trees to screen if bootstrap analysis was performed
    
    # v1.2 2020-11-11; set and export LC_NUMERIC=en_US.UTF-8, to avoid problems with locales that use 1,32 instead of 1.32
    
    # v1.1  2020-11-11; added function extract_tree_from_oufile and calls it on NJ|UPGMA bootRepl consensus trees
                          if nw_display is not available in PATH
    
    #v1.0  2020-11-09; This version has thorough error checking and reporting with improved code flow
         *  re-ingeneered the main code block. It now first runs the standard distance-matrix + clustering computations
	        before doing bootstrapping. This helps in rapidly detecting problematic model parameters, as reflected in negative distances
	 *  displays NJ|UPGMA tree without bootstrap values, if no bootstrapping is requested	
	 *  added functions:
	      - check_matrix to check if they contain negative values, issuing an error message
	      - display_treeOI display trees with nw_display, only if they do not contain negative branch lengths, issuing an error message
	      - print_dev_history to clean up the script\'s header section
	 *  made 100% shellcheck-compliant; prints all error messages to >&2 (STDERR)
         *  changed/fixed defaults for CV=1 & gamma="0.0"; 
	 *  use outfiles+=() to capture outfiles for -R 1 and -R 2 when boot -eq 0; now runs with -i -R1|2 as single opts
	 *  Changed default boot=100
	 *  now accepts TiTv=real like 6.7 and checks for real numbers passed to -t
	 *  checks if outtree contains negative branches before attempting to display with nw_display & prints proper message

    # v0.5 2020-11-09; added functions:
    	 *  rdm_odd_int to compute random odd integers to pass to PHYLIP programs
    	 *  check_output (silently checks for expected output files, dying if not found)

    #v0.4 2020-11-08; released to the public domain @ https://github.com/vinuesa/intro2linux
    	   * added nw_support and nw_display calls to map and display bottstrap support values onto NJ or UPGMA trees
    	   * fixed several warnings issued by shellcheck
    	   * added strict bash interpreter calling with set -e; set -u; set -o pipefail
    	   * localized variables received in functions
    	   * explicitly set default_DNA_model=F84 and default_prot_model=JTT
    	   * added option -v to print version and exit
    		
    #v0.3 September 15, 2017; improved checking of user input
    #v0.2 September 16, 2015; basic options and parameter checking for protdist phylogenies
    #v0.1 October 16, 2013; basic options and parameter checking for dnadist phylogenies
EOF_

exit 0
}
#---------------------------------------------------------------------------------#

function print_install_notes()
{
   cat <<INSTALL
   
    #1. PHYLIP PHYLogeny Inference Package by Joe Felensetein https://evolution.genetics.washington.edu/phylip.html
    
    The bin/ directory of the intro2linux distribution provides precompiled Linux-x86_64 binaries 
          for seqboot dnadist protdist neighbor consense. Copy them into your \$HOME/bin directory or any other in your PATH
    
    If you need executables for a different architecture, visit https://evolution.genetics.washington.edu/phylip/executables.html
    
    #2. Newick Utilities http://cegg.unige.ch/newick_utils
    The Newick Utilities have been described in an Open-Access paper (that is, available online for everyone):

    The Newick Utilities: High-throughput Phylogenetic tree Processing in the UNIX Shell
    Thomas Junier and Evgeny M. Zdobnov
    Bioinformatics 2010 26:1669-1670
    http://bioinformatics.oxfordjournals.org/content/26/13/1669.full
    doi:10.1093/bioinformatics/btq243
    
    The src/ dir contains the newick-utils-1.6-Linux-x86_64-disabled-extra.tar.gz source file
    
    Visit http://cegg.unige.ch/newick_utils if you need other pre-compiled versions
    
    Brief compilation notes
    
    1. Copy the newick-utils-1.6-Linux-x86_64-disabled-extra.tar.gz to a suitable directory in your \$HOME
    2. cd into the directory hoding the source *tar.gz file
    3. unpack and compile with the following commands 
          tar -xvzf newick-utils-1.6-Linux-x86_64-disabled-extra.tar.gz
	  cd newick-utils-1.6
	  ./configure --prefix=$HOME  # if you do NOT have administrator privileges
          #./configure                # if you have superuser privileges
	  make
	  make check
	  make install                # if you do NOT have administrator privileges
	  # sudo make install         # if you have superuser privileges

INSTALL

exit 0

}
#---------------------------------------------------------------------------------#

function check_output()
{
    local outfile=$1
    if [ ! -s "$outfile" ]
    then
        echo
	echo " >>> ERROR! The expected output file $outfile was not produced, will exit now!" >&2
        echo
	exit 2
    fi
}
#---------------------------------------------------------------------------------#

function check_phylip_ok()
{
    # implements a rudimentary format validation check
    local phyfile=$1
    local num_tax=
    local num_char=

    num_tax=$(awk 'NR == 1 {print $1}' "$phyfile")
    num_char=$(awk 'NR == 1 {print $2}' "$phyfile")
    
    if [[ ! "$num_tax" =~ ^([0-9]+)$ ]] ||  [[ ! "$num_char" =~ ^([0-9]+)$ ]]
    then
         echo "ERROR: $phyfile is not a phylip-formatted input file!"
	 echo
	 exit 1  
    elif [[ "$num_tax" =~ ^([0-9]+)$ ]] && [ "$num_tax" -lt 4 ]
    then
         echo "ERROR: $phyfile should contain at least 4 taxa"
	 exit 1
    elif [[ "$num_char" =~ ^([0-9]+)$ ]] && [ "$num_char" -lt 5 ]
    then
         echo "ERROR: $phyfile should contain at least 5 aligned residues"
	 exit 1
    else
        echo "# File validation OK: $phyfile seems to be a standard phylip file ..."
        echo '-------------------------------------------------------------------------------------------'
    fi
}
#---------------------------------------------------------------------------------#

function check_dependencies
{
    # check that the following PHYLIP binaries are in $PATH; die if not
    # optional: the Newick_utilities src is found in the src/ dir
    dependencies=(seqboot dnadist protdist neighbor consense) 
    for programname in "${dependencies[@]}"
    do
       local bin=
       bin=$(type -P "$programname")
       if [ -z "$bin" ]; then
          echo
          echo "# ERROR: $programname not in place!" >&2
          echo "# ... you will need to install \"$programname\" first or include it in \$PATH" >&2
          echo "# ... exiting" >&2
          exit 1
       fi
    done

    echo
    echo '# Run check_dependencies() ... looks good: all required external dependencies are in place!'
    echo '-------------------------------------------------------------------------------------------'
}
#---------------------------------------------------------------------------------#

function check_nw_utils_ok()
{
    # This function checks that nw_display and nw_support can be executed
    #  as they may be in path but cannot find /lib64/libc.so.6: version GLIBC_2.14
    #  when the binaries are compiled with dynamic linking to the libs
    
    # The function returns 0 when nw_support nw_display cannot be executed and 1 when they can
    # The function does not run and return anything if nw_support nw_display are not in PATH
    
    dependencies=(nw_support nw_display)
    
    local nw_utils_ok=
    local bin=
    
    for programname in "${dependencies[@]}"
    do
        bin=$(type -P "$programname")
        if [ -n "$bin" ]; then
	   # prints something like nw_support: /lib64/libc.so.6: version `GLIBC_2.14' not found
	   "$programname" 2>> "nw_utils_check.out.tmp" 
        fi
    done
    
    # check the contents of nw_utils_check.out.tmp and set proper flag
    if [ -s nw_utils_check.out.tmp ]
    then
           # set the nw_utils_ok flag to 0 if the lib is not found, to check later in the code when nw_support and nw_display are called
           grep -i error nw_utils_check.out.tmp &> /dev/null
	   nw_utils_ok=$? 
	   [ "$nw_utils_ok" -eq 0 ] && echo "# WARNING: ${dependencies[*]} are in PATH but cannot find /lib64/libc.so.6: version GLIBC_2.14" >&2 
           rm nw_utils_check.out.tmp 
           echo  "$nw_utils_ok"
    fi
}
#---------------------------------------------------------------------------------#

function rdm_odd_int()
{
  # generate a random odd integer for seqboot and neighbor
  while true
  do 
     i=$RANDOM
     if (( "$i" % 2 )) # true when odd
     then 
         echo "$i" 
	 break
     fi
done

}
#---------------------------------------------------------------------------------#

function check_matrix()
{
    # check that the distance matrix does not contain negative values due to wrong model parameters
    local matrix=$1
    
    if grep -El ' \-[0-9\.]+' "$matrix"
    then
        echo "ERROR: computed negative distances!" >&2
	     [ "$runmode" -eq 1 ] && echo "You may need to ajust the model and|or gamma value; try lowering TiTv if > 6" >&2
	     [ "$runmode" -eq 2 ] && echo "You may need to ajust the matrix and|or gamma value" >&2
	exit 3
    fi
}
#---------------------------------------------------------------------------------#

function display_treeOK()
{
     # checks that tree does not contain branches with negative lengths, 
     #    as these cannot be displayed  with nw_dsiplay
     local tree=$1
     
     # check that there are no negative branch lengths in the nj_tree
     if ! grep -El '\-[0-9\.]+' "$tree"
     then 
     	 echo "# displaying $tree ..."
     	 echo
     	 nw_display -w 80 -Ir "$tree"
     	 echo
     else
     	 echo "ERROR: cannot display $nj_tree with nw_dsiplay, as it contains branches with negative lengths!" >&2
	     [ "$runmode" -eq 1 ] && echo "You may need to adjust (lower?) TiTv=$TiTv" >&2
	     [ "$runmode" -eq 2 ] && echo "You may need to use another matrix or adjunst gamma" >&2
     	 exit 4
     fi
}
#---------------------------------------------------------------------------------#

function extract_tree_from_outfile()
{
    # grep out lines containing tree characters | - and print to STDOUT
    local outfile=$1
   
    grep --color=never -E '[[:blank:]]+\||-' "$outfile" | grep -Ev 'Neighbor-|^---'
}

#------------------------------------- PHYLIP FUNCTIONS ---------------------------------------#
# >>> these are fucntions to write the command files to pass parameters to PHYLIP programs <<< # 

function remove_phylip_param_files()
{
   [ -s dnadist.params ] && rm dnadist.params
   [ -s protdist.params ] && rm protdist.params
   [ -s seqboot.params ] && rm seqboot.params
   [ -s neighbor.params ] && rm neighbor.params
   [ -s consense.params ] && rm consense.params
   return 0
}
#------------------------------------- PHYLIP FUNCTIONS ---------------------------------------#

function write_dnadist_params
{
    # writes a parameter file to run dnadist, based on provided arguments
    local model=$1
    local boot=$2
    local TiTv=$3
    local sequential=$4
    local gamma=$5
    local CV=$6

    # Runmode 1 = dnadist 
    if [ "$model" = 'F84' ] || [ -z "$model" ]
    then 
        {
          [ "$model" = 'Kimura' ]       && echo "D"                 
          [ "$model" = 'Jukes-Cantor' ] && echo -ne "D\nD\n"        
          [ "$model" = 'LogDet' ]       && echo -ne "D\nD\nD\n"     
          [ "$TiTv" != "2" ]            && echo -ne "T\n$TiTv\n"    
          [ "$gamma" != "0" ]           && echo -ne "G\n"           
          [ "$boot" -gt 0 ]             && echo -ne "M\nD\n$boot\n" 
          [ "$sequential" -eq 1 ]       && echo -ne "I\n"           
                                           echo -ne "Y\n"             
          [ "$gamma" != "0" ]           && echo -ne "$CV\n"         
    
        } > dnadist.params
   fi
}
#---------------------------------------------------------------------------------#

function write_protdist_params
{
    # write_protdist_params "$model" "$boot" "$sequential" "$gamma" "$CV"
    # writes a parameter file to run protdist, based on provided arguments
    # Models: JTT, PMB, PAM, Kimura
    local model=$1
    local boot=$2
    local sequential=$3
    local gamma=$4
    local CV=$5

    if [ "$model" = 'JTT' ] || [ -n "$model" ]
    then
       { # Runmode 2 = protdist 
         [ "$model" = 'PMB' ]     && echo "P"                 
         [ "$model" = 'PAM' ]     && echo -ne "P\nP\n"        
         [ "$model" = 'Kimura' ]  && echo -ne "P\nP\nP\n"     
         [ "$gamma" != "0" ]      && echo -ne "G\n"           
         [ "$boot" -gt 0 ]        && echo -ne "M\nD\n$boot\n" 
         [ "$sequential" -eq 1 ]  && echo -ne "I\n"    
                                     echo -ne "Y\n"       
         [ "$gamma" != "0" ]      && echo -ne "$CV\n"   			   
       } > protdist.params
    fi  
}
#---------------------------------------------------------------------------------#

function write_seqboot_params
{
    # writes a parameter file to run seqboot, based on provided arguments
    local boot=$1
    local sequential=$2
    
    # Write Seqboot params
    {
      [ "$boot" -gt 0 ]       && echo -ne "R\n$boot\n" 
      [ "$sequential" -eq 1 ] && echo -ne "I\n"        
                                 echo -ne "Y\n"
                                 echo -ne "$ROI\n"
    } > seqboot.params 			      
}
#---------------------------------------------------------------------------------#

function write_neighbor_params
{
    # writes a parameter file to run neighbor, based on provided arguments
    local boot=$1
    local upgma=$2
    local outgroup=$3
    
    {
       # Write $ROI params
       [ "$upgma" -gt 0 ]    && echo -ne "N\n"                 
       [ "$outgroup" -gt 1 ] && echo -ne "O\n$outgroup\n"      
       [ "$boot" -gt 0 ]     && echo -ne "M\n$boot\n$ROI\n"    
                                echo -ne "Y\n"        
    } > neighbor.params
}
#---------------------------------------------------------------------------------#

function write_consense_params
{
    #writes a parameter file to run consense; very difficult ;)
    [ -s consense.params ] && rm consense.params
    {
      [ "$outgroup" -gt 1 ]  && echo -ne "O\n$outgroup\n"  
                                echo -ne "Y\n"           
    } > consense.params				
}
#---------------------------------------------------------------------------------#

# don't leave litter behind ... remove intermediate input/output files
function cleanup_dir 
{
   [ -s infile ]  && rm infile
   [ -s outfile ] && rm outfile
   [ -s intree ]  && rm intree
   [ -s outtree ] && rm outtree
   
   for file in *.params
   do
       [ -s "$file" ] && rm "$file"
   done    
}

#-------------------------------- END PHYLIP FUNCTIONS----------------------------#

function print_help
{
   #':b:c:m:I:o:R:t:hHIDsuv'
   cat<<EOF
   $progname v.$VERSION [OPTIONS]
   
   -h  prints this help message
   
 REQUIRED
   -i <input_phylip_file> (if sequential, use also flag -s)
   -R <integer>  RUNMODE
        1    run dnadist  (expects DNA  alignment in phylip format)
	2    run protdist (expects PROT alignment in phylip format)

 OPTIONAL
   -b <integer> No. of bootstrap pseudoreplicates [default: $boot]
   -g <real number> alpha for gamma distribution  [default: $gamma]
   -o <integer> outgroup sequence no.             [default: $outgroup]
   -s sequential format; -s (flag; no val)!       [default: interleaved]
   -t <digit> Transition/Transversion ratio       [default: $TiTv]
   -m <model name>  NOTE: TYPE AS SHOWN!
        DNA:  F84, Kimura, Jukes-Cantor, LogDet   [default: $def_DNA_model]
        PROT: JTT, PMB, PAM, Kimura               [default: $def_prot_model]
   -u <flag> use UPGMA for clustering             [default:  NJ]
   -v <flag> print program version and exit
   -D <flag> Activate debugging to keep cmd files [default: $DEBUG]
   -H <flag> print development history and exit   
   -I <flag> print installation notes and exit   
   
 AIM: run PHYLIP\'s distance methods [NJ|UPGMA] for DNA and proteins with bootstrapping
      This code was written to teach basic Bash scripting to my students at the 
      Bachelor\'s Program in Genome Sciences, Center for Genome Sciences, UNAM, Mexico
      https://www.lcg.unam.mx/      

 OUTPUT: [NJ|UPGMA] phylogenies and, if requested, bootstrap consensus trees 
         and  [NJ|UPGMA] phylogenies with bootstrap support values mappend on bipartitions
 
 EXTERNAL DEPENDENCIES: 
   * PHYLIP (https://evolution.genetics.washington.edu/phylip.html) programs:
      seqboot dnadist protdist neighbor consense 
   * Newick utilities programs (http://cegg.unige.ch/newick_utils) programs:
      optional: nw_support and nw_display; need to install separately
   - Notes: PHYLIP Linux 64bit binaries available in the bin/ dir
            Copy them to your \$HOME/bin directory or any other in your PATH
  
 LICENSING & SOURCE
       Author: Pablo Vinuesa | https://www.ccg.unam.mx/~vinuesa/ | twitter: @pvinmex
       Released under the GNU General Public License version 3 (GPLv3)
       http://www.gnu.org/copyleft/gpl.html
       source: https://github.com/vinuesa/intro2linux

EOF

exit 1   
}

#-------------------------------------------------------------------------------------------------#
#---------------------------------------- GET OPTIONS --------------------------------------------#
#-------------------------------------------------------------------------------------------------#

# GETOPTS
while getopts ':b:g:i:m:R:o:t:hDHIsuv' OPTIONS
do
   case $OPTIONS in

   b)   boot=$OPTARG
        ;;
   g)   gamma=$OPTARG
	[ "$gamma" != 0 ] && CV=$(echo "1/sqrt($gamma)" | bc -l) && printf -v CV "%.4f\n" "$CV"
        ;;
   h)   print_help
        ;;
   i)   input_phylip=$OPTARG
        ;;
   s)   sequential=1
        ;;
   m)   model=$OPTARG
        ;;
   o)   outgroup=$OPTARG
        ;;
   R)   runmode=$OPTARG
        ;;
   t)   TiTv=$OPTARG
        ;;
   u)   upgma=1
        ;;
   v)   echo "$progname v$VERSION" && exit
        ;;
   D)   DEBUG=1
        ;;
   H)   print_dev_history
        ;;
   I)   print_install_notes
        ;;
   :)   printf "argument missing from -%s option\n" "$OPTARG"
   	print_help
     	exit 2 
     	;;
   \?)  echo "need the following args: "
   	print_help
        exit 3
	;;
   esac >&2   # print the ERROR MESSAGES to STDERR

done

shift $((OPTIND - 1))

#--------------------------#
# >>> Check User Input <<< #
#--------------------------#

if [ -z "$input_phylip" ]
then
       echo
       echo "# ERROR: no input phylip file defined!"
       print_help
       exit 1    
fi

if [ -z "$runmode" ]
then
       echo
       echo "# ERROR: no runmode defined!"
       print_help
       exit 1    
fi

# automatically set TiTv=0 when running with protein matrices or the LogDet DNA model
[ "$runmode" -eq 2 ] && TiTv=0
[ "$model" == LogDet ] && TiTv=0


# check that bootstrap value is an integer
re='^[0-9]+$'
if [[ ! "$boot" =~ $re ]]
then
   echo
   echo "# ERROR: boot:$boot is not a positive integer >= 0; provide a value between 100 and 1000" >&2 
   echo
   print_help
   echo
   exit 3
fi

# check that Ti/Tv is a real number, integer or decimal
re_TiTv='^[0-9.]+$'
if [[ ! "$TiTv" =~ $re_TiTv ]] && [ "$runmode" -eq 1 ]
then
   echo
   echo "# ERROR: TiTv:$TiTv is not an integer >= 0; provide a value between 0-10" >&2 
   echo
   print_help
   echo
   exit 3
elif [[ "$TiTv" =~ $re_TiTv ]] && [ "$runmode" -eq 1 ] && { [ "$model" == Jukes-Cantor ] || [ "$model" == LogDet ]; }
then
       echo
       echo "# ERROR: $model is only valid when analyzing DNA alignments under the F84 or K2P models" >&2
       echo
       print_help
       exit 1    
fi

if [ "$gamma" != 0 ] # note, need to treat as string, since gamma will be generalle a float like 0.5
then
    # this is to avoid having dots within filename (converts gamma==0.2 to gammaf=02)
    gammaf=${gamma/\./} #gammaf=${gamma%.*}${gamma##*.}
else
    gammaf="$gamma"
fi    

# check model vs. runmode compatibility and suitable bootstrap values are provided
#   using two alternative sintaxes for educational purposes
if [ "$runmode" -eq 1 ] && { [ "$model" = JTT ] || [ "$model" = PMB ] || [ "$model" = PAM ] || [ "$model" = "$def_prot_model" ]; }
then
       echo
       echo "# ERROR: $model is only valid when analyzing protein alignments under -R 2" >&2 
       echo
       print_help
       exit 1    
elif [ "$runmode" -eq 2 ] && [[ "$model" =~ ^(F84|Jukes-Cantor|LogDet|"$def_DNA_model")$ ]]
then
       echo
       echo "# ERROR: $model is only valid when analyzing DNA alignments under -R 1" >&2
       echo
       print_help
       exit 1    
elif [ "$boot" -lt 0 ]
then
       echo
       echo "# ERROR: bootstrap value=$boot is < 0 and not permitted" >&2
       echo
       print_help
       exit 1    
elif [ "$boot" -gt 1000 ]
then
       echo
       echo "# WARNING: bootstrap value=$boot is > 1000. This may take a long time to run."
       echo -n "# Are you sure you want to continue anyway? Type Y|N "
       read -r answer
       if [ "$answer" = N ] || [ "$answer" = n ] 
       then
           echo
	   echo "# Ok, will exit then!"
	   echo
	   exit 2
       else
           echo
	   echo "# Ok, will proceed running with $boot bootstrap pseudoreplicates ..."
	   echo
       fi 
elif [ "$boot" -gt 0 ] && [ "$boot" -lt 100 ]
then
       echo
       echo "# WARNING: bootstrap value=$boot is a bit low. Use a value >= 100."
       echo -n "# Are you sure you want to continue anyway? Type Y|N "
       read -r answer
       if [ "$answer" = N ] || [ "$answer" = n ] 
       then
           echo
	   echo "# Ok, will exit then!"
	   echo
	   exit 2
       else
           echo
	   echo "# Ok, will proceed ..."
	   echo
       fi 
fi

# set the default DNA or PROT models, if not defined by the user
if [ "$runmode" -eq 1 ] && [ -z "$model" ]
then
       model="$def_DNA_model"
       #echo "# DNA model set to default: $def_DNA_model..."       
fi

if [ "$runmode" -eq 2 ] && [ -z "$model" ]

then
       model="$def_prot_model"
       #echo "# Protein matrix set to default: $def_prot_model ..."       
fi

# check that the user provided a valid DNA substitution model
if [ "$runmode" -eq 1 ] && [[ ! "$model" =~ ^(F84|Kimura|Jukes-Cantor|LogDet|"$def_DNA_model")$ ]]
then
   echo
   echo "# ERROR: $model is not a recognized substitution model for DNA sequences used by PHYLIP" >&2 
   echo
   print_help
   echo
   exit 3
fi

# check that the user provided a valid name for empirical substitution matrix
#    note the much shorter and cleaner notation of this test using extended regexes, than the previous one
if [ "$runmode" -eq 2 ] && [[ ! "$model" =~ ^(JTT|PMB|PAM|Kimura)$ ]]
then
   echo
   echo "# ERROR: $model is not a recognized substitution matrix for protein sequences used by PHYLIP" >&2 
   echo
   print_help
   echo
   exit 3
fi

#>>> Set the script's run random odd number required by diverse PHYLIP programs <<<
ROI=$(rdm_odd_int)

# print the run settings
echo
echo "### $progname v.$VERSION run on $start_time with the following parameters:"
echo "# work_directory=$wkdir"
echo "# input_phylip=$input_phylip"
echo "# model=$model | gamma=$gamma | CV=$CV | gammaf=$gammaf | ti/tv ratio=$TiTv |
     outgroup=$outgroup | UPGMA=$upgma | bootstrap no.=$boot | ROI=$ROI"
echo "# runmode=$runmode"
echo "# DEBUG=$DEBUG" 
echo "# command: $progname ${args[*]}"
echo

#-------------------------------------------------------------------------------------------------#
#------------------------------------------ MAIN CODE --------------------------------------------#
#-------------------------------------------------------------------------------------------------#

# 1. make sure the external dependencies are found in PATH
#     and that nw_display and nw_support find the required GLIBC_2.14 in /lib64/libc.so.6
check_dependencies
nw_utils_ok=$(check_nw_utils_ok) # nw_utils_ok -eq 1 when OK, -eq 0 when not

# 2. Start processing the input file
# make sure there are no old outfile or outtree files lying around from previous runs
[ -s infile ] && rm infile
[ -s outfile ] && rm outfile
[ -s outtree ] && rm outtree
# nor old params files
remove_phylip_param_files

# 2.1) make sure we have an input file or die
if [ -s "$input_phylip" ]
then
     # make basic format validation check
     check_phylip_ok "$input_phylip"
     cp "$input_phylip" infile
else
     echo
     echo "# FATAL ERROR: input phylip file $input_phylip does not exist or is empty" >&2 
     echo "# Make sure $input_phylip is in $wkdir. Exiting ..." >&2
     echo
     exit 1
fi

# ----------------------------------------------------------- #
# >>>>>>>>>>>>>>>> Compute distance matrices <<<<<<<<<<<<<<<< #
# ----------------------------------------------------------- #

# 3. In any case (with or without bootstrap) compute the NJ tree for the original phylip file
#     Run dnadist or protdist, as required   
echo
echo ">>> Computing distance matrix for $input_phylip ..."

if [ "$runmode" -eq 1 ]		 
then
    echo "# running write_dnadist_params $model 0 $TiTv $sequential $gamma $CV"
    write_dnadist_params "$model" "0" "$TiTv" "$sequential" "$gamma" "$CV"
    echo "# running dnadist < dnadist.params"
    dnadist < dnadist.params &> /dev/null
    check_output outfile
    
    # https://fvue.nl/wiki/Bash:_Error_%60Unbound_variable%27_when_appending_to_empty_array
    # Set last item specifically
    # nstead of appending one element, set the last item specifically, without any "unbound variable" error
    # t[${#t[*]}]=foo
    dnadist_outfile="${input_phylip%.*}_${model}${gammaf}gamma_distMat.out"
    cp outfile "$dnadist_outfile" && \
      outfiles+=("$dnadist_outfile")
    
    # check that the distance matrix does not contain negative values due to wrong model parameters
    check_matrix "$dnadist_outfile"    
    
    mv outfile infile
elif [ "$runmode" -eq 2 ]
then
    echo "# running write_protdist_params $model 0 $sequential $gamma $CV"
    write_protdist_params "$model" "0" "$sequential" "$gamma" "$CV"
    echo "# running protdist < protdist.params"
    protdist < protdist.params &> /dev/null
    check_output outfile
    
    protdist_outfile="${input_phylip%.*}_${model}${gammaf}gamma_distMat.out"
    cp outfile "$protdist_outfile" && \
      outfiles+=("$protdist_outfile")
    
    # check that the distance matrix does not contain negative values due to wrong model parameters
    check_matrix "$protdist_outfile"
    
    mv outfile infile
fi


# ------------------------------------------------------------------------------------ #
# >>>>>>>>>>>>>>>> Computing NJ|UPGMA trees from original alignments <<<<<<<<<<<<<<<<< #
# ------------------------------------------------------------------------------------ #

echo
echo ">>> Computing distance tree for $input_phylip ..."  

# 4. now that we have the dist matrix, do the clustering with NJ or UPGMA
echo "# running write_neighbor_params 0 $upgma"
write_neighbor_params "0" "$upgma" "$outgroup"
echo "# running neighbor < neighbor.params"
neighbor < neighbor.params &> /dev/null
check_output outfile
check_output outtree

# 5.1 rename outtrees and tree outfiles; remap bootstrap values to bipartitions and display tree to screen
if [ "$upgma" -gt 0 ]
then
     # https://fvue.nl/wiki/Bash:_Error_%60Unbound_variable%27_when_appending_to_empty_array
     # Set last item specifically
     # instead of appending one element, set the last item specifically, without any "unbound variable" error
     # t[${#t[*]}]=foo
     upgma_tree=
     upgma_outfile=
     
     if [ -s outtree ] 
     then
          upgma_tree="${input_phylip%.*}_${model}${gammaf}gamma_UPGMA.ph"
          mv outtree "$upgma_tree"
	  echo "# wrote tree $upgma_tree to disk" 
	  outfiles[${#outfiles[*]}]="$upgma_tree"
     fi 
     
     if [ -s outfile ]
     then
         upgma_outfile="${input_phylip%.*}_${model}${gammaf}gamma_UPGMA.outfile"
         mv outfile "$upgma_outfile"
         outfiles[${#outfiles[*]}]="$upgma_outfile"
     fi 
       
     # check that there are no negative branch lengths in the nj_tree
     #  and display with nw_display, only if no bootstrapping is requested
     if [ "$boot" -eq 0 ] && [[ $(type -P nw_display) ]] && [[ "$nw_utils_ok" -eq 1 ]]
     then 
         display_treeOK "$upgma_tree"
     elif [ "$boot" -eq 0 ] && { [[ ! $(type -P nw_display) ]] || [[ "$nw_utils_ok" -ne 1 ]]; }
     then
          echo "# extract_tree_from_outfile $upgma_outfile"
	  echo
	  extract_tree_from_outfile "$upgma_outfile"
	  echo
     fi	 
else
     nj_tree=
     nj_outfile=

     if [ -s outtree ] 
     then
	  nj_tree="${input_phylip%.*}_${model}${gammaf}gamma_NJ.ph"
          mv outtree "$nj_tree"
          echo "# wrote tree $nj_tree to disk"
          outfiles[${#outfiles[*]}]="$nj_tree"
     fi
     
     if [ -s outfile ] 
     then
          nj_outfile="${input_phylip%.*}_${model}${gammaf}gamma_NJ.outfile"
          mv outfile "$nj_outfile"
          outfiles[${#outfiles[*]}]="$nj_outfile"
     fi

     # check that there are no negative branch lengths in the nj_tree
     #  and display with nw_display, only if no bootstrapping is requested
     if [ "$boot" -eq 0 ] && [[ $(type -P nw_display) ]] && [[ "$nw_utils_ok" -eq 1 ]]
     then
         display_treeOK "$nj_tree"
     elif [ "$boot" -eq 0 ] && { [[ ! $(type -P nw_display) ]] || [[ "$nw_utils_ok" -ne 1 ]]; }
     then
         echo "# extract_tree_from_outfile $nj_outfile"
	 echo
	 extract_tree_from_outfile "$nj_outfile"
	 echo
     fi	 
fi

echo "# > finished computing distance matrix and tree for $input_phylip!"

if [ "$boot" -eq 0 ]
then
    # 6. Print final output summary message
    echo
    echo '===================== OUTPUT SUMMARY ====================='
    no_outfiles=${#outfiles[@]}
    echo "# $no_outfiles output files were generated:"
    printf "%s\n" "${outfiles[@]}"

    echo
    echo -n "# FINISHED run at: "; date 
    echo "   ==> Exiting now ..."
    echo
else
    echo "=================================================================="
    echo
fi


# ----------------------------------------------------------- #
# >>>>>>>>>>>>>>>>>>>> Bootstrap Analysis <<<<<<<<<<<<<<<<<<< #
# ----------------------------------------------------------- #

# Run seqboot if requested
# run seqboot if -b > 0

if [ "$boot" -gt 0 ] 
then
     # 1. restore the original infile for the bootstrapping and standard NJ/UPGMA analysis below
     cp "$input_phylip" infile
     check_output infile

     echo
     echo ">>> Bootstrap Analysis based on $boot pseudoreplicates for $input_phylip ..."
     echo "# running  write_seqboot_params $boot $sequential"
     write_seqboot_params "$boot" "$sequential"
     echo "# running seqboot < seqboot.params &> /dev/null"
     seqboot < seqboot.params &> /dev/null
     check_output outfile
     mv outfile infile

     echo "# > computing distance matrices on $boot bootstrapped alignments ..."
     # 2. if bootstrapping, then compute consensus tree
     # we need to run dnadist or protdist, depending on runmode
     if [ "$runmode" -eq 1 ]
     then
	 echo "# running write_dnadist_params $model $boot $TiTv $sequential $gamma $CV"
    	 write_dnadist_params "$model" "$boot" "$TiTv" "$sequential" "$gamma" "$CV"
	 echo "# running dnadist < dnadist.params &> /dev/null"
    	 dnadist < dnadist.params &> /dev/null
	 check_output outfile
	 
	 # check that matrix does not contain negative values
	 check_matrix outfile
	 
	 mv outfile infile
     elif [ "$runmode" -eq 2 ]
     then
	 echo
	 echo "# running write_protdist_params $model $boot $sequential $gamma $CV"
    	 write_protdist_params "$model" "$boot" "$sequential" "$gamma" "$CV"
	 echo "# running protdist < protdist.params &> /dev/null"
    	 protdist < protdist.params &> /dev/null
	 check_output outfile
	 
	 # check that matrix does not contain negative values
	 check_matrix outfile

	 mv outfile infile
     fi

     # 3. Now we have the distance matrices and we can proceed equaly for runmodes 1 and 2
     #    >>> Run neighbor
     echo "# > Computing distance trees from bootstrapped data ..."	  
     echo "# running write_neighbor_params $boot $upgma $outgroup"
     write_neighbor_params "$boot" "$upgma" "$outgroup"
     echo "# running neighbor < neighbor.params &> /dev/null"
     neighbor < neighbor.params &> /dev/null
     check_output outtree
     
     boot_trees=
     if [ -s outtree ]
     then
	 # this is the file holding the trees for the n-distance matrices for n-boot replicated alignments
	 boot_trees="${input_phylip%.*}_${model}${gammaf}gamma_${boot}bootRepl_trees.nwk"
	 cp outtree "$boot_trees"
	 
	 # append to array with +=, otherwise will complain as unset with set -u 
	 # https://fvue.nl/wiki/Bash:_Error_%60Unbound_variable%27_when_appending_to_empty_array
	 outfiles+=("$boot_trees")
     fi
     mv outtree intree
     rm outfile

     # 4. Compute consensus tree with consense
     echo "# > Computing MJR consensus tree from trees reconstructed from bootstrap pseudoreplicates ..."
     echo "# running write_consense_params"
     write_consense_params
     echo "# running consense < consense.params &> /dev/null"
     consense < consense.params &> /dev/null
     check_output outtree
     check_output outfile
     
     # variables holding consensus trees
     upgma_consensus_tree=
     upgma_consensus_outfile=
     nj_consensus_tree=
     nj_consensus_outfile=

     if [ "$upgma" -gt 0 ]
     then
	 # https://fvue.nl/wiki/Bash:_Error_%60Unbound_variable%27_when_appending_to_empty_array
         # Set last item specifically
	 # instead of appending one element, set the last item specifically, without any "unbound variable" error
	 # t[${#t[*]}]=foo
	 upgma_consensus_tree="${input_phylip%.*}_UPGMAconsensus_${model}${gammaf}gamma_${boot}bootRepl.ph"
	 mv outtree "$upgma_consensus_tree"
	 outfiles[${#outfiles[*]}]="$upgma_consensus_tree"
    	 
	 upgma_consensus_outfile="${input_phylip%.*}_UPGMAconsensus_${model}${gammaf}gamma_${boot}bootRepl.outfile"
	 mv outfile "$upgma_consensus_outfile"
	 outfiles[${#outfiles[*]}]="$upgma_consensus_outfile"
     else
         nj_consensus_tree="${input_phylip%.*}_NJconsensus_${model}${gammaf}gamma_${boot}bootRepl.ph"
	 mv outtree "$nj_consensus_tree"
	 outfiles[${#outfiles[*]}]="$nj_consensus_tree"
    	 
	 nj_consensus_outfile="${input_phylip%.*}_NJconsensus_${model}${gammaf}gamma_${boot}bootRepl.outfile"
	 mv outfile "$nj_consensus_outfile"
	 outfiles[${#outfiles[*]}]="$nj_consensus_outfile"
     fi 

     # 5. Rename outtrees and tree outfiles; remap bootstrap values to bipartitions and display tree on screen
     if [ "$upgma" -gt 0 ]
     then
          # if we requested bootstrapping, map bootstrap values onto UPGMA tree using
          #   nw_support upgma.ph bootRepl_tree.ph > UPGMA_with_boot_support.ph
          if [ -s "$upgma_tree" ] && [ -s "$boot_trees" ] && [[ $(type -P nw_support) ]] && [ "$nw_utils_ok" -eq 1 ]
          then
	      upgma_tree_with_boot="${input_phylip%.*}_${model}${gammaf}gamma_UPGMA_with_${boot}boot_support.ph"
	      echo "# mapping bootstrap values on UPGMA tree with nw_support ..."
	      nw_support "$upgma_tree" "$boot_trees" > "$upgma_tree_with_boot"

	      if [ -s "$upgma_tree_with_boot" ] && [[ $(type -P nw_display) ]] && [ "$nw_utils_ok" -eq 1 ]
	      then
	          outfiles[${#outfiles[*]}]="$upgma_tree_with_boot"
	          
		  # check that there are no negative branch lengths in the nj_tree 
		  #   before displaying with nw_display
		  display_treeOK "$upgma_tree_with_boot"
	      fi
	  elif [ -s "$upgma_outfile" ] && [ -s "$upgma_consensus_outfile" ] && { [[ ! $(type -P nw_support) ]] || [ "$nw_utils_ok" -ne 1 ]; }
	  then
               echo "# extract_tree_from_outfile $upgma_outfile"
	       echo
	       extract_tree_from_outfile "$upgma_outfile"
	       echo
	       echo "# extract_tree_from_outfile $upgma_consensus_outfile"
	       echo
	       extract_tree_from_outfile "$upgma_consensus_outfile"
	       echo
	  fi
     else
         # if we requested bootstrapping, map bootstrap values onto NJ tree using
         #   nw_support NJ.ph bootRepl_tree.ph > NJ_with_boot_support.ph
         if [ -s "$nj_tree" ] && [ -s "$boot_trees" ] && [[ $(type -P nw_support) ]] && [ "$nw_utils_ok" -eq 1 ]
         then
             nj_tree_with_boot="${input_phylip%.*}_${model}${gammaf}gamma_NJ_with_${boot}boot_support.ph"
	     echo "# mapping bootstrap values on NJ tree with nw_support ..."
	     nw_support "$nj_tree" "$boot_trees" > "$nj_tree_with_boot"
	 
	     check_output "$nj_tree_with_boot"
	 
	     if [ -s "$nj_tree_with_boot" ] && [[ $(type -P nw_display) ]] && [ "$nw_utils_ok" -eq 1 ]
	     then
	         outfiles+=("$nj_tree_with_boot")
	         
		 # check that there are no negative branch lengths in the nj_tree 
		 #   before displaying with nw_display
		 display_treeOK "$nj_tree_with_boot"
	     fi
    	 elif [ -s "$nj_outfile" ] && [ -s "$nj_consensus_outfile" ] && { [[ ! $(type -P nw_support) ]] || [ "$nw_utils_ok" -ne 1 ]; }
	 then

                echo "# extract_tree_from_outfile $nj_outfile"
	        echo
	        extract_tree_from_outfile "$nj_outfile"
		echo

		echo "# extract_tree_from_outfile $nj_consensus_outfile"
		echo
		extract_tree_from_outfile "$nj_consensus_outfile"
		echo
         fi
    fi
fi

# 6. Tidy up: remove the *params files and other temporary files 
# that could interfere with future runs and litter the directory

[ "$DEBUG" -eq 0 ] && cleanup_dir

# 7. Print final output summary message
echo
echo '===================== OUTPUT SUMMARY ====================='

no_outfiles=${#outfiles[@]}
echo "# $no_outfiles output files were generated:"
printf "%s\n" "${outfiles[@]}"

echo
echo -n "# FINISHED run at: "; date 
echo "   ==> Exiting now ..."
echo