treetime/wrappers.py
def timetree(params):
    """
    implementeing treetime tree
    """
    if params.relax is None:
        relaxed_clock_params = None
    elif params.relax==[]:
        relaxed_clock_params=True
    elif len(params.relax)==2:
        relaxed_clock_params={'slack':params.relax[0], 'coupling':params.relax[1]}


    dates = utils.parse_dates(params.dates)
    if len(dates)==0:
        print("No valid dates -- exiting.")
        return 1

    if assure_tree(params, tmp_dir='timetree_tmp'):
        print("No tree -- exiting.")
        return 1

    outdir = get_outdir(params, '_treetime')

    gtr = create_gtr(params)
    infer_gtr = params.gtr=='infer'

    ###########################################################################
    ### READ IN VCF
    ###########################################################################
    #sets ref and fixed_pi to None if not VCF
    aln, ref, fixed_pi = read_if_vcf(params)
    is_vcf = True if ref is not None else False
    branch_length_mode = params.branch_length_mode
    #variable-site-only trees can have big branch lengths, the auto setting won't work.
    if is_vcf or (params.aln and params.sequence_length):
        if branch_length_mode == 'auto':
            branch_length_mode = 'joint'



    ###########################################################################
    ### SET-UP and RUN
    ###########################################################################
    if params.aln is None and params.sequence_length is None:
        print("one of arguments '--aln' and '--sequence-length' is required.", file=sys.stderr)
        return 1
    myTree = TreeTime(dates=dates, tree=params.tree, ref=ref,
                      aln=aln, gtr=gtr, seq_len=params.sequence_length,
                      verbose=params.verbose)
    myTree.tip_slack=params.tip_slack
    if not myTree.one_mutation:
        print("TreeTime setup failed, exiting")
        return 1

    # coalescent model options
    try:
        coalescent = float(params.coalescent)
        if coalescent<10*myTree.one_mutation:
            coalescent = None
    except:
        if params.coalescent in ['opt', 'const', 'skyline']:
            coalescent = params.coalescent
        else:
            print("unknown coalescent model specification, has to be either "
                  "a float, 'opt', 'const' or 'skyline' -- exiting")
            return 1

    # determine whether confidence intervals are to be computed and how the
    # uncertainty in the rate estimate should be treated
    calc_confidence = params.confidence
    if params.clock_std_dev:
        vary_rate = params.clock_std_dev if calc_confidence else False
    elif params.confidence and params.covariation:
        vary_rate = True
    elif params.confidence:
        print("\nOutside of covariance aware mode TreeTime cannot estimate confidence intervals "
                "without specified standard deviation of the clock rate Please specify '--clock-std-dev' "
                "or rerun with '--covariance'. Will proceed without confidence estimation")
        vary_rate = False
        calc_confidence = False
    else:
        vary_rate = False

    # RUN
    root = None if params.keep_root else params.reroot
    success = myTree.run(root=root, relaxed_clock=relaxed_clock_params,
               resolve_polytomies=(not params.keep_polytomies),
               Tc=coalescent, max_iter=params.max_iter,
               fixed_clock_rate=params.clock_rate,
               n_iqd=params.clock_filter,
               time_marginal="assign" if calc_confidence else False,
               vary_rate = vary_rate,
               branch_length_mode = branch_length_mode,
               fixed_pi=fixed_pi,
               use_covariation = params.covariation)
    if success==ttconf.ERROR: # if TreeTime.run failed, exit
        print("\nTreeTime run FAILED: please check above for errors and/or rerun with --verbose 4.\n")
        return 1

    ###########################################################################
    ### OUTPUT and saving of results
    ###########################################################################
    if infer_gtr:
        print('\nInferred GTR model:')
        print(myTree.gtr)

    print(myTree.date2dist)

    basename = get_basename(params, outdir)
    if coalescent in ['skyline', 'opt', 'const']:
        print("Inferred coalescent model")
        if coalescent=='skyline':
            print_save_plot_skyline(myTree, plot=basename+'skyline.pdf', save=basename+'skyline.tsv', screen=True)
        else:
            Tc = myTree.merger_model.Tc.y[0]
            print(" --T_c: \t %1.2e \toptimized inverse merger rate in units of substitutions"%Tc)
            print(" --T_c: \t %1.2e \toptimized inverse merger rate in years"%(Tc/myTree.date2dist.clock_rate))
            print(" --N_e: \t %1.2e \tcorresponding 'effective population size' assuming 50 gen/year\n"%(Tc/myTree.date2dist.clock_rate*50))

    # plot
    import matplotlib.pyplot as plt
    from .treetime import plot_vs_years
    leaf_count = myTree.tree.count_terminals()
    label_func = lambda x: (x.name if x.is_terminal() and ((leaf_count<30
                                        and (not params.no_tip_labels))
                                      or params.tip_labels) else '')

    plot_vs_years(myTree, show_confidence=False, label_func=label_func,
                  confidence=0.9 if params.confidence else None)
    tree_fname = (outdir + params.plot_tree)
    plt.savefig(tree_fname)
    print("--- saved tree as \n\t %s\n"%tree_fname)

    plot_rtt(myTree, outdir + params.plot_rtt)
    if params.relax:
        fname = outdir+'substitution_rates.tsv'
        print("--- wrote branch specific rates to\n\t %s\n"%fname)
        with open(fname, 'w') as fh:
            fh.write("#node\tclock_length\tmutation_length\trate\tfold_change\n")
            for n in myTree.tree.find_clades(order="preorder"):
                if n==myTree.tree.root:
                    continue
                g = n.branch_length_interpolator.gamma
                fh.write("%s\t%1.3e\t%1.3e\t%1.3e\t%1.2f\n"%(n.name, n.clock_length, n.mutation_length, myTree.date2dist.clock_rate*g, g))

    export_sequences_and_tree(myTree, basename, is_vcf, params.zero_based,
                              timetree=True, confidence=calc_confidence)

    return 0
