comparison phylorelatives.py @ 35:015d1fb47ec4 draft

Updated script with png plot and --major-only option

34:6935429eef05

    ape = importr('ape')
    dist = ape.dist_dna(data,model="raw",pairwise=missing_info_option)
    nj_tree = ape.nj(dist)
    return nj_tree

def ape_consensus(tree,tree_function,data,iters=1000):
    """Return majority rule consensus tree"""
    ape = importr('ape')
    tree_func = robjects.r(tree_function)
    bootstrap = ape.boot_phylo(tree, data, tree_func,
                              B=iters, quiet=1, trees=1)
    bs_trees  = bootstrap.rx('trees')[0]
    cons_tree = ape.consensus(bs_trees,p=0.5)
    return cons_tree


def dendro_relatives(tree,minor):
    """Return minor allele sequence relatives in tree"""
    ape = importr('ape')

    t = dendropy.Tree.get_from_string(newick,"newick")
    ascii_tree = t.as_ascii_plot()
    return ascii_tree


def write_nwk(tree):
    "Write proper Newick string"
    ape = importr('ape')
    nwk = list(ape.write_tree(tree))[0]
    return nwk

    parser = argparse.ArgumentParser(description='Constructs relatedness of a set of sequences based on the pairwise proportion of different sites. It reports the test sequences relatives, tree plots and trees in Newick format. One or more test sequences are accepted as long as their name includes the strict suffix "_minor" or "_test" (i.e. >seq1_minor). IMPORTANT: Sequences must have the same length!', epilog='Boris Rebolledo-Jaramillo (boris-at-bx.psu.edu)')
    parser.add_argument('-i', '--input', metavar='FASTA', action='append', type=str, help='This option can be specified multiple times. Sequences will be added to "multi-fasta.fa". (e.g. -i major1.fa -i major2.fa -i minor1.fa)')
    parser.add_argument('-b', '--bootstrap', type=int, metavar='INT',default=1000, help='Change number of replicas. 0 to deactivate. (Default: 1000)')
    parser.add_argument('-p', '--pairwise', action='store_true', help='Use pairwise deletion of gaps/missing data. (Default: Complete deletion)')
    parser.add_argument('-r', '--root', type=str, metavar='FASTA', default=False, help='Root trees using FASTA sequence as outgroup. (Default: Display unrooted trees)')
    parser.add_argument('--relatives-out', type=str, metavar='FILE', default=None, help=argparse.SUPPRESS)
    parser.add_argument('--newick-out', type=str, metavar='FILE', default=None, help=argparse.SUPPRESS)
    parser.add_argument('--trees-out', type=str, metavar='FILE', default=None, help=argparse.SUPPRESS)
    parser.add_argument('-m', '--multi-fasta', type=str, default='multi-fasta.fa', help=argparse.SUPPRESS)
    args = parser.parse_args()
    
    
    if args.input:
        for fasta in args.input:
            try:
                open(fasta)

    
    # Get sequence ids in alignment and identify test sequences
    fasta_ids = [seqname.strip()[1:] for seqname in list(open(args.multi_fasta)) if seqname.startswith('>')]
    minor_ids = [seqname for seqname in fasta_ids if seqname.endswith('minor') or seqname.endswith('test')]

    if len(minor_ids) == 0:
        sys.exit("\nERROR: No test sequences found. _minor or _test suffixes are required in the sequence name!\n")
    else:
        pass

    if args.pairwise:

    else:
        nj_tree = ape_nj(data)
        nj_func = 'function (xx) nj(dist.dna(xx, model="raw"))'
    
    if args.bootstrap == 0:
        cons_tree = nj_tree
    elif args.bootstrap !=1000:
        cons_tree = ape_consensus(nj_tree,nj_func,data,iters=args.bootstrap)
    else:
        cons_tree = ape_consensus(nj_tree,nj_func,data)
    
    # Generate report, trees, and Newick strings
    if args.relatives_out is not None:
        relatives = open(args.relatives_out,'w+')
    else:
        relatives = open(args.multi_fasta+'-relatives.tab','w+')
    if args.newick_out is not None:
        newick = open(args.newick_out,'w+')
    else:
        newick = open(args.multi_fasta+'-newick.txt','w+')
    if args.trees_out is not None:
        tree_plot_file = open(args.trees_out, 'w+')
    else:
        tree_plot_file = open(args.multi_fasta+'-tree-ascii.txt', 'w+')

    newick.write('%s\n' % (write_nwk(cons_tree)))
    relatives.write('#source\tsample\trelatives\n')
    for node in minor_ids:
        nj_relatives = [relative for relative in dendro_relatives(cons_tree,node) if relative != node]
        relatives.write( 'Neighbor_Joining_tree\t%s\t%s\n' % (node,','.join(sorted(nj_relatives))) )

    tree_plot_file.write(dendro_plot(cons_tree,root=root_id))

    newick.close()
    relatives.close()

