From fdc28d8fe43af63d74044167fbb91c1245620593 Mon Sep 17 00:00:00 2001 From: Suman Neupane Date: Tue, 9 May 2017 16:39:03 -0400 Subject: [PATCH] some modifications to avoid redundancy --- brnlenMCMC.py | 15 ++++++++++----- readSeq.py | 6 +++--- treeLike.py | 11 +++++++++-- treeLike_JCGAMMA.py | 15 +++++++++++++-- 4 files changed, 35 insertions(+), 12 deletions(-) diff --git a/brnlenMCMC.py b/brnlenMCMC.py index 2034c52..87c6160 100644 --- a/brnlenMCMC.py +++ b/brnlenMCMC.py @@ -14,6 +14,8 @@ from math import exp, log ########################################################################################## +tree_file_name = 'tree.tre' +sequence_file = 'example3.nex' alpha = 0.5 #gamma shape parameter for rate categories n_gen = 50000 save_every = 50 @@ -129,7 +131,7 @@ def allocatePartial(node, patterns, rates): node.partial[i*16+l*4 + 1] = pfromC_lchild*pfromC_rchild ####################################################### - +# pGA = pdiff*(node.lchild.partial[i*16+l*4 + 0]) pGC = pdiff*(node.lchild.partial[i*16+l*4 + 1]) pGG = psame*(node.lchild.partial[i*16+l*4 + 2]) @@ -166,6 +168,7 @@ def allocatePartial(node, patterns, rates): log_like = sum(like_list) return log_like + def mcmcbrn(postorder, patterns, rates): nodes = readnewick(treenewick()) @@ -195,6 +198,7 @@ def mcmcbrn(postorder, patterns, rates): print '**************************' newf.write('\n') +# print for r in range(n_gen): for i in range(len(postorder)): preedgelen = nodes[i].edgelen @@ -241,14 +245,14 @@ def mcmcbrn(postorder, patterns, rates): print print '**************************' -# newf.flush() + newf.flush() mcmc+=1 def treenewick(): script_dir = os.path.dirname(os.path.realpath(sys.argv[0])) - path = os.path.join(script_dir, 'tree.tre') + path = os.path.join(script_dir, tree_file_name) with open(path, 'r') as content: newick = content.read() return newick @@ -485,5 +489,6 @@ if __name__ == '__main__': # yuletree = '(((1:0.54019,(5:0.40299,10:0.40299):0.1372):0.72686,(6:0.10576,4:0.10576):1.16129):0.42537,(2:0.58122,(9:0.21295,(7:0.16691,(8:0.14622,3:0.14622):0.02069):0.04604):0.36827):1.1112)' rates_list = gammaRates(alpha) postorder = readnewick(treenewick()) - result = prepareTree(postorder, readSeq.patterns(), rates_list) - result2 = mcmcbrn(postorder, readSeq.patterns(), rates_list) + result = prepareTree(postorder, readSeq.patterns(sequence_file), rates_list) +# try1 = readSeq.patterns() + result2 = mcmcbrn(postorder, readSeq.patterns(sequence_file), rates_list) diff --git a/readSeq.py b/readSeq.py index e7ece02..91c58a5 100644 --- a/readSeq.py +++ b/readSeq.py @@ -1,4 +1,4 @@ -def patterns(): +def patterns(sequence_file): # import re, os, glob, itertools, fnmatch, sys, shutil from itertools import combinations @@ -10,7 +10,7 @@ def patterns(): genes = [] data = {} # print 'Reading nexus files...' - for filename in glob.glob(os.path.join(script_dir, '*.nex')): + for filename in glob.glob(os.path.join(script_dir, sequence_file)): m = re.match('(.+).nex', os.path.basename(filename)) gene_name = m.group(1) @@ -24,7 +24,7 @@ def patterns(): m = re.search('nchar\s*=\s*(\d+)', f, re.M | re.S) nchar = int(m.group(1)) - print 'nchar=', nchar +# print 'nchar=', nchar m = re.search('Matrix\s+(.+?);', f, re.M | re.S) matrix = m.group(1).strip() diff --git a/treeLike.py b/treeLike.py index c3007fa..1142e5e 100644 --- a/treeLike.py +++ b/treeLike.py @@ -9,6 +9,13 @@ import random import re, os, itertools, sys, glob from itertools import chain from math import exp, log + +########################################################################################## +tree_file_name = 'tree.tre' +sequence_file = 'example3.nex' +########################################################################################## + + class node(object): def __init__(self, ndnum): # initialization function self.rsib = None # right sibling @@ -140,7 +147,7 @@ def allocatePartial(node, patterns): def treenewick(): script_dir = os.path.dirname(os.path.realpath(sys.argv[0])) - path = os.path.join(script_dir, 'tree.tre') + path = os.path.join(script_dir, tree_file_name) with open(path, 'r') as content: newick = content.read() return newick @@ -364,4 +371,4 @@ if __name__ == '__main__': # yuletree = '(((1:0.03915,5:0.03915):0.387,(4:0.42253,2:0.42253):0.004):0.118,3:0.54433)' postorder = readnewick(treenewick()) - result = prepareTree(postorder, readSeq.patterns()) \ No newline at end of file + result = prepareTree(postorder, readSeq.patterns(sequence_file)) \ No newline at end of file diff --git a/treeLike_JCGAMMA.py b/treeLike_JCGAMMA.py index b5bc82f..6a777cf 100644 --- a/treeLike_JCGAMMA.py +++ b/treeLike_JCGAMMA.py @@ -10,6 +10,17 @@ import re, os, itertools, sys, glob from itertools import chain from scipy.stats import gamma from math import exp, log + +########################################################################################## +tree_file_name = 'tree.tre' +sequence_file = 'example3.nex' +alpha = 0.5 #gamma shape parameter for rate categories +########################################################################################## + + + + + class node(object): def __init__(self, ndnum): # initialization function self.rsib = None # right sibling @@ -160,7 +171,7 @@ def allocatePartial(node, patterns, rates): def treenewick(): script_dir = os.path.dirname(os.path.realpath(sys.argv[0])) - path = os.path.join(script_dir, 'treecopy.tre') + path = os.path.join(script_dir, tree_file_name) with open(path, 'r') as content: newick = content.read() return newick @@ -399,4 +410,4 @@ if __name__ == '__main__': # yuletree = '(((1:0.54019,(5:0.40299,10:0.40299):0.1372):0.72686,(6:0.10576,4:0.10576):1.16129):0.42537,(2:0.58122,(9:0.21295,(7:0.16691,(8:0.14622,3:0.14622):0.02069):0.04604):0.36827):1.1112)' rates_list = gammaRates(alpha) postorder = readnewick(treenewick()) - result = prepareTree(postorder, readSeq.patterns(), rates_list) + result = prepareTree(postorder, readSeq.patterns(sequence_file), rates_list)