def timetree(params):
    """
    implementeing treetime tree
    """
    if params.relax is None:
        relaxed_clock_params = None
    elif params.relax==[]:
        relaxed_clock_params=True
    elif len(params.relax)==2:
        relaxed_clock_params={'slack':params.relax[0], 'coupling':params.relax[1]}


    dates = utils.parse_dates(params.dates)
    if len(dates)==0:
        print("No valid dates -- exiting.")
        return 1

    if assure_tree(params, tmp_dir='timetree_tmp'):
        print("No tree -- exiting.")
        return 1

    outdir = get_outdir(params, '_treetime')

    gtr = create_gtr(params)
    infer_gtr = params.gtr=='infer'

    ###########################################################################
    ### READ IN VCF
    ###########################################################################
    #sets ref and fixed_pi to None if not VCF
    aln, ref, fixed_pi = read_if_vcf(params)
    is_vcf = True if ref is not None else False
    branch_length_mode = params.branch_length_mode
    #variable-site-only trees can have big branch lengths, the auto setting won't work.
    if is_vcf or (params.aln and params.sequence_length):
        if branch_length_mode == 'auto':
            branch_length_mode = 'joint'



    ###########################################################################
    ### SET-UP and RUN
    ###########################################################################
    if params.aln is None and params.sequence_length is None:
        print("one of arguments '--aln' and '--sequence-length' is required.", file=sys.stderr)
        return 1
    myTree = TreeTime(dates=dates, tree=params.tree, ref=ref,
                      aln=aln, gtr=gtr, seq_len=params.sequence_length,
                      verbose=params.verbose)
    myTree.tip_slack=params.tip_slack
    if not myTree.one_mutation:
        print("TreeTime setup failed, exiting")
        return 1

    # coalescent model options
    try:
        coalescent = float(params.coalescent)
        if coalescent<10*myTree.one_mutation:
            coalescent = None
    except:
        if params.coalescent in ['opt', 'const', 'skyline']:
            coalescent = params.coalescent
        else:
            print("unknown coalescent model specification, has to be either "
                  "a float, 'opt', 'const' or 'skyline' -- exiting")
            return 1

    # determine whether confidence intervals are to be computed and how the
    # uncertainty in the rate estimate should be treated
    calc_confidence = params.confidence
    if params.clock_std_dev:
        vary_rate = params.clock_std_dev if calc_confidence else False
    elif params.confidence and params.covariation:
        vary_rate = True
    elif params.confidence:
        print("\nOutside of covariance aware mode TreeTime cannot estimate confidence intervals "
                "without specified standard deviation of the clock rate Please specify '--clock-std-dev' "
                "or rerun with '--covariance'. Will proceed without confidence estimation")
        vary_rate = False
        calc_confidence = False
    else:
        vary_rate = False

    # RUN
    root = None if params.keep_root else params.reroot
    success = myTree.run(root=root, relaxed_clock=relaxed_clock_params,
               resolve_polytomies=(not params.keep_polytomies),
               Tc=coalescent, max_iter=params.max_iter,
               fixed_clock_rate=params.clock_rate,
               n_iqd=params.clock_filter,
               time_marginal="assign" if calc_confidence else False,
               vary_rate = vary_rate,
               branch_length_mode = branch_length_mode,
               fixed_pi=fixed_pi,
               use_covariation = params.covariation)
    if success==ttconf.ERROR: # if TreeTime.run failed, exit
        print("\nTreeTime run FAILED: please check above for errors and/or rerun with --verbose 4.\n")
        return 1

    ###########################################################################
    ### OUTPUT and saving of results
    ###########################################################################
    if infer_gtr:
        print('\nInferred GTR model:')
        print(myTree.gtr)

    print(myTree.date2dist)

    basename = get_basename(params, outdir)
    if coalescent in ['skyline', 'opt', 'const']:
        print("Inferred coalescent model")
        if coalescent=='skyline':
            print_save_plot_skyline(myTree, plot=basename+'skyline.pdf', save=basename+'skyline.tsv', screen=True)
        else:
            Tc = myTree.merger_model.Tc.y[0]
            print(" --T_c: \t %1.2e \toptimized inverse merger rate in units of substitutions"%Tc)
            print(" --T_c: \t %1.2e \toptimized inverse merger rate in years"%(Tc/myTree.date2dist.clock_rate))
            print(" --N_e: \t %1.2e \tcorresponding 'effective population size' assuming 50 gen/year\n"%(Tc/myTree.date2dist.clock_rate*50))

    # plot
    import matplotlib.pyplot as plt
    from .treetime import plot_vs_years
    leaf_count = myTree.tree.count_terminals()
    label_func = lambda x: (x.name if x.is_terminal() and ((leaf_count<30
                                        and (not params.no_tip_labels))
                                      or params.tip_labels) else '')

    plot_vs_years(myTree, show_confidence=False, label_func=label_func,
                  confidence=0.9 if params.confidence else None)
    tree_fname = (outdir + params.plot_tree)
    plt.savefig(tree_fname)
    print("--- saved tree as \n\t %s\n"%tree_fname)

    plot_rtt(myTree, outdir + params.plot_rtt)
    if params.relax:
        fname = outdir+'substitution_rates.tsv'
        print("--- wrote branch specific rates to\n\t %s\n"%fname)
        with open(fname, 'w') as fh:
            fh.write("#node\tclock_length\tmutation_length\trate\tfold_change\n")
            for n in myTree.tree.find_clades(order="preorder"):
                if n==myTree.tree.root:
                    continue
                g = n.branch_length_interpolator.gamma
                fh.write("%s\t%1.3e\t%1.3e\t%1.3e\t%1.2f\n"%(n.name, n.clock_length, n.mutation_length, myTree.date2dist.clock_rate*g, g))

    export_sequences_and_tree(myTree, basename, is_vcf, params.zero_based,
                              timetree=True, confidence=calc_confidence)

    return 0