if __name__ == "__main__": main()

35:015d1fb47ec4

    ape = importr('ape')
    dist = ape.dist_dna(data,model="raw",pairwise=missing_info_option)
    nj_tree = ape.nj(dist)
    return nj_tree

def ape_bootstrap(tree,tree_function,data,iters=1000):
    """Return bootstrap tree with bootstrap values on nodes"""
    ape = importr('ape')
    tree_func = robjects.r(tree_function)
    bootstrap = ape.boot_phylo(tree, data, tree_func,
                              B=iters, quiet=1, trees=1)
    robjects.r('''
                add_bs_value = function(tree, boot, iters) {
                tree$node.label = ((boot$BP)/iters)*100
                return(tree)
                }''')
    add_bs_value = robjects.r['add_bs_value']
    bs_tree =  add_bs_value(tree, bootstrap, iters)
    #bs_trees  = bootstrap.rx('trees')[0]
    #bs_tree = ape.consensus(bs_trees,p=0.5)
    return bs_tree


def dendro_relatives(tree,minor):
    """Return minor allele sequence relatives in tree"""
    ape = importr('ape')

    t = dendropy.Tree.get_from_string(newick,"newick")
    ascii_tree = t.as_ascii_plot()
    return ascii_tree


def ape_plot_tree(outfile, tree1, tree2=None, root=False):
    """Plot tree to png file"""
    ape = importr('ape')
    graphics  = importr('graphics')
    grdevices = importr('grDevices')
    if tree2 is None:
        grdevices.png(file=outfile, width=1024, height=768)
        if root:
            ape.plot_phylo(ape.root(tree1,root),use_edge_length=0,
                          show_node_label=1,edge_width=2, font=2,
                          cex=1,underscore=0, no_margin=1)
        else:
            ape.plot_phylo(tree1,use_edge_length=0,
                          show_node_label=1,edge_width=2, font=2,
                          cex=1,underscore=0, no_margin=1)
        #graphics.title(main='Neighbor Joining')
        grdevices.dev_off()
    elif tree2 is not None:
        grdevices.png(file=outfile, width=1024, height=768)
        graphics.par(mfcol=array.array('i',[1,2]))
        if root:
            ape.plot_phylo(ape.root(tree1,root),use_edge_length=0,
                          show_node_label=1,edge_width=2, font=2,
                          cex=1,underscore=0, no_margin=1)
        else:
            ape.plot_phylo(tree1,use_edge_length=0,
                          show_node_label=1,edge_width=2, font=2,
                          cex=1,underscore=0, no_margin=1)
        graphics.title(main='Neighbor Joining', cex=1.5, font=2)
        if root:
            ape.plot_phylo(ape.root(tree2,root),use_edge_length=0,
                          show_node_label=1,edge_width=2, font=2,
                          cex=1,underscore=0, no_margin=1)
        else:
            ape.plot_phylo(tree2,use_edge_length=0,
                          show_node_label=1,edge_width=2, font=2,
                          cex=1,underscore=0, no_margin=1)
        graphics.title(main='Maximum Parsimony',cex=1.5, font=2)
        grdevices.dev_off()
    return


