Python LoadTree.bifurcating方法代码示例

# 需要导入模块: from cogent import LoadTree [as 别名]
# 或者: from cogent.LoadTree import bifurcating [as 别名]
def TreeAlign(model, seqs, tree=None, indel_rate=0.01, indel_length=0.01,
    ui = None, ests_from_pairwise=True, param_vals=None):
    """Returns a multiple alignment and tree.
    
    Uses the provided substitution model and a tree for determining the
    progressive order. If a tree is not provided a Neighbour Joining tree is
    constructed from pairwise distances estimated from pairwise aligning the
    sequences. If running in parallel, only the distance estimation is
    parallelised and only the master CPU returns the alignment and tree, other
    CPU's return None, None.
    
    Arguments:
        - model: a substitution model
        - seqs: a sequence collection
        - indel_rate, indel_length: parameters for the progressive pair-HMM
        - ests_from_pairwise: if no tree provided and True, the median value
          of the substitution model parameters are used
        - param_vals: named key, value pairs for model parameters. These
          override ests_from_pairwise.
    """
    _exclude_params = ['mprobs', 'rate', 'bin_switch']
    if param_vals:
        param_vals = dict(param_vals)
    else:
        param_vals = {}
    if isinstance(seqs, dict):
        seq_names = list(seqs.keys())
    else:
        seq_names = seqs.getSeqNames()
    
    two_seqs = len(seq_names) == 2
    
    if tree:
        tip_names = tree.getTipNames()
        tip_names.sort()
        seq_names.sort()
        assert tip_names == seq_names, \
            "names don't match between seqs and tree: tree=%s; seqs=%s" % \
            (tip_names, seq_names)
        ests_from_pairwise = False
    elif two_seqs:
        tree = LoadTree(tip_names=seqs.getSeqNames())
        ests_from_pairwise = False
    else:
        if ests_from_pairwise:
            est_params = [param for param in model.getParamList() \
                                    if param not in _exclude_params]
        else:
            est_params = None
        
        dcalc = EstimateDistances(seqs, model, do_pair_align=True,
                                    est_params=est_params)
        dcalc.run()
        dists = dcalc.getPairwiseDistances()
        tree = NJ.nj(dists)
    
    LF = model.makeLikelihoodFunction(tree.bifurcating(name_unnamed=True), aligned=False)
    if ests_from_pairwise and not param_vals:
        # we use the Median to avoid the influence of outlier pairs
        param_vals = {}
        for param in est_params:
            numbers = dcalc.getParamValues(param)
            print("Param Estimate Summary Stats: %s" % param)
            print(numbers.summarize())
            param_vals[param] = numbers.Median
    
    ui.display("Doing %s alignment" % ["progressive", "pairwise"][two_seqs])
    with LF.updatesPostponed():
        for param, val in list(param_vals.items()):
            LF.setParamRule(param, value=val, is_constant=True)
        LF.setParamRule('indel_rate', value=indel_rate, is_constant=True)
        LF.setParamRule('indel_length', value=indel_length, is_constant=True)
        LF.setSequences(seqs)
    edge = LF.getLogLikelihood().edge
    align = edge.getViterbiPath().getAlignment()
    info = Info()
    info["AlignParams"] = param_vals
    info["AlignParams"].update(dict(indel_length=indel_length, indel_rate=indel_rate))
    align.Info = info
    return align, tree

# 需要导入模块: from cogent import LoadTree [as 别名]
# 或者: from cogent.LoadTree import bifurcating [as 别名]
def main():

    # Parse input to get parameters 
    option_parser, opts, args =\
        parse_command_line_parameters(**script_info)
    
    tree_file = opts.input_tree
    trait_table_fp = opts.input_trait_table
    verbose = opts.verbose 
    
    #Handle parameters with more complex defaults
    if opts.output_table:
        output_table_fp = opts.output_table  
    else:
        output_table_fp = add_to_filename(trait_table_fp,"reformatted")
    
    if opts.output_tree:
        output_tree_fp = opts.output_tree  
    else:
        output_tree_fp = add_to_filename(trait_table_fp,"reformatted")

     

    if verbose:
        print "Running with options:"
        print "%s:%s" %("Tree file",tree_file)
        print "%s:%s" %("Trait table",trait_table_fp)
        print "%s:%s" %("Output tree",output_tree_fp)
        print "%s:%s" %("Output trait table",output_table_fp)
        print "%s:%s" %("Add branch length to root",opts.add_branch_length_to_root)
        print "%s:%s" %("Convert to NEXUS?",opts.convert_to_nexus)
         
    # Open output files
    output_trait_table_file = open(output_table_fp,"w+")
    output_tree_file  = open(output_tree_fp,"w+")
    

    # Begin reformatting
    
    # TODO: convert parameter candidates to params
    #delimiter = "\t"
    
    
    delimiter = " "
    min_branch_length = 0.0001
    root_name = "root"
    
    if format_for_bayestraits:
        convert_to_nexus = True
        convert_to_bifurcating = True
        filter_table_by_tree_tips = True
        filter_tree_by_table_entries = True
        enforce_min_branch_length = True
        convert_trait_floats_to_ints = True

    #Load inputs
    input_tree = LoadTree(tree_file)
    
    trait_table = open(trait_table_fp,"U")
    trait_table_lines = trait_table.readlines()
    
    #if opts.verbose:
        #print "Tree prior to reformatting:"
        #print input_tree.getNewick(with_distances=True)
        #print "Tree nodes prior to reformatting:"
        #print_node_summary_table(input_tree)
    
    # Remove taxa missing in tree tips
    if filter_table_by_tree_tips: 
        trait_table_lines = filter_table_by_presence_in_tree(input_tree,trait_table_lines,delimiter=delimiter)
    
    #Convert floating point values to ints
    if convert_trait_floats_to_ints:
        trait_table_lines = convert_trait_values(\
            trait_table_lines,conversion_fn = int,delimiter=delimiter)
   
    #Write out results
    output_trait_table_file.writelines(trait_table_lines)
    trait_table.close()
    output_trait_table_file.close()
   
    # Use the reformatted trait table to filter tree
    trait_table = open(output_table_fp,"U")
    trait_table_lines = trait_table.readlines()
    

    if filter_tree_by_table_entries:
        input_tree = filter_tree_tips_by_presence_in_table(input_tree,trait_table_lines,delimiter=delimiter) 
    
    # Tree reformatting
    
    
    if convert_to_bifurcating:
        input_tree = input_tree.bifurcating() # Required by most ancSR programs
    
   
    if convert_to_bifurcating:
        input_tree = ensure_root_is_bifurcating(input_tree)
        # The below nutty-looking re-filtering step is necessary
        # When ensuring the root is bifurcating, internal nodes can get moved to the tips
#.........这里部分代码省略.........