#!/usr/bin/env python3

# Copyright (C) 2017, Weizhi Song.
# songwz03@gmail.com

# TreeSAK is free software: you can redistribute it and/or modify
# it under the terms of the GNU Affero General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.

# TreeSAK is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU Affero General Public License for more details.

# You should have received a copy of the GNU Affero General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.

import sys
import warnings
import argparse
from TreeSAK.TreeSAK_config import config_dict


def version(config_dict):
    version_file = open('%s/VERSION' % config_dict['config_file_path'])
    return version_file.readline().strip()


def print_main_help():

    help_message = ''' 
                 ...::: TreeSAK v%s :::...
                 
    Marker-related
       ExtractMarkerSeq       ->  Extract marker by blastn  
       deltall                ->  Parse stdout of deltaLL.rb
       get_arCOG_seq          ->  Retrieve arCOG sequences
       AssessMarkerPA         ->  Assess Markers by P/A among groups
       SplitScore             ->  Assess markers by split score
       AssessMarkerDeltaLL    ->  Assess Markers by DeltaLL
       OMA                    ->  Prepare input files for running OMA 
       OMA2                   ->  Filter OMA predicted OGs
       filter_rename_ar53     ->  Filter rename GTDB markers
       
    Multiple Sequence Alignment
       BMGE                   ->  Run BMGE
       pruneMSA               ->  Prune MSA with alignment_pruner.pl
       recode                 ->  Recode amino acids to Dayoff 4, Dayoff 6 or SR4 categories
       fa2phy                 ->  Convert MSA format (fasta to phylip)
       phy2fa                 ->  Convert MSA format (phylip to fasta)
       SliceMSA               ->  Slice MSA by column
       ConcateMSA             ->  Concatenate MSAs
       ConvertMSA             ->  Convert MSA format
       OneLineAln             ->  One-line fasta format alignments
       SingleLinePhy          ->  Put sequences in single line in phylip format
       CS_trim                ->  (to be added) perform chi-squared trimming to reduce compositional heterogeneity
       gap_stats              ->  The percentage of gap in each sequence of a MSA

    Tree-related
       iTOL                   ->  Prepare iTOL files
       iTOL_Tanglegram        ->  Prepare iTOL Tanglegram files
       batch_itol             ->  Batch access iTOL
       iTOL_gene_tree         ->  Genome metadata to gene metadata
       iTOL_msa_stats         ->  iTOL_msa_stats
       PB                     ->  Infer tree with PhyloBayes-MPI 
       assessPB               ->  Compare PhyloBayes chains
       supertree              ->  Infer species tree from multiple gene trees 
       MarkerSeq2Tree         ->  Marker sequence to tree
       MarkerRef2Tree         ->  Marker (reference sequence) to Tree
       GTDB_tree              ->  Get GTDB tree
       subset                 ->  Subset tree
       TaxonTree              ->  Subset GTDB tree by taxon
       compare_trees          ->  Compare trees with Mantel test
       rename_leaves          ->  Rename tree leaves
       print_leaves           ->  print out tree leaves
       FLN                    ->  Format leaf names (e.g. remove spaces in names)
       ModifyTopo             ->  Modify tree topology
       GeneRax                ->  (to be added) Run GeneRax    
       ALE                    ->  Modules for running ALE
       cogTree                ->  Infer tree for individual COG function
       koTree                 ->  Infer tree for individual KEGG function
       RootTree               ->  Root tree with outgroup leaves
       RootTreeGTDB           ->  Root tree by GTDB taxonomy
       LcaToLeaves            ->  Get two leaves that define an internal node
       replace_clade          ->  Replace tree clade
       GeneTree               ->  Infer gene tree
       guide_tree             ->  Prepare guide tree for iqtree
       get_lca_id             ->  Get LCA id
       tree_fmt               ->  Change tree format
       rm_leaf                ->  Remove leaf/leaves from tree
       
    Model-related
       PMSF                   ->  run iqtree with PMSF
       PPA                    ->  (to be added) Perform Posterior Predictive Analysis (across-site)
       
    Dating-related
       dating                 ->  Perform molecular dating
       CompareMCMC            ->  Compare MCMCTree outputs
       PlotMcmcNode           ->  Distribution of node's age estimation 
       VisHPD95               ->  HPD95 of estimated node age
       pRTC                   ->  Perform probabilistic RTC dating
       mcmcTC                 ->  Adding time constraints to mcmctree tree
       mcmc2tree              ->  Get the tree with internal node from mcmctree output
       mcmctree_out           ->  Parse MCMCTree produced .out file
       reltime                ->  Parse RelTime produced .out file
       
    Phylo-related stats
       PhyloBiAssoc            ->  A wrapper for binaryPGLMM test

    # Upgrade with: pip3 install --upgrade TreeSAK
    ''' % version(config_dict)

    print(help_message)