def write_nwk(tree):
    "Write proper Newick string"
    ape = importr('ape')
    nwk = list(ape.write_tree(tree))[0]
    return nwk

    parser = argparse.ArgumentParser(description='Constructs relatedness of a set of sequences based on the pairwise proportion of different sites. It reports the test sequences relatives, tree plots and trees in Newick format. One or more test sequences are accepted as long as their name includes the strict suffix "_minor" or "_test" (i.e. >seq1_minor). IMPORTANT: Sequences must have the same length!', epilog='Boris Rebolledo-Jaramillo (boris-at-bx.psu.edu)')
    parser.add_argument('-i', '--input', metavar='FASTA', action='append', type=str, help='This option can be specified multiple times. Sequences will be added to "multi-fasta.fa". (e.g. -i major1.fa -i major2.fa -i minor1.fa)')
    parser.add_argument('-b', '--bootstrap', type=int, metavar='INT',default=1000, help='Change number of replicas. 0 to deactivate. (Default: 1000)')
    parser.add_argument('-p', '--pairwise', action='store_true', help='Use pairwise deletion of gaps/missing data. (Default: Complete deletion)')
    parser.add_argument('-r', '--root', type=str, metavar='FASTA', default=False, help='Root trees using FASTA sequence as outgroup. (Default: Display unrooted trees)')
    parser.add_argument('-j', '--major-only', action='store_true', help='In major-only mode no minor allele sequence is required and only a NJ major alleles tree is generated (Default: Require minor allele sequences)')
    parser.add_argument('--relatives-out', type=str, metavar='FILE', default=None, help=argparse.SUPPRESS)
    parser.add_argument('--newick-out', type=str, metavar='FILE', default=None, help=argparse.SUPPRESS)
    parser.add_argument('--trees-out', type=str, metavar='FILE', default=None, help=argparse.SUPPRESS)
    parser.add_argument('-m', '--multi-fasta', type=str, default='multi-fasta.fa', help=argparse.SUPPRESS)
    args = parser.parse_args()
    #parser.print_help()
    
    if args.input:
        for fasta in args.input:
            try:
                open(fasta)

    
    # Get sequence ids in alignment and identify test sequences
    fasta_ids = [seqname.strip()[1:] for seqname in list(open(args.multi_fasta)) if seqname.startswith('>')]
    minor_ids = [seqname for seqname in fasta_ids if seqname.endswith('minor') or seqname.endswith('test')]

    if len(minor_ids) == 0 and not args.major_only:
        sys.exit("\nERROR: No test sequences found. _minor or _test suffixes are required in the sequence name!\n")
    else:
        pass

    if args.pairwise:

    else:
        nj_tree = ape_nj(data)
        nj_func = 'function (xx) nj(dist.dna(xx, model="raw"))'
    
    if args.bootstrap == 0:
        bs_tree = nj_tree
    elif args.bootstrap !=1000:
        bs_tree = ape_bootstrap(nj_tree,nj_func,data,iters=args.bootstrap)
    else:
        bs_tree = ape_bootstrap(nj_tree,nj_func,data)
    
    # Generate report, trees, and Newick strings
    if args.relatives_out is not None:
        relatives = open(args.relatives_out,'w+')
    else:
        relatives = open(args.multi_fasta+'-relatives.tab','w+')
    if args.newick_out is not None:
        newick = open(args.newick_out,'w+')
    else:
        newick = open(args.multi_fasta+'-newick.nwk','w+')
    if args.trees_out is not None:
        tree_plot_file = open(args.trees_out, 'w+')
    else:
        #tree_plot_file = open(args.multi_fasta+'-tree-ascii.txt', 'w+')
        tree_plot_file = args.multi_fasta+'-tree.png'

    newick.write('%s\n' % (write_nwk(bs_tree)))
    #tree_plot_file.write(dendro_plot(bs_tree,root=root_id))
    ape_plot_tree(tree_plot_file,bs_tree,root=root_id)
    
    if args.major_only:
        relatives.write('Major allele only mode cannot generate a report')
    else:
        relatives.write('#source\tsample\trelatives\n')
        for node in minor_ids:
            nj_relatives = [relative for relative in dendro_relatives(bs_tree,node) if relative != node]
            relatives.write( 'Neighbor_Joining_tree\t%s\t%s\n' % (node,','.join(sorted(nj_relatives))) )

    newick.close()
    relatives.close()

if __name__ == "__main__": main()