def timetree ( params ) : if params . relax is None : relaxed_clock_params = None elif params . relax == [ ] : relaxed_clock_params = True elif len ( params . relax ) == 2 : relaxed_clock_params = { 'slack' : params . relax [ 0 ] , 'coupling' : params . relax [ 1 ] } dates = utils . parse_dates ( params . dates ) if len ( dates ) == 0 : print ( "No valid dates -- exiting." ) return 1 if assure_tree ( params , tmp_dir = 'timetree_tmp' ) : print ( "No tree -- exiting." ) return 1 outdir = get_outdir ( params , '_treetime' ) gtr = create_gtr ( params ) infer_gtr = params . gtr == 'infer' ########################################################################### ### READ IN VCF ########################################################################### #sets ref and fixed_pi to None if not VCF aln , ref , fixed_pi = read_if_vcf ( params ) is_vcf = True if ref is not None else False branch_length_mode = params . branch_length_mode #variable-site-only trees can have big branch lengths, the auto setting won't work. if is_vcf or ( params . aln and params . sequence_length ) : if branch_length_mode == 'auto' : branch_length_mode = 'joint' ########################################################################### ### SET-UP and RUN ########################################################################### if params . aln is None and params . sequence_length is None : print ( "one of arguments '--aln' and '--sequence-length' is required." , file = sys . stderr ) return 1 myTree = TreeTime ( dates = dates , tree = params . tree , ref = ref , aln = aln , gtr = gtr , seq_len = params . sequence_length , verbose = params . verbose ) myTree . tip_slack = params . tip_slack if not myTree . one_mutation : print ( "TreeTime setup failed, exiting" ) return 1 # coalescent model options try : coalescent = float ( params . coalescent ) if coalescent < 10 * myTree . one_mutation : coalescent = None except : if params . coalescent in [ 'opt' , 'const' , 'skyline' ] : coalescent = params . coalescent else : print ( "unknown coalescent model specification, has to be either " "a float, 'opt', 'const' or 'skyline' -- exiting" ) return 1 # determine whether confidence intervals are to be computed and how the # uncertainty in the rate estimate should be treated calc_confidence = params . confidence if params . clock_std_dev : vary_rate = params . clock_std_dev if calc_confidence else False elif params . confidence and params . covariation : vary_rate = True elif params . confidence : print ( "\nOutside of covariance aware mode TreeTime cannot estimate confidence intervals " "without specified standard deviation of the clock rate Please specify '--clock-std-dev' " "or rerun with '--covariance'. Will proceed without confidence estimation" ) vary_rate = False calc_confidence = False else : vary_rate = False # RUN root = None if params . keep_root else params . reroot success = myTree . run ( root = root , relaxed_clock = relaxed_clock_params , resolve_polytomies = ( not params . keep_polytomies ) , Tc = coalescent , max_iter = params . max_iter , fixed_clock_rate = params . clock_rate , n_iqd = params . clock_filter , time_marginal = "assign" if calc_confidence else False , vary_rate = vary_rate , branch_length_mode = branch_length_mode , fixed_pi = fixed_pi , use_covariation = params . covariation ) if success == ttconf . ERROR : # if TreeTime.run failed, exit print ( "\nTreeTime run FAILED: please check above for errors and/or rerun with --verbose 4.\n" ) return 1 ########################################################################### ### OUTPUT and saving of results ########################################################################### if infer_gtr : print ( '\nInferred GTR model:' ) print ( myTree . gtr ) print ( myTree . date2dist ) basename = get_basename ( params , outdir ) if coalescent in [ 'skyline' , 'opt' , 'const' ] : print ( "Inferred coalescent model" ) if coalescent == 'skyline' : print_save_plot_skyline ( myTree , plot = basename + 'skyline.pdf' , save = basename + 'skyline.tsv' , screen = True ) else : Tc = myTree . merger_model . Tc . y [ 0 ] print ( " --T_c: \t %1.2e \toptimized inverse merger rate in units of substitutions" % Tc ) print ( " --T_c: \t %1.2e \toptimized inverse merger rate in years" % ( Tc / myTree . date2dist . clock_rate ) ) print ( " --N_e: \t %1.2e \tcorresponding 'effective population size' assuming 50 gen/year\n" % ( Tc / myTree . date2dist . clock_rate * 50 ) ) # plot import matplotlib . pyplot as plt from . treetime import plot_vs_years leaf_count = myTree . tree . count_terminals ( ) label_func = lambda x : ( x . name if x . is_terminal ( ) and ( ( leaf_count < 30 and ( not params . no_tip_labels ) ) or params . tip_labels ) else '' ) plot_vs_years ( myTree , show_confidence = False , label_func = label_func , confidence = 0.9 if params . confidence else None ) tree_fname = ( outdir + params . plot_tree ) plt . savefig ( tree_fname ) print ( "--- saved tree as \n\t %s\n" % tree_fname ) plot_rtt ( myTree , outdir + params . plot_rtt ) if params . relax : fname = outdir + 'substitution_rates.tsv' print ( "--- wrote branch specific rates to\n\t %s\n" % fname ) with open ( fname , 'w' ) as fh : fh . write ( "#node\tclock_length\tmutation_length\trate\tfold_change\n" ) for n in myTree . tree . find_clades ( order = "preorder" ) : if n == myTree . tree . root : continue g = n . branch_length_interpolator . gamma fh . write ( "%s\t%1.3e\t%1.3e\t%1.3e\t%1.2f\n" % ( n . name , n . clock_length , n . mutation_length , myTree . date2dist . clock_rate * g , g ) ) export_sequences_and_tree ( myTree , basename , is_vcf , params . zero_based , timetree = True , confidence = calc_confidence ) return 0
implementeing treetime tree
implementeing treetime tree