if __name__ == '__main__':

    ########################################################################################### initialize subparsers ############################################################################################

    # initialize the options parser
    parser = argparse.ArgumentParser()
    subparsers = parser.add_subparsers(help="--", dest='subparser_name')
    warnings.filterwarnings('ignore')

    # parse options
    if (len(sys.argv) == 1) or (sys.argv[1] in ['-h', '-H', '-help', '-Help', '--help', '--Help']):
        print_main_help()
        sys.exit(0)

    elif sys.argv[1] == 'subset':
        from TreeSAK import subset
        subset_parser = subparsers.add_parser('subset', description='Subset tree', usage=subset.subset_usage)
        subset_parser.add_argument('-i',   required=True,                           help='input tree file')
        subset_parser.add_argument('-fi',  required=False, default=1, type=int,     help='input tree format, default: 1')
        subset_parser.add_argument('-k',   required=False, default=None,            help='leaves to keep')
        subset_parser.add_argument('-r',   required=False, default=None,            help='leaves to remove')
        subset_parser.add_argument('-o',   required=True,                           help='output tree file')
        subset_parser.add_argument('-fo',  required=False, default=1, type=int,     help='output tree format, default: 1')
        args = vars(parser.parse_args())
        subset.subset(args)

    elif sys.argv[1] == 'label_tree':
        from TreeSAK import label_tree
        label_tree_parser = subparsers.add_parser('label_tree', description='Add labels to tree leaves', usage=label_tree.label_tree_usage)
        label_tree_parser.add_argument('-tree',             required=True,                                  help='tree file in newick format')
        label_tree_parser.add_argument('-label',            required=False,  default=None,                  help='label file (label,leaf)')
        label_tree_parser.add_argument('-taxon',            required=False,  default=None,                  help='taxonomic classification')
        label_tree_parser.add_argument('-rank',             required=False,  default=None,                  help='taxonomic rank to label')
        args = vars(parser.parse_args())
        label_tree.label_tree(args, config_dict)

    elif sys.argv[1] == 'OneLineAln':
        from TreeSAK import OneLineAln
        OneLineAln_parser = subparsers.add_parser('OneLineAln', description='One-line fasta format alignments', usage=OneLineAln.OneLineAln_usage)
        OneLineAln_parser.add_argument('-in',               required=True,                       help='input MSA in fasta format')
        OneLineAln_parser.add_argument('-out',              required=False, default=None,        help='output file')
        OneLineAln_parser.add_argument('-upper',            required=False, action='store_true', help='turn to uppercase')
        args = vars(parser.parse_args())
        OneLineAln.OneLineAln(args)

    elif sys.argv[1] == 'SliceMSA':
        from TreeSAK import SliceMSA
        SliceMSA_parser = subparsers.add_parser('SliceMSA', description='Slice MSA by column', usage=SliceMSA.SliceMSA_usage)
        SliceMSA_parser.add_argument('-i',                  required=True,                         help='input MSA in fasta format')
        SliceMSA_parser.add_argument('-fi',                 required=False, default='fasta',       help='format (NOT file extension) of input MSA, default: fasta')
        SliceMSA_parser.add_argument('-s',                  required=True,                         help='columns to export, e.g. 200-300, -100, 50-')
        SliceMSA_parser.add_argument('-o',                  required=True,                         help='output file or folder')
        SliceMSA_parser.add_argument('-fo',                 required=False, default='fasta',       help='format of output MSA, select from fasta and phylip-relaxed, default: fasta')
        SliceMSA_parser.add_argument('-force',              required=False, action="store_true",   help='force overwrite existing output folder')
        args = vars(parser.parse_args())
        SliceMSA.SliceMSA(args)

    elif sys.argv[1] == 'compare_trees':
        from TreeSAK import compare_trees
        compare_trees_parser = subparsers.add_parser('compare_trees', usage=compare_trees.compare_trees_usage)
        compare_trees_parser.add_argument('-o',   required=True,                       help='output directory')
        compare_trees_parser.add_argument('-t1',  required=True,                       help='tree (folder) 1')
        compare_trees_parser.add_argument('-t2',  required=True,                       help='tree (folder) 2')
        compare_trees_parser.add_argument('-tx',  required=False, default='newick',    help='extention of tree files, default: newick')
        compare_trees_parser.add_argument('-dm',  required=False, action="store_true", help='export distance-alike matrix, obtained by subtract the similarity value from 1')
        compare_trees_parser.add_argument('-t',   required=False, type=int, default=1, help='number of threads')
        compare_trees_parser.add_argument('-tmp', required=False, action="store_true", help='keep tmp files')
        compare_trees_parser.add_argument('-f',   required=False, action="store_true", help='force overwrite')
        args = vars(parser.parse_args())
        compare_trees.compare_trees(args)

    elif sys.argv[1] == 'rename_leaves':
        from TreeSAK import rename_leaves
        rename_leaves_parser = subparsers.add_parser('rename_leaves', usage=rename_leaves.rename_leaves_usage)
        rename_leaves_parser.add_argument('-i', required=True,                       help='input tree')
        rename_leaves_parser.add_argument('-r', required=True,                       help='rename file')
        rename_leaves_parser.add_argument('-f', required=False, default=1, type=int, help='tree format, default: 1')
        rename_leaves_parser.add_argument('-o', required=True,                       help='output tree')
        args = vars(parser.parse_args())
        rename_leaves.rename_leaves(args)

    elif sys.argv[1] == 'FLN':
        from TreeSAK import format_leaf_name
        format_leaf_name_parser = subparsers.add_parser('FLN', usage=format_leaf_name.format_leaf_name_usage)
        format_leaf_name_parser.add_argument('-i',                  required=True,                          help='input tree')
        format_leaf_name_parser.add_argument('-fmt',                required=False, default=1,              help='tree format, default: 1')
        format_leaf_name_parser.add_argument('-o',                  required=True,                          help='output tree')
        format_leaf_name_parser.add_argument('-s2u',                required=False, action="store_true",    help='change space in tree leaves to underscore')
        format_leaf_name_parser.add_argument('-ns',                 required=False, action="store_true",    help='remove space from leaf names')
        format_leaf_name_parser.add_argument('-nsqm',               required=False, action="store_true",    help='remove single quotation marks from leaf names')
        format_leaf_name_parser.add_argument('-ndqm',               required=False, action="store_true",    help='remove double quotation marks from leaf names')
        args = vars(parser.parse_args())
        format_leaf_name.format_leaf_name(args)

    elif sys.argv[1] == 'AssessCVG':
        from TreeSAK import AssessCVG
        AssessCVG_parser = subparsers.add_parser('AssessCVG', usage=AssessCVG.AssessCVG_usage)
        AssessCVG_parser.add_argument('-m1', required=True, help='mcmc.txt from run 1')
        AssessCVG_parser.add_argument('-m2', required=True, help='mcmc.txt from run 2')
        AssessCVG_parser.add_argument('-o',  required=True, help='output convergence plot')
        args = vars(parser.parse_args())
        AssessCVG.AssessCVG(args)

    elif sys.argv[1] == 'deltall':
        from TreeSAK import deltall
        parse_deltall_stdout_parser = subparsers.add_parser('parse_deltall_stdout', usage=deltall.deltall_usage)
        parse_deltall_stdout_parser.add_argument('-i', required=True, help='input file (e.g., nohup.out)')
        parse_deltall_stdout_parser.add_argument('-o', required=True, help='output summary')
        args = vars(parser.parse_args())
        deltall.deltall(args)

    elif sys.argv[1] == 'get_arCOG_seq':
        from TreeSAK import get_arCOG_seq
        get_arCOG_seq_parser = subparsers.add_parser('get_arCOG_seq', usage=get_arCOG_seq.get_arCOG_seq_usage)
        get_arCOG_seq_parser.add_argument('-i',      required=True,                         help='arCOD id file, one id per line')
        get_arCOG_seq_parser.add_argument('-db_dir', required=True,                         help='database folder')
        get_arCOG_seq_parser.add_argument('-o',      required=True,                         help='output folder')
        get_arCOG_seq_parser.add_argument('-f',      required=False, action="store_true",   help='force overwrite existing output folder')
        args = vars(parser.parse_args())
        get_arCOG_seq.get_arCOG_seq(args)

    elif sys.argv[1] == 'ConvertMSA':
        from TreeSAK import ConvertMSA
        ConvertMSA_parser = subparsers.add_parser('ConvertMSA', usage=ConvertMSA.ConvertMSA_usage)
        ConvertMSA_parser.add_argument('-i',       required=True,                       help='input alignment')
        ConvertMSA_parser.add_argument('-xi',      required=False, default='aln',       help='input alignment extension')
        ConvertMSA_parser.add_argument('-fi',      required=True,                       help='input alignment format, e.g., fasta, phylip')
        ConvertMSA_parser.add_argument('-o',       required=True,                       help='output alignment')
        ConvertMSA_parser.add_argument('-xo',      required=False, default='aln',       help='output alignment extension')
        ConvertMSA_parser.add_argument('-fo',      required=True,                       help='output alignment format, e.g., fasta, phylip')
        ConvertMSA_parser.add_argument('-oneline', required=False, action="store_true", help='put sequence in single line, available if -fo is fasta')
        ConvertMSA_parser.add_argument('-nogap',   required=False, action="store_true", help='remove gaps from alignment, available if -fo is fasta')
        ConvertMSA_parser.add_argument('-f',       required=False, action="store_true", help='force overwrite existing output folder')
        args = vars(parser.parse_args())
        ConvertMSA.ConvertMSA(args)

    elif sys.argv[1] == 'MarkerRef2Tree':
        from TreeSAK import MarkerRef2Tree
        MarkerRef2Tree_parser = subparsers.add_parser('MarkerRef2Tree', usage=MarkerRef2Tree.MarkerRef2Tree_usage)
        MarkerRef2Tree_parser.add_argument('-i',               required=True,                           help='faa dir')
        MarkerRef2Tree_parser.add_argument('-x',               required=False, default='faa',           help='faa file extension, default: faa')
        MarkerRef2Tree_parser.add_argument('-m',               required=True,                           help='marker seq dir, file extension need to be faa')
        MarkerRef2Tree_parser.add_argument('-mx',              required=False, default='faa',           help='marker seq file extension, default: faa')
        MarkerRef2Tree_parser.add_argument('-g',               required=False, default=None,            help='genome group')
        MarkerRef2Tree_parser.add_argument('-c',               required=False, default='85',            help='presence absence cutoffs, default: 85')
        MarkerRef2Tree_parser.add_argument('-o',               required=True,                           help='output dir')
        MarkerRef2Tree_parser.add_argument('-e',               required=False, default='1e-30',         help='e-value cutoff, default: 1e-30')
        MarkerRef2Tree_parser.add_argument('-t',               required=True,  type=int,                help='num of threads')
        MarkerRef2Tree_parser.add_argument('-mmn',             required=False, default=1, type=int,     help='minimal marker number, default: 1')
        MarkerRef2Tree_parser.add_argument('-psiblast',        required=False, action="store_true",     help='run psiblast')
        MarkerRef2Tree_parser.add_argument('-bmge',            required=False, action="store_true",     help='trim MSA with BMGE, default is trimal')
        MarkerRef2Tree_parser.add_argument('-bmge_m',          required=False, default='BLOSUM30',      help='BMGE trim model, default: BLOSUM30')
        MarkerRef2Tree_parser.add_argument('-bmge_esc',        required=False, default='0.55',          help='BMGE entropy score cutoff, default: 0.55')
        MarkerRef2Tree_parser.add_argument('-f',               required=False, action="store_true",     help='force overwrite')
        args = vars(parser.parse_args())
        MarkerRef2Tree.MarkerRef2Tree(args)

    elif sys.argv[1] == 'AssessMarkerPA':
        from TreeSAK import AssessMarkerPA
        AssessMarkerPA_parser = subparsers.add_parser('AssessMarkerPA', usage=AssessMarkerPA.AssessMarkerPA_usage)
        AssessMarkerPA_parser.add_argument('-ta',         required=True,                               help='trimmed alignments')
        AssessMarkerPA_parser.add_argument('-tax',        required=True,                               help='extension of trimmed alignments')
        AssessMarkerPA_parser.add_argument('-aa',         required=True,                               help='faa file dir')
        AssessMarkerPA_parser.add_argument('-aax',        required=True,                               help='faa file ext')
        AssessMarkerPA_parser.add_argument('-g',          required=True,                               help='genome group')
        AssessMarkerPA_parser.add_argument('-c',          required=False, default='25-50-75-85-100',   help='cutoffs, default: 25-50-75-85-100')
        AssessMarkerPA_parser.add_argument('-o',          required=True,                               help='output dir')
        AssessMarkerPA_parser.add_argument('-f',          required=False, action="store_true",         help='force overwrite existing output folder')
        args = vars(parser.parse_args())
        AssessMarkerPA.AssessMarkerPA(args)

    elif sys.argv[1] == 'AssessMarkerDeltaLL':
        from TreeSAK import AssessMarkerDeltaLL
        AssessMarkerDeltaLL_parser = subparsers.add_parser('AssessMarkerDeltaLL', usage=AssessMarkerDeltaLL.AssessMarkerDeltaLL_usage)
        AssessMarkerDeltaLL_parser.add_argument('-deltall',     required=True,                          help='DeltaLL stdout')
        AssessMarkerDeltaLL_parser.add_argument('-o',           required=True,                          help='output dir')
        AssessMarkerDeltaLL_parser.add_argument('-c',           required=False, default='25-50-75-100', help='cutoffs, default: 25-50-75-100')
        AssessMarkerDeltaLL_parser.add_argument('-mmn',         required=False, default=20, type=int,   help='minimal marker number, default: 20')
        AssessMarkerDeltaLL_parser.add_argument('-aln',         required=True,                          help='faa file dir')
        AssessMarkerDeltaLL_parser.add_argument('-jst',         required=False, default='6',            help='threads to request in job script, for running iqtree')
        AssessMarkerDeltaLL_parser.add_argument('-qsub',        required=False, action="store_true",    help='submit job scripts')
        AssessMarkerDeltaLL_parser.add_argument('-f',           required=False, action="store_true",    help='force overwrite')
        args = vars(parser.parse_args())
        AssessMarkerDeltaLL.AssessMarkerDeltaLL(args)

    elif sys.argv[1] == 'dating_ss':
        from TreeSAK import dating_ss
        Dating_parser = subparsers.add_parser('Dating_ss', usage=dating_ss.Dating_usage)
        Dating_parser.add_argument('-deltall', required=True,                          help='DeltaLL stdout')
        Dating_parser.add_argument('-aod',     required=True,                          help='AssessMarkerDeltaLL output dir')
        Dating_parser.add_argument('-og',      required=True,                          help='outgroup leaves, one leaf id per line')
        Dating_parser.add_argument('-eu',      required=True,                          help='EU tree with time constraints')
        Dating_parser.add_argument('-o',       required=True,                          help='dating wd')
        Dating_parser.add_argument('-c',       required=False, default='25-50-75-100', help='cutoffs, default: 25-50-75-100')
        Dating_parser.add_argument('-mmn',     required=False, default=20, type=int,   help='minimal marker number, default: 20')
        Dating_parser.add_argument('-ra',      required=False, default=45, type=int,   help='root age, default: 45')
        Dating_parser.add_argument('-qsub',    required=False, action="store_true",    help='submit job scripts for getting in.BV')
        Dating_parser.add_argument('-f',       required=False, action="store_true",    help='force overwrite')
        Dating_parser.add_argument('-to_test', required=True,                          help='Settings to test')
        args = vars(parser.parse_args())
        dating_ss.dating_ss(args)

    elif sys.argv[1] == 'dating':
        from TreeSAK import dating
        dating_parser = subparsers.add_parser('dating', usage=dating.dating_usage)
        dating_parser.add_argument('-i',            required=True,                          help='input tree file')
        dating_parser.add_argument('-m',            required=True,                          help='sequence alignments')
        dating_parser.add_argument('-o',            required=True,                          help='output directory')
        dating_parser.add_argument('-p',            required=True,                          help='output prefix')
        dating_parser.add_argument('-s',            required=True,                          help='settings to compare')
        dating_parser.add_argument('-BDparas',      required=True,                          help='BDparas, either M/m or C/c')
        dating_parser.add_argument('-st',           required=False, default='2',            help='sequence type, 0 for nucleotides, 1 for codons, 2 for AAs, default: 2')
        dating_parser.add_argument('-hpc4',         required=False, action="store_true",    help='submit job on hpc4')
        dating_parser.add_argument('-hpc4_wt',      required=False, default='119:59:59',    help='walltime, default: 119:59:59')
        dating_parser.add_argument('-hpc4_a',       required=False, default='marmolecol',   help='hpc4 account, select from marmolecol and spongeholobiont, default: marmolecol')
        dating_parser.add_argument('-hpc4_q',       required=False, default='amd',          help='hpc4 queue, select from amd, intel, gpu-a30, gpu-l20 and gpu-rtx5880, default: amd')
        dating_parser.add_argument('-hpc4_conda',   required=False, default='mybase2',      help='conda environment, default: mybase2')
        dating_parser.add_argument('-f',            required=False, action="store_true",    help='force overwrite')
        args = vars(parser.parse_args())
        dating.dating(args)

    elif sys.argv[1] == 'fa2phy':
        from TreeSAK import fa2phy
        fa2phy_parser = subparsers.add_parser('fa2phy', usage=fa2phy.fa2phy_usage)
        fa2phy_parser.add_argument('-i', required=True, help='input MSA in fasta format')
        fa2phy_parser.add_argument('-o', required=True, help='output MSA in phylip format')
        args = vars(parser.parse_args())
        fa2phy.fa2phy(args)

    elif sys.argv[1] == 'print_leaves':
        from TreeSAK import print_leaves
        print_leaves_parser = subparsers.add_parser('print_leaves', usage=print_leaves.print_leaves_usage)
        print_leaves_parser.add_argument('-i',   required=True,                         help='input tree file')
        print_leaves_parser.add_argument('-fmt', required=False, default=1, type=int,   help='tree format, default is 1')
        print_leaves_parser.add_argument('-o',   required=True,                         help='output text file')
        args = vars(parser.parse_args())
        print_leaves.print_leaves(args)

    elif sys.argv[1] == 'CompareMCMC':
        from TreeSAK import CompareMCMC
        CompareMCMC_parser = subparsers.add_parser('CompareMCMC', usage=CompareMCMC.CompareMCMC_usage)
        CompareMCMC_parser.add_argument('-mx',      required=True,                          help='mcmc.txt for x axis')
        CompareMCMC_parser.add_argument('-my',      required=True,                          help='mcmc.txt for y axis')
        CompareMCMC_parser.add_argument('-lx',      required=False, default=None,           help='label for x axis')
        CompareMCMC_parser.add_argument('-ly',      required=False, default=None,           help='label for y axis')
        CompareMCMC_parser.add_argument('-max',     required=False, default=None, type=int, help='maximum axis value')
        CompareMCMC_parser.add_argument('-fs',      required=False, default=16, type=int,   help='label font size, default: 16')
        CompareMCMC_parser.add_argument('-o',       required=True,                          help='output plot')
        args = vars(parser.parse_args())
        CompareMCMC.CompareMCMC(args)

    elif sys.argv[1] == 'PlotMcmcNode':
        from TreeSAK import PlotMcmcNode
        PlotMcmcNode_parser = subparsers.add_parser('PlotMcmcNode', usage=PlotMcmcNode.PlotMcmcNode_usage)
        PlotMcmcNode_parser.add_argument('-i',  required=True,                  help='input txt file')
        PlotMcmcNode_parser.add_argument('-o',  required=True,                  help='output plot')
        args = vars(parser.parse_args())
        PlotMcmcNode.PlotMcmcNode(args)

    elif sys.argv[1] == 'VisHPD95':
        from TreeSAK import VisHPD95
        VisHPD95_parser = subparsers.add_parser('VisHPD95', usage=VisHPD95.VisHPD95_usage)
        VisHPD95_parser.add_argument('-i',      required=True,                          help='input file')
        VisHPD95_parser.add_argument('-n',      required=True,                          help='nodes to plot')
        VisHPD95_parser.add_argument('-x',      required=False, default=10,type=int,    help='plot width, default is 10')
        VisHPD95_parser.add_argument('-y',      required=False, default=6,type=int,     help='plot height, default is 6')
        VisHPD95_parser.add_argument('-o',      required=True,                          help='output plot')
        args = vars(parser.parse_args())
        VisHPD95.VisHPD95(args)

    elif sys.argv[1] == 'PMSF':
        from TreeSAK import PMSF
        PMSF_parser = subparsers.add_parser('PMSF', usage=PMSF.PMSF_usage)
        PMSF_parser.add_argument('-i',      required=True,                          help='input MSA file')
        PMSF_parser.add_argument('-gm',     required=False, default='LG+F+G',       help='iqtree model for guide tree, default: LG+F+G')
        PMSF_parser.add_argument('-m',      required=False, default='LG+C60+F+G',   help='iqtree model, default: LG+C60+F+G')
        PMSF_parser.add_argument('-mdef',   required=False, default=None,           help='model definition NEXUS file')
        PMSF_parser.add_argument('-o',      required=True,                          help='output plot')
        PMSF_parser.add_argument('-p',      required=False, default='PMSF',         help='tree prefix, default: PMSF')
        PMSF_parser.add_argument('-topo',   required=False, default=None,           help='topological constraint tree, pass to -g, default is None')
        PMSF_parser.add_argument('-t',      required=False, type=int, default=1,    help='num of threads')
        PMSF_parser.add_argument('-f',      required=False, action="store_true",    help='force overwrite')
        args = vars(parser.parse_args())
        PMSF.PMSF(args)

    elif sys.argv[1] == 'SingleLinePhy':
        from TreeSAK import SingleLinePhy
        SingleLinePhy_parser = subparsers.add_parser('SingleLinePhy', usage=SingleLinePhy.SingleLinePhy_usage)
        SingleLinePhy_parser.add_argument('-i', required=True, help='input file')
        SingleLinePhy_parser.add_argument('-o', required=True, help='output file')
        args = vars(parser.parse_args())
        SingleLinePhy.SingleLinePhy(args)

    elif sys.argv[1] == 'SplitScore':
        from TreeSAK import SplitScore
        print(SplitScore.SplitScore_usage)
        exit()

    elif sys.argv[1] == 'SplitScore1OMA':
        from TreeSAK import SplitScore1OMA
        SplitScore1OMA_parser = subparsers.add_parser('SplitScore1OMA', usage=SplitScore1OMA.SplitScore1OMA_usage)
        SplitScore1OMA_parser.add_argument('-i',   required=True,                        help='OrthologousGroups.txt, produced by OMA')
        SplitScore1OMA_parser.add_argument('-s',   required=True,                        help='OrthologousGroupsFasta, produced by OMA')
        SplitScore1OMA_parser.add_argument('-u',   required=False, default= None,        help='ID of interested genomes, no file extension')
        SplitScore1OMA_parser.add_argument('-o',   required=True,                        help='output directory')
        SplitScore1OMA_parser.add_argument('-m',   required=False, default='LG+G+I',     help='iqtree_model, default: LG+G+I')
        SplitScore1OMA_parser.add_argument('-c',   required=False, type=int, default=80, help='coverage cutoff, default: 80')
        SplitScore1OMA_parser.add_argument('-f',   required=False, action="store_true",  help='force overwrite')
        SplitScore1OMA_parser.add_argument('-jst', required=False, type=int, default=3,  help='num of threads for inferring gene tree, default: 3')
        args = vars(parser.parse_args())
        SplitScore1OMA.SplitScore1OMA(args)

    elif sys.argv[1] == 'SplitScore1':
        from TreeSAK import SplitScore1
        SplitScore1_parser = subparsers.add_parser('SplitScore1', usage=SplitScore1.SplitScore1_usage)
        SplitScore1_parser.add_argument('-i',   required=True,                          help='orthologous gene sequence')
        SplitScore1_parser.add_argument('-x',   required=True,                          help='fasta file extension')
        SplitScore1_parser.add_argument('-o',   required=True,                          help='output directory')
        SplitScore1_parser.add_argument('-u',   required=False, default=None,           help='interested genomes, no file extension')
        SplitScore1_parser.add_argument('-m',   required=False, default='LG+G+I',       help='iqtree_model, default: LG+G+I')
        SplitScore1_parser.add_argument('-c',   required=False, type=int, default=75,   help='coverage cutoff, default: 75')
        SplitScore1_parser.add_argument('-f',   required=False, action="store_true",    help='force overwrite')
        SplitScore1_parser.add_argument('-seq_named_by_gnm', required=False, action="store_true", help='specify if sequence named by genome id')
        SplitScore1_parser.add_argument('-jst', required=False, type=int, default=1,    help='num of threads for iqtree2, default: 1')
        args = vars(parser.parse_args())
        SplitScore1.SplitScore1(args)

    elif sys.argv[1] == 'SplitScore2':
        from TreeSAK import SplitScore2
        SplitScore2_parser = subparsers.add_parser('SplitScore2', usage=SplitScore2.SplitScore2_usage)
        SplitScore2_parser.add_argument('-i',                   required=True,                        help='outputs of iqtree from step 1, only needs the contree and ufboot files')
        SplitScore2_parser.add_argument('-g',                   required=True,                        help='genome group')
        SplitScore2_parser.add_argument('-k',                   required=True,                        help='genome taxon, GTDB format')
        SplitScore2_parser.add_argument('-f',                   required=False, action="store_true",  help='force overwrite')
        SplitScore2_parser.add_argument('-t',                   required=False, type=int, default=1,  help='num of threads, default: 1')
        SplitScore2_parser.add_argument('-c',                   required=False, default='25,50,75',   help='marker ranking cutoffs, default: 25,50,75')
        SplitScore2_parser.add_argument('-seq_named_by_gnm',    required=False, action="store_true",  help='named_by_gnm, specify if sequence named by gnm')
        SplitScore2_parser.add_argument('-o',                   required=True,                        help='output directory')
        args = vars(parser.parse_args())
        SplitScore2.SplitScore2(args)

    elif sys.argv[1] == 'MarkerSeq2Tree':
        from TreeSAK import MarkerSeq2Tree
        MarkerSeq2Tree_parser = subparsers.add_parser('MarkerSeq2Tree', usage=MarkerSeq2Tree.MarkerSeq2Tree_usage)
        MarkerSeq2Tree_parser.add_argument('-i',        required=True,                          help='marker seq dir')
        MarkerSeq2Tree_parser.add_argument('-x',        required=True,                          help='marker seq ext')
        MarkerSeq2Tree_parser.add_argument('-o',        required=True,                          help='output dir')
        MarkerSeq2Tree_parser.add_argument('-t',        required=False, type=int, default=1,    help='num of threads')
        MarkerSeq2Tree_parser.add_argument('-bmge',     required=False, action="store_true",    help='perform BMGE trimming on concatenated MSA')
        MarkerSeq2Tree_parser.add_argument('-bmge_m',   required=False, default='BLOSUM30',     help='BMGE trim model, default: BLOSUM30')
        MarkerSeq2Tree_parser.add_argument('-bmge_esc', required=False, default='0.55',         help='BMGE entropy score cutoff, default: 0.55')
        MarkerSeq2Tree_parser.add_argument('-prune',    required=False, default=None,           help='conservation cutoffs for alignment_pruner.pl')
        MarkerSeq2Tree_parser.add_argument('-f',        required=False, action="store_true",    help='force overwrite')
        args = vars(parser.parse_args())
        MarkerSeq2Tree.MarkerSeq2Tree(args)

    elif sys.argv[1] == 'GTDB_tree':
        from TreeSAK import GTDB_tree
        GTDB_tree_parser = subparsers.add_parser('GTDB_tree', usage=GTDB_tree.GTDB_tree_usage)
        GTDB_tree_parser.add_argument('-p', required=True,                         help='output prefix')
        GTDB_tree_parser.add_argument('-i', required=True,                         help='genome folder')
        GTDB_tree_parser.add_argument('-x', required=True,                         help='genome file extension')
        GTDB_tree_parser.add_argument('-t', required=False, type=int, default=1,   help='number of threads')
        args = vars(parser.parse_args())
        GTDB_tree.GTDB_tree(args)

    elif sys.argv[1] == 'ExtractMarkerSeq':
        from TreeSAK import ExtractMarkerSeq
        ExtractMarkerSeq_parser = subparsers.add_parser('ExtractMarkerSeq', usage=ExtractMarkerSeq.ExtractMarkerSeq_usage)
        ExtractMarkerSeq_parser.add_argument('-m',               required=True,                          help='marker seq dir')
        ExtractMarkerSeq_parser.add_argument('-mx',              required=True,                          help='marker seq ext')
        ExtractMarkerSeq_parser.add_argument('-aa',              required=True,                          help='faa file dir')
        ExtractMarkerSeq_parser.add_argument('-aax',             required=True,                          help='faa file ext')
        ExtractMarkerSeq_parser.add_argument('-o',               required=True,                          help='output dir')
        ExtractMarkerSeq_parser.add_argument('-e',               required=True,  default=1e-30,          help='e-value cutoff, default: 1e-30')
        ExtractMarkerSeq_parser.add_argument('-t',               required=True,  type=int,               help='num of threads')
        ExtractMarkerSeq_parser.add_argument('-f',               required=False, action="store_true",    help='force overwrite')
        args = vars(parser.parse_args())
        ExtractMarkerSeq.ExtractMarkerSeq(args)

    elif sys.argv[1] == 'OMA':
        from TreeSAK import OMA
        OMA_parser = subparsers.add_parser('OMA', usage=OMA.OMA_usage)
        OMA_parser.add_argument('-i',   required=True,                       help='genome folder')
        OMA_parser.add_argument('-x',   required=True,                       help='genome file extension')
        OMA_parser.add_argument('-st',  required=False, default='AA',        help='sequence type, AA or DNA, default: AA')
        OMA_parser.add_argument('-og',  required=False, default=None,        help='outgroup genomes, without file extension, default is None')
        OMA_parser.add_argument('-o',   required=True,  default=None,        help='output dir, i.e., OMA working directory')
        OMA_parser.add_argument('-f',   required=False, action="store_true", help='force overwrite')
        OMA_parser.add_argument('-t',   required=False, type=int, default=6, help='number of threads for running OMA, default: 6')
        args = vars(parser.parse_args())
        OMA.OMA(args)

    elif sys.argv[1] == 'OMA2':
        from TreeSAK import OMA2
        OMA2_parser = subparsers.add_parser('OMA2', usage=OMA2.OMA2_usage)
        OMA2_parser.add_argument('-i',  required=True,                          help='OrthologousGroups.txt')
        OMA2_parser.add_argument('-s',  required=True,                          help='sequence dir, OrthologousGroupsFasta')
        OMA2_parser.add_argument('-g',  required=False, default=None,           help='interested genomes')
        OMA2_parser.add_argument('-o',  required=True,  default=None,           help='output directory')
        OMA2_parser.add_argument('-n', required=False, default=None,            help='minimal number of gene in a OG, not compatible with -c')
        OMA2_parser.add_argument('-c', required=False, default=None,            help='minimal genome coverage cutoff, not compatible with -n')
        OMA2_parser.add_argument('-f',  required=False, action="store_true",    help='force overwrite')
        args = vars(parser.parse_args())
        OMA2.OMA2(args)

    elif sys.argv[1] in ['ALE', 'ale', 'Ale']:
        from TreeSAK import ALE
        print(ALE.ALE_usage)
        exit()

    elif sys.argv[1] == 'ALE1':
        from TreeSAK import ALE1
        ALE1_parser = subparsers.add_parser('ALE1', usage=ALE1.ALE1_usage)
        ALE1_parser.add_argument('-i',          required=False, default=None,          help='orthologous groups, either from orthofinder or oma')
        ALE1_parser.add_argument('-s',          required=False, default=None,          help='sequence file, e.g., combined.faa')
        ALE1_parser.add_argument('-ms',         required=False, default=None,          help='input is a folder holds the sequence of each marker')
        ALE1_parser.add_argument('-msx',        required=False, default='fa',          help='file extension of marker sequence file, default: fa')
        ALE1_parser.add_argument('-p',          required=True,                         help='orthologous identification program, orthofinder or oma')
        ALE1_parser.add_argument('-m',          required=False, type=int, default=50,  help='min_og_genome_num, default: 50')
        ALE1_parser.add_argument('-bmge',       required=False, action="store_true",   help='trim MSA with BMGE, default no trimming')
        ALE1_parser.add_argument('-bmge_m',     required=False, default='BLOSUM30',    help='BMGE trim model, default: BLOSUM30')
        ALE1_parser.add_argument('-bmge_esc',   required=False, default='0.55',        help='BMGE entropy score cutoff, default: 0.55')
        ALE1_parser.add_argument('-o',          required=True,                         help='output dir, i.e., OMA working directory')
        ALE1_parser.add_argument('-jst',        required=False, type=int, default=3,   help='number of threads specified in job script, default: 3')
        ALE1_parser.add_argument('-f',          required=False, action="store_true",   help='force overwrite')
        args = vars(parser.parse_args())
        ALE1.ALE1(args)

    elif sys.argv[1] == 'ALE2':
        from TreeSAK import ALE2
        ALE2_parser = subparsers.add_parser('ALE2', usage=ALE2.ALE2_usage)
        ALE2_parser.add_argument('-1',      required=True,                         help='ALE1 output directory')
        ALE2_parser.add_argument('-s',      required=True,                         help='rooted species tree')
        ALE2_parser.add_argument('-o',      required=True,                         help='output dir, i.e., OMA working directory')
        ALE2_parser.add_argument('-runALE', required=False, action="store_true",   help='run ALE')
        ALE2_parser.add_argument('-docker', required=False, default=None,          help='Docker image, if ALE was installed with Docker, e.g., gregmich/alesuite_new')
        ALE2_parser.add_argument('-f',      required=False, action="store_true",   help='force overwrite')
        ALE2_parser.add_argument('-t',      required=False, type=int, default=6,   help='number of threads, default: 6')
        args = vars(parser.parse_args())
        ALE2.ALE2(args)

    elif sys.argv[1] == 'ALE3':
        from TreeSAK import ALE3
        ALE3_parser = subparsers.add_parser('ALE3', usage=ALE3.ALE3_usage)
        ALE3_parser.add_argument('-2',   required=True,                             help='Folder with uml_rec files')
        ALE3_parser.add_argument('-c',   required=False, type=float, default=75,    help='gene family presence cutoff in percentage, default: 75')
        ALE3_parser.add_argument('-a',   required=False, default=None,              help='OG functional description')
        ALE3_parser.add_argument('-o',   required=True,                             help='output dir')
        ALE3_parser.add_argument('-f',   required=False, action="store_true",       help='force overwrite')
        args = vars(parser.parse_args())
        ALE3.ALE3(args)

    elif sys.argv[1] == 'ALE4':
        from TreeSAK import ALE4
        ALE4_parser = subparsers.add_parser('ALE4', usage=ALE4.ALE4_usage)
        ALE4_parser.add_argument('-1',              required=True,                              help='ALE1 output directory')
        ALE4_parser.add_argument('-2',              required=True,                              help='ALE2 output directory')
        ALE4_parser.add_argument('-c',              required=True,                              help='genome_taxon, GTDB format')
        ALE4_parser.add_argument('-color',          required=True,                              help='phylum color code')
        ALE4_parser.add_argument('-o',              required=True,                              help='output dir, i.e., ALE4_op_dir')
        ALE4_parser.add_argument('-f',              required=False, action="store_true",        help='force overwrite')
        ALE4_parser.add_argument('-api',            required=True,                              help='iTOL API key')
        ALE4_parser.add_argument('-fc',             required=False, type=float, default=0.5,    help='hgt_freq_cutoff, default: 0.5')
        ALE4_parser.add_argument('-mld',            required=False, type=int, default=5,        help='donor_node_min_leaf_num, default: 5')
        ALE4_parser.add_argument('-mlr',            required=False, type=int, default=5,        help='recipient_node_min_leaf_num, default: 5')
        ALE4_parser.add_argument('-itol',           required=False, default='batch_access_tmp', help='iTOL project_name, default: batch_access_tmp')
        args = vars(parser.parse_args())
        ALE4.ALE4(args)

    elif sys.argv[1] == 'ALE6':
        from TreeSAK import ALE6
        ALE6_parser = subparsers.add_parser('ALE6', usage=ALE6.ALE6_usage)
        ALE6_parser.add_argument('-1',              required=True,                          help='ALE1 output directory')
        ALE6_parser.add_argument('-3',              required=True,                          help='ALE3 output directory')
        ALE6_parser.add_argument('-s',              required=True,                          help='rooted species tree')
        ALE6_parser.add_argument('-n',              required=False, default=None,           help='interested internal node(s)')
        ALE6_parser.add_argument('-cog',            required=False, default=None,           help='COG annotation results')
        ALE6_parser.add_argument('-kegg',           required=False, default=None,           help='KEGG annotation results')
        ALE6_parser.add_argument('-o',              required=True,                          help='output directory')
        ALE6_parser.add_argument('-f',              required=False, action="store_true",    help='force overwrite')
        args = vars(parser.parse_args())
        ALE6.ALE6(args)

    elif sys.argv[1] == 'ALE7':
        from TreeSAK import ALE7
        ALE7_parser = subparsers.add_parser('ALE7', usage=ALE7.ALE7_usage)
        ALE7_parser.add_argument('-6',      required=True,  help='ALE6 output directory')
        ALE7_parser.add_argument('-fun',    required=True,  help='interested functions')
        ALE7_parser.add_argument('-node',   required=True,  help='interested internal nodes')
        ALE7_parser.add_argument('-o',      required=True,  help='output directory')
        args = vars(parser.parse_args())
        ALE7.ALE7(args)

    elif sys.argv[1] == 'RootTree':
        from TreeSAK import RootTree
        RootTree_parser = subparsers.add_parser('RootTree', usage=RootTree.RootTree_usage)
        RootTree_parser.add_argument('-i',          required=True,                          help='input tree')
        RootTree_parser.add_argument('-og',         required=True,                          help='out group leaves')
        RootTree_parser.add_argument('-o',          required=True,                          help='output tree')
        RootTree_parser.add_argument('-add_root',   required=False, action='store_true',    help='add the root branch')
        RootTree_parser.add_argument('-fmt',        required=False, default=1, type=int,    help='tree format, default: 1')
        args = vars(parser.parse_args())
        RootTree.RootTree(args)

    elif sys.argv[1] == 'BMGE':
        from TreeSAK import BMGE
        BMGE_parser = subparsers.add_parser('BMGE', usage=BMGE.BMGE_usage)
        BMGE_parser.add_argument('-p',              required=True,                         help='output prefix')
        BMGE_parser.add_argument('-i',              required=True,                         help='input MSA')
        BMGE_parser.add_argument('-m',              required=False, default='BLOSUM30',    help='trim model, default: BLOSUM30')
        BMGE_parser.add_argument('-esc',            required=False, default='0.55',        help='entropy score cutoff, default: 0.55')
        args = vars(parser.parse_args())
        BMGE.BMGE(args)

    elif sys.argv[1] == 'pRTC':
        from TreeSAK import pRTC
        pRTC_parser = subparsers.add_parser('pRTC', usage=pRTC.pRTC_usage)
        pRTC_parser.add_argument('-i',              required=True,                   help='the file "out" generated by MCMCTree')
        pRTC_parser.add_argument('-m',              required=True,                   help='the file "mcmc.txt" generated by MCMCTree')
        pRTC_parser.add_argument('-r',              required=True,                   help='the folder that contains RTCs')
        pRTC_parser.add_argument('-o',              required=True,                   help='output txt file')
        pRTC_parser.add_argument('-ruby',           required=False, default='ruby',  help='path to ruby executable file, default: ruby')
        args = vars(parser.parse_args())
        pRTC.pRTC(args)

    elif sys.argv[1] == 'ConcateMSA':
        from TreeSAK import ConcateMSA
        ConcateMSA_parser = subparsers.add_parser('ConcateMSA', usage=ConcateMSA.ConcateMSA_usage)
        ConcateMSA_parser.add_argument('-i',        required=True,                          help='input MSA folder')
        ConcateMSA_parser.add_argument('-x',        required=True,                          help='input file extension')
        ConcateMSA_parser.add_argument('-p',        required=True,                          help='output prefix')
        ConcateMSA_parser.add_argument('-gene2gnm', required=False, action="store_true",    help='gene id to gnm id, split sequence id before the last _')
        args = vars(parser.parse_args())
        ConcateMSA.ConcateMSA(args)

    elif sys.argv[1] == 'SingleAleHGT':
        from TreeSAK import SingleAleHGT
        SingleAleHGT_parser = subparsers.add_parser('SingleAleHGT', usage=SingleAleHGT.SingleAleHGT_usage)
        SingleAleHGT_parser.add_argument('-faa',    required=False, default=None,               help='input aa file, e.g., OMA0001.faa')
        SingleAleHGT_parser.add_argument('-msa',    required=False, default=None,               help='input MSA file, e.g., OMA0001.aln')
        SingleAleHGT_parser.add_argument('-o',      required=True,                              help='output dir, e.g., SingleAleHGT_wd')
        SingleAleHGT_parser.add_argument('-s',      required=True,                              help='rooted species tree')
        SingleAleHGT_parser.add_argument('-c',      required=True,                              help='genome_taxon, GTDB format')
        SingleAleHGT_parser.add_argument('-color',  required=True,                              help='phylum color code')
        SingleAleHGT_parser.add_argument('-fc',     required=False, type=float, default=0.5,    help='hgt_freq_cutoff, default: 0.5')
        SingleAleHGT_parser.add_argument('-mld',    required=False, type=int, default=5,        help='donor_node_min_leaf_num, default: 5')
        SingleAleHGT_parser.add_argument('-mlr',    required=False, type=int, default=5,        help='recipient_node_min_leaf_num, default: 5')
        SingleAleHGT_parser.add_argument('-trim',   required=False, action="store_true",        help='trim MSA')
        SingleAleHGT_parser.add_argument('-docker', required=False, default=None,               help='Docker image, if ALE was installed with Docker, e.g., gregmich/alesuite_new')
        SingleAleHGT_parser.add_argument('-itol',   required=False, default='batch_access_tmp', help='iTOL project_name, default: batch_access_tmp')
        SingleAleHGT_parser.add_argument('-api',    required=True,                              help='iTOL API key')
        SingleAleHGT_parser.add_argument('-t',      required=False, type=int, default=6,        help='number of threads, default: 6')
        SingleAleHGT_parser.add_argument('-f',      required=False, action="store_true",        help='force overwrite')
        args = vars(parser.parse_args())
        SingleAleHGT.SingleAleHGT(args)

    elif sys.argv[1] == 'LcaToLeaves':
        from TreeSAK import LcaToLeaves
        LcaToLeaves_parser = subparsers.add_parser('LcaToLeaves', usage=LcaToLeaves.LcaToLeaves_usage)
        LcaToLeaves_parser.add_argument('-s',       required=True, help='the .stree file from ALE')
        LcaToLeaves_parser.add_argument('-n',       required=True, help='internal node(s)')
        args = vars(parser.parse_args())
        LcaToLeaves.LcaToLeaves(args)

    elif sys.argv[1] == 'PhyloBiAssoc':
        from TreeSAK import PhyloBiAssoc
        PhyloBiAssoc_parser = subparsers.add_parser('PhyloBiAssoc', usage=PhyloBiAssoc.PhyloBiAssoc_usage)
        PhyloBiAssoc_parser.add_argument('-i',      required=True,                       help='tree file')
        PhyloBiAssoc_parser.add_argument('-d',      required=True,                       help='data file')
        PhyloBiAssoc_parser.add_argument('-o',      required=True,                       help='output directory')
        PhyloBiAssoc_parser.add_argument('-t',      required=False, type=int, default=1, help='number of threads, default: 1')
        PhyloBiAssoc_parser.add_argument('-f',      required=False, action="store_true", help='force overwrite')
        args = vars(parser.parse_args())
        PhyloBiAssoc.PhyloBiAssoc(args)

    elif sys.argv[1] == 'replace_clade':
        from TreeSAK import replace_clade
        replace_clade_parser = subparsers.add_parser('replace_clade', usage=replace_clade.replace_clade_usage)
        replace_clade_parser.add_argument('-m',  required=True,                         help='main tree file')
        replace_clade_parser.add_argument('-mf', required=False, default=1, type=int,   help='main tree format, default: 1')
        replace_clade_parser.add_argument('-s',  required=True,                         help='subtree file')
        replace_clade_parser.add_argument('-b',  required=False, default=None,          help='replace multiple clades in batch manner')
        replace_clade_parser.add_argument('-sf', required=False, default=1, type=int,   help='subtree format, default: 1')
        replace_clade_parser.add_argument('-l',  required=True,                         help='leaves on main tree to be replaced')
        replace_clade_parser.add_argument('-o',  required=True,                         help='output tree')
        replace_clade_parser.add_argument('-of', required=False, default=9, type=int,   help='output tree format, default is 9')
        args = vars(parser.parse_args())
        replace_clade.replace_clade(args)

    elif sys.argv[1] == 'RootTreeGTDB':
        from TreeSAK import RootTreeGTDB
        RootTreeGTDB_parser = subparsers.add_parser('RootTreeGTDB', usage=RootTreeGTDB.RootTreeGTDB_usage)
        RootTreeGTDB_parser.add_argument('-tree',     required=True,                       help='input unrooted tree file or folder')
        RootTreeGTDB_parser.add_argument('-x',        required=False, default='treefile',  help='tree file extension, default is: treefile')
        RootTreeGTDB_parser.add_argument('-tax',      required=False, default='fna',       help='leaf taxon')
        RootTreeGTDB_parser.add_argument('-db',       required=True,                       help='GTDB database files')
        RootTreeGTDB_parser.add_argument('-d',        required=False, default=None,        help='domain, either ar or bac')
        RootTreeGTDB_parser.add_argument('-add_root', required=False, action='store_true', help='add the root branch')
        RootTreeGTDB_parser.add_argument('-o',        required=True,                       help='output folder')
        RootTreeGTDB_parser.add_argument('-f',        required=False, action="store_true", help='force overwrite')
        RootTreeGTDB_parser.add_argument('-r',        required=True,                       help='GTDB release, e.g., r220, r226')
        args = vars(parser.parse_args())
        RootTreeGTDB.RootTreeGTDB(args)

    elif sys.argv[1] in ['PB', 'Pb', 'pb']:
        from TreeSAK import PB
        PB_parser = subparsers.add_parser('PB', usage=PB.PB_usage)
        PB_parser.add_argument('-i',            required=True,                              help='input MSA file')
        PB_parser.add_argument('-o',            required=True,                              help='output directory')
        PB_parser.add_argument('-p',            required=True,                              help='output prefix')
        PB_parser.add_argument('-fa2plp',       required=False, action="store_true",        help='convert MSA format from fasta to phylip')
        PB_parser.add_argument('-n',            required=False, type=int, default=4,        help='number of chains, default: 4')
        PB_parser.add_argument('-t',            required=False, type=int, default=36,       help='num of cores per mpirun, default: 36')
        PB_parser.add_argument('-f',            required=False, action="store_true",        help='force overwrite')
        PB_parser.add_argument('-hpc4',         required=False, action="store_true",        help='submit to HKUST hpc4')
        PB_parser.add_argument('-hpc4_wt',      required=False, default='119:59:59',        help='hpc4 walltime, default: 119:59:59')
        PB_parser.add_argument('-hpc4_q',       required=False, default='amd',              help='hpc4 queue, default: amd')
        PB_parser.add_argument('-hpc4_a',       required=False, default='spongeholobiont',  help='hpc4 account, default: spongeholobiont')
        PB_parser.add_argument('-hpc4_conda',   required=False, default='mybase2',          help='hpc4 conda env, default: mybase2')
        args = vars(parser.parse_args())
        PB.PB(args)

    elif sys.argv[1] == 'assessPB':
        from TreeSAK import assessPB
        assessPB_parser = subparsers.add_parser('assessPB', usage=assessPB.assessPB_usage)
        assessPB_parser.add_argument('-i',      required=True,                          help='input txt file containing all the chains')
        assessPB_parser.add_argument('-o',      required=True,                          help='output directory')
        assessPB_parser.add_argument('-p',      required=False, default=None,           help='output prefix')
        assessPB_parser.add_argument('-bi',     required=False, default=1000,           help='burn-in, default: 1000')
        assessPB_parser.add_argument('-si',     required=False, default=10,             help='sample interval, default: 10')
        assessPB_parser.add_argument('-f',      required=False, action="store_true",    help='force overwrite')
        args = vars(parser.parse_args())
        assessPB.assessPB(args)

    elif sys.argv[1] == 'gap_stats':
        from TreeSAK import gap_stats
        gap_stats_parser = subparsers.add_parser('gap_stats', usage=gap_stats.gap_stats_usage)
        gap_stats_parser.add_argument('-i',  required=True,                             help='MSA in fasta format')
        gap_stats_parser.add_argument('-c',  required=False, default=100, type=float,   help='maximum gap allowed in MSAs, default is 100, i.e., no filtering')
        gap_stats_parser.add_argument('-o',  required=True,                             help='output stats file')
        gap_stats_parser.add_argument('-o2', required=False, default=None,              help='filtered MSA')
        gap_stats_parser.add_argument('-o3', required=False, default=None,              help='iTOL file')
        args = vars(parser.parse_args())
        gap_stats.gap_stats(args)

    elif sys.argv[1] == 'recode':
        from TreeSAK import recode
        recode_parser = subparsers.add_parser('recode', usage=recode.recode_usage)
        recode_parser.add_argument('-i',          required=True, help='input file')
        recode_parser.add_argument('-m',          required=True, help='recoding scheme, choose from d4, d6 or s4')
        recode_parser.add_argument('-o',          required=True, help='output file')
        args = vars(parser.parse_args())
        recode.recode(args)

    elif sys.argv[1] == 'pruneMSA':
        from TreeSAK import pruneMSA
        pruneMSA_parser = subparsers.add_parser('pruneMSA', usage=pruneMSA.pruneMSA_usage)
        pruneMSA_parser.add_argument('-i',          required=True, help='input MSA file')
        pruneMSA_parser.add_argument('-c',          required=True, help='conservation cutoffs, comma separated')
        args = vars(parser.parse_args())
        pruneMSA.pruneMSA(args)

    elif sys.argv[1] == 'iTOL':
        from TreeSAK import iTOL
        iTOL_parser = subparsers.add_parser('iTOL', description='Decorate tree with iTOL', usage=iTOL.iTOL_usage)
        iTOL_parser.add_argument('-Labels',             required=False, action='store_true',    help='Labels')
        iTOL_parser.add_argument('-ColoredLabel',       required=False, action='store_true',    help='ColoredLabel')
        iTOL_parser.add_argument('-MultiStyleLabel',    required=False, action='store_true',    help='MultiStyleLabel')
        iTOL_parser.add_argument('-ColorStrip',         required=False, action='store_true',    help='ColorStrip')
        iTOL_parser.add_argument('-ColorRange',         required=False, action='store_true',    help='ColorRange')
        iTOL_parser.add_argument('-ColorClade',         required=False, action='store_true',    help='ColorClade')
        iTOL_parser.add_argument('-ColorLabel',         required=False, action='store_true',    help='ColorLabel')
        iTOL_parser.add_argument('-SimpleBar',          required=False, action='store_true',    help='SimpleBar')
        iTOL_parser.add_argument('-Heatmap',            required=False, action='store_true',    help='Heatmap')
        iTOL_parser.add_argument('-ExternalShape',      required=False, action='store_true',    help='ExternalShape')
        iTOL_parser.add_argument('-Binary',             required=False, action='store_true',    help='Binary')
        iTOL_parser.add_argument('-BinaryID',           required=False, action='store_true',    help='Binary specified IDs as 1')
        iTOL_parser.add_argument('-BinaryShape',        required=False, default='2',            help='Binary Shape, choose from 1(rectangle), 2(circle), 3(star), 4, 5 and 6, default is 2')
        iTOL_parser.add_argument('-BinaryColor',        required=False, default='red',          help='Binary Color, default is red')
        iTOL_parser.add_argument('-Connection',         required=False, action='store_true',    help='Connection')
        iTOL_parser.add_argument('-PieChart',           required=False, action='store_true',    help='PieChart')
        iTOL_parser.add_argument('-Collapse',           required=False, action='store_true',    help='Collapse')
        iTOL_parser.add_argument('-TimeScale',          required=False, action='store_true',    help='TimeScale')
        iTOL_parser.add_argument('-id',                 required=False, default=None,           help='File contains leaf id')
        iTOL_parser.add_argument('-ll',                 required=False, default=None,           help='Leaf Label')
        iTOL_parser.add_argument('-gc',                 required=False, default=None,           help='Specify Group/column Color (optional)')
        iTOL_parser.add_argument('-cc',                 required=False, default=None,           help='Specify Column Color (for ExternalShape format) (optional)')
        iTOL_parser.add_argument('-lg',                 required=False, default=None,           help='Leaf to Group')
        iTOL_parser.add_argument('-lv',                 required=False, default=None,           help='Leaf to Value')
        iTOL_parser.add_argument('-lm',                 required=False, default=None,           help='Leaf to data Matrix')
        iTOL_parser.add_argument('-dr',                 required=False, default=None,           help='Donor to Recipient')
        iTOL_parser.add_argument('-scale',              required=False, default=None,           help='Scale Values, in format 0-3-6-9')
        iTOL_parser.add_argument('-lt',                 required=False, default=None,           help='Legend Title')
        iTOL_parser.add_argument('-hide_legend',        required=False, action='store_true',    help='hide legend for ColorStrip')
        iTOL_parser.add_argument('-show_label',         required=False, action='store_true',    help='show label on ColorStrip')
        iTOL_parser.add_argument('-strip_width',        required=False, default=100,            help='width of ColorStrip, default is 100')
        iTOL_parser.add_argument('-o',                  required=False, default=None,           help='Output filename')
        iTOL_parser.add_argument('-ct',                 required=False, default=None,           help='TimeScale, chart.ttl file')
        iTOL_parser.add_argument('-rg',                 required=False, default=None,           help='TimeScale, time range, e.g., 3.5-0')
        iTOL_parser.add_argument('-u',                  required=False, default='Ga',           help='TimeScale, time unit, choose from Ga and Ma, default is Ga')
        iTOL_parser.add_argument('-rk',                 required=False, default='Era,Eon',      help='TimeScale, interested rank, default is Era,Eon')
        args = vars(parser.parse_args())
        iTOL.iTOL(args)

    elif sys.argv[1] == 'supertree':
        from TreeSAK import supertree
        supertree_parser = subparsers.add_parser('supertree', usage=supertree.supertree_usage)
        supertree_parser.add_argument('-i',         required=True,                          help='orthologous gene sequence')
        supertree_parser.add_argument('-x',         required=True,                          help='faa file extension')
        supertree_parser.add_argument('-o',         required=True,                          help='output directory')
        supertree_parser.add_argument('-bmge',      required=False, action="store_true",    help='trim with BMGE, default is trimal')
        supertree_parser.add_argument('-bmge_m',    required=False, default='BLOSUM30',     help='trim model, default: BLOSUM30')
        supertree_parser.add_argument('-bmge_esc',  required=False, default='0.55',         help='entropy score cutoff, default: 0.55')
        supertree_parser.add_argument('-iqtree_m',  required=False, default='LG+G+I',       help='iqtree_model, default: LG+G+I')
        supertree_parser.add_argument('-pb',        required=False, action="store_true",    help='infer tree with PhyloBayes-MPI, default is iqtree')
        supertree_parser.add_argument('-f',         required=False, action="store_true",    help='force overwrite')
        supertree_parser.add_argument('-t',         required=False, type=int, default=1,    help='num of threads, default: 1')
        args = vars(parser.parse_args())
        supertree.supertree(args)

    elif sys.argv[1] == 'TaxonTree':
        from TreeSAK import TaxonTree
        TaxonTree_parser = subparsers.add_parser('TaxonTree', usage=TaxonTree.TaxonTree_usage)
        TaxonTree_parser.add_argument('-i',     required=True,  help='input tree file')
        TaxonTree_parser.add_argument('-tax',   required=True,  help='interested taxon')
        TaxonTree_parser.add_argument('-o',     required=True,  help='output tree file')
        args = vars(parser.parse_args())
        TaxonTree.TaxonTree(args)

    elif sys.argv[1] == 'mcmcTC':
        from TreeSAK import mcmcTC
        mcmcTC_parser = subparsers.add_parser('mcmcTC', usage=mcmcTC.mcmcTC_usage)
        mcmcTC_parser.add_argument('-i',    required=True, help='input tree')
        mcmcTC_parser.add_argument('-o',    required=True, help='output tree')
        mcmcTC_parser.add_argument('-tc',   required=True, help='time constraint file')
        args = vars(parser.parse_args())
        mcmcTC.mcmcTC(args)

    elif sys.argv[1] == 'mcmc2tree':
        from TreeSAK import mcmc2tree
        mcmc2tree_parser = subparsers.add_parser('mcmc2tree', usage=mcmc2tree.mcmc2tree_usage)
        mcmc2tree_parser.add_argument('-i', required=True, help='the .out file from mcmctree')
        args = vars(parser.parse_args())
        mcmc2tree.mcmc2tree(args)

    elif sys.argv[1] == 'cogTree':
        from TreeSAK import cogTree
        cogTree_parser = subparsers.add_parser('cogTree', usage=cogTree.cogTree_usage)
        cogTree_parser.add_argument('-i',         required=True,                          help='orthologous gene sequence')
        cogTree_parser.add_argument('-fun',       required=True,                          help='interested functions')
        cogTree_parser.add_argument('-cog',       required=False, default=None,           help='COG annotation results')
        cogTree_parser.add_argument('-o',         required=True,                          help='output directory')
        cogTree_parser.add_argument('-bmge',      required=False, action="store_true",    help='trim with BMGE, default is trimal')
        cogTree_parser.add_argument('-bmge_m',    required=False, default='BLOSUM30',     help='trim model, default: BLOSUM30')
        cogTree_parser.add_argument('-bmge_esc',  required=False, default='0.55',         help='entropy score cutoff, default: 0.55')
        cogTree_parser.add_argument('-iqtree_m',  required=False, default='LG+G+I',       help='iqtree_model, default: LG+G+I')
        cogTree_parser.add_argument('-f',         required=False, action="store_true",    help='force overwrite')
        cogTree_parser.add_argument('-t',         required=False, type=int, default=1,    help='num of threads, default: 1')
        args = vars(parser.parse_args())
        cogTree.cogTree(args)

    elif sys.argv[1] == 'koTree':
        from TreeSAK import koTree
        koTree_parser = subparsers.add_parser('koTree', usage=koTree.koTree_usage)
        koTree_parser.add_argument('-i',         required=True,                          help='orthologous gene sequence')
        koTree_parser.add_argument('-fun',       required=True,                          help='interested functions')
        koTree_parser.add_argument('-kegg',      required=False, default=None,           help='KEGG annotation results')
        koTree_parser.add_argument('-o',         required=True,                          help='output directory')
        koTree_parser.add_argument('-bmge',      required=False, action="store_true",    help='trim with BMGE, default is trimal')
        koTree_parser.add_argument('-bmge_m',    required=False, default='BLOSUM30',     help='trim model, default: BLOSUM30')
        koTree_parser.add_argument('-bmge_esc',  required=False, default='0.55',         help='entropy score cutoff, default: 0.55')
        koTree_parser.add_argument('-iqtree_m',  required=False, default='LG+G+I',       help='iqtree_model, default: LG+G+I')
        koTree_parser.add_argument('-f',         required=False, action="store_true",    help='force overwrite')
        koTree_parser.add_argument('-t',         required=False, type=int, default=1,    help='num of threads, default: 1')
        args = vars(parser.parse_args())
        koTree.koTree(args)

    elif sys.argv[1] == 'iTOL_gene_tree':
        from TreeSAK import iTOL_gene_tree
        iTOL_gene_tree_parser = subparsers.add_parser('iTOL_gene_tree', usage=iTOL_gene_tree.iTOL_gene_tree_usage)
        iTOL_gene_tree_parser.add_argument('-i',    required=True,                          help='input metadata')
        iTOL_gene_tree_parser.add_argument('-tree', required=False, default=None,           help='gene id, in tree file')
        iTOL_gene_tree_parser.add_argument('-txt',  required=False, default=None,           help='gene id, in txt file')
        iTOL_gene_tree_parser.add_argument('-o',    required=True,                          help='output metadata')
        iTOL_gene_tree_parser.add_argument('-na',   required=False, action='store_true',    help='include leaves with na values')
        args = vars(parser.parse_args())
        iTOL_gene_tree.iTOL_gene_tree(args)

    elif sys.argv[1] == 'GeneTree':
        from TreeSAK import GeneTree
        GeneTree_parser = subparsers.add_parser('GeneTree', usage=GeneTree.GeneTree_usage)
        GeneTree_parser.add_argument('-i',          required=False, default=None,           help='sequence file')
        GeneTree_parser.add_argument('-o',          required=True,                          help='output dir')
        GeneTree_parser.add_argument('-t',          required=False, type=int, default=1,    help='number of threads, default is 1')
        GeneTree_parser.add_argument('-trimal',     required=False, action="store_true",    help='trim with trimal, default is BMGE')
        GeneTree_parser.add_argument('-f',          required=False, action="store_true",    help='force overwrite')
        GeneTree_parser.add_argument('-max_gap',    required=False, default='40',           help='maximum percentage of gap, default is 40')
        args = vars(parser.parse_args())
        GeneTree.GeneTree(args)

    elif sys.argv[1] == 'reltime':
        from TreeSAK import reltime
        reltime_parser = subparsers.add_parser('reltime', usage=reltime.reltime_usage)
        reltime_parser.add_argument('-i',       required=True,                          help='reltime output file')
        reltime_parser.add_argument('-n',       required=False, default=None,           help='interested node txt')
        reltime_parser.add_argument('-scale',   required=False, default=1, type=int,    help='scale the value, default is 1 (no scaling)')
        reltime_parser.add_argument('-o',       required=True,                          help='output txt file')
        args = vars(parser.parse_args())
        reltime.reltime(args)

    elif sys.argv[1] == 'mcmctree_vs_reltime':
        from TreeSAK import mcmctree_vs_reltime
        mcmctree_vs_reltime_parser = subparsers.add_parser('mcmctree_vs_reltime', usage=mcmctree_vs_reltime.mcmctree_vs_reltime_usage)
        mcmctree_vs_reltime_parser.add_argument('-m', required=True, help='.out file from MCMCTree')
        mcmctree_vs_reltime_parser.add_argument('-r', required=True, help='output from elTime')
        mcmctree_vs_reltime_parser.add_argument('-n', required=True, help='interested nodes txt file')
        mcmctree_vs_reltime_parser.add_argument('-o', required=True, help='output pdf')
        args = vars(parser.parse_args())
        mcmctree_vs_reltime.mcmctree_vs_reltime(args)

    elif sys.argv[1] == 'iTOL_msa_stats':
        from TreeSAK import iTOL_msa_stats
        iTOL_msa_stats_parser = subparsers.add_parser('iTOL_msa_stats', usage=iTOL_msa_stats.iTOL_msa_stats_usage)
        iTOL_msa_stats_parser.add_argument('-i', required=True, help='MSA file')
        args = vars(parser.parse_args())
        iTOL_msa_stats.iTOL_msa_stats(args)

    elif sys.argv[1] == 'filter_rename_ar53':
        from TreeSAK import filter_rename_ar53
        filter_rename_ar53_parser = subparsers.add_parser('filter_rename_ar53', usage=filter_rename_ar53.filter_rename_ar53_usage)
        filter_rename_ar53_parser.add_argument('-i', required=True,                       help='sequence folder')
        filter_rename_ar53_parser.add_argument('-x', required=True,                       help='file extension')
        filter_rename_ar53_parser.add_argument('-g', required=False, default=None,        help='interested genome, no ext, one id per line')
        filter_rename_ar53_parser.add_argument('-m', required=False, default=None,        help='interested marker, no ext, one id per line')
        filter_rename_ar53_parser.add_argument('-o', required=True,                       help='output folder')
        filter_rename_ar53_parser.add_argument('-f', required=False, action="store_true", help='force overwrite')
        args = vars(parser.parse_args())
        filter_rename_ar53.filter_rename_ar53(args)

    elif sys.argv[1] == 'batch_itol':
        from TreeSAK import batch_itol
        batch_itol_parser = subparsers.add_parser('batch_itol', usage=batch_itol.batch_itol_usage)
        batch_itol_parser.add_argument('-i',        required=True,                              help='input tree file or folder')
        batch_itol_parser.add_argument('-x',        required=False, default=None,               help='file extension')
        batch_itol_parser.add_argument('-o',        required=True,                              help='output file or folder')
        batch_itol_parser.add_argument('-a',        required=False, default=None,               help='a txt file contain absolute to all annotation files')
        batch_itol_parser.add_argument('-para',     required=False, default=None,               help='parameter file')
        batch_itol_parser.add_argument('-api',      required=True,                              help='iTOL API key')
        batch_itol_parser.add_argument('-ip',       required=False, default='batch_access_tmp', help='iTOL project name, default: batch_access_tmp')
        batch_itol_parser.add_argument('-f',        required=False, action="store_true",        help='force overwrite')
        args = vars(parser.parse_args())
        batch_itol.batch_itol(args)

    elif sys.argv[1] == 'guide_tree':
        from TreeSAK import guide_tree
        guide_tree_parser = subparsers.add_parser('guide_tree', usage=guide_tree.guide_tree_usage)
        guide_tree_parser.add_argument('-i',    required=True,                  help='input tree')
        guide_tree_parser.add_argument('-g',    required=True,                  help='gene/genome group/taxonomy')
        guide_tree_parser.add_argument('-id',   required=False, default=None,   help='interested genes/genomes')
        guide_tree_parser.add_argument('-o',    required=True,                  help='output tree')
        args = vars(parser.parse_args())
        guide_tree.guide_tree(args)

    elif sys.argv[1] == 'get_lca_id':
        from TreeSAK import get_lca_id
        get_lca_id_parser = subparsers.add_parser('get_lca_id', usage=get_lca_id.get_lca_id_usage)
        get_lca_id_parser.add_argument('-i', required=True, help='input tree file')
        get_lca_id_parser.add_argument('-l', required=True, help='leaves file')
        args = vars(parser.parse_args())
        get_lca_id.get_lca_id(args)

    elif sys.argv[1] == 'mcmctree_out':
        from TreeSAK import mcmctree_out
        mcmctree_out_parser = subparsers.add_parser('mcmctree_out', usage=mcmctree_out.mcmctree_out_usage)
        mcmctree_out_parser.add_argument('-i', required=True, help='MCMCTree out file/dir')
        mcmctree_out_parser.add_argument('-n', required=True, help='interested nodes')
        mcmctree_out_parser.add_argument('-o', required=True, help='output table')
        args = vars(parser.parse_args())
        mcmctree_out.mcmctree_out(args)

    elif sys.argv[1] == 'tree_fmt':
        from TreeSAK import tree_fmt
        tree_fmt_parser = subparsers.add_parser('tree_fmt', usage=tree_fmt.tree_fmt_usage)
        tree_fmt_parser.add_argument('-i',   required=True,                          help='input tree')
        tree_fmt_parser.add_argument('-fi',  required=False, default=1, type=int,    help='input tree format, default: 1')
        tree_fmt_parser.add_argument('-o',   required=True,                          help='output tree')
        tree_fmt_parser.add_argument('-fo',  required=False, default=1, type=int,    help='output tree format, default: 1')
        args = vars(parser.parse_args())
        tree_fmt.tree_fmt(args)

    elif sys.argv[1] == 'rm_leaf':
        from TreeSAK import rm_leaf
        rm_leaf_parser = subparsers.add_parser('rm_leaf', usage=rm_leaf.rm_leaf_usage)
        rm_leaf_parser.add_argument('-i',   required=True,                        help='input tree file')
        rm_leaf_parser.add_argument('-if',  required=False, default=0, type=int,  help='input tree format, default is 0')
        rm_leaf_parser.add_argument('-l',   required=True,                        help='leaf/leaves to remove')
        rm_leaf_parser.add_argument('-o',   required=True,                        help='output tree file')
        rm_leaf_parser.add_argument('-of',  required=False, default=0, type=int,  help='output tree format, default is 0')
        args = vars(parser.parse_args())
        rm_leaf.rm_leaf(args)

    else:
        print('Unrecognized module: %s, program exited!' % sys.argv[1])
        exit()


upload_to_pypi_cmd = '''

cd /Users/songweizhi/PycharmProjects/TreeSAK
rm -r build dist TreeSAK.egg-info
python3 setup.py sdist bdist_wheel
twine upload dist/*

__token__

pip3 install --upgrade TreeSAK

'''
