From 2bd003a8434d859a6101301e1fb3cd1916892250 Mon Sep 17 00:00:00 2001 From: Suman Neupane Date: Sat, 6 May 2017 23:19:19 -0400 Subject: [PATCH] cleared redudant parts --- brnlenMCMC.py | 22 ++++++++++++++-------- readSeq.py | 2 +- simSequences.py | 6 +++--- tree.tre | 2 +- treeLike_JCGAMMA.py | 4 +++- 5 files changed, 22 insertions(+), 14 deletions(-) diff --git a/brnlenMCMC.py b/brnlenMCMC.py index 26a5d96..e439ca3 100644 --- a/brnlenMCMC.py +++ b/brnlenMCMC.py @@ -169,7 +169,7 @@ def replaedglen(): def mcmcbrn(postorder, patterns, rates): nodes = readnewick(treenewick()) - run = 2 + run = 1 mcmc = 0 output = os.path.join('brnlenMCMC_results.txt') newf = open(output, 'w') @@ -189,19 +189,21 @@ def mcmcbrn(postorder, patterns, rates): for r in range(run): b = 0.1 - for i,nd in enumerate(postorder): +# for i,nd in enumerate(postorder): + for i in range(len(postorder)): + print 'mcmc gen=', mcmc+1, print - # print i, nd.number # print 'current', print 'current_node.edgelen=', nodes[i].number, nodes[i].edgelen # print preedgelen = nodes[i].edgelen # print 'nodes[i]=', nodes[i].edgelen - # prelike = prepareTree(postorder, patterns, rates) currentlike = prepareTree(nodes, patterns, rates) +# currentlike = 0.0 currentprior = (-nodes[i].edgelen/b)-(log(b)) +# newf.write('%s\t'%(currentprior)) current_ln_posterior = currentlike + currentprior @@ -211,7 +213,10 @@ def mcmcbrn(postorder, patterns, rates): print 'proposed_node.edgelen=',nodes[i].number, nodes[i].edgelen # print 'new_edgelen=', nodes[i].edgelen proposedprior = (-nodes[i].edgelen/b)-(log(b)) +# newf.write('%s\t'%(proposedprior)) + proposedtlike = prepareTree(nodes, patterns, rates) +# proposedtlike = 0.0 proposed_ln_posterior = proposedtlike + proposedprior print print 'current prior, currentlike=', currentprior, currentlike @@ -248,7 +253,7 @@ def mcmcbrn(postorder, patterns, rates): def treenewick(): script_dir = os.path.dirname(os.path.realpath(sys.argv[0])) - path = os.path.join(script_dir, 'tree2.tre') + path = os.path.join(script_dir, 'tree.tre') with open(path, 'r') as content: newick = content.read() return newick @@ -462,9 +467,9 @@ def calcExpectedHeight(num_species, mu_over_s): if __name__ == '__main__': - random_seed = 146549 # 7632557, 12345 + random_seed = 346549 # 7632557, 12345 number_of_species = 5 - mutation_speciation_rate_ratio = 0.4 # 0.689655172 # yields tree height 1 for 6 species + mutation_speciation_rate_ratio = 0.689655172 # 0.689655172 # yields tree height 1 for 6 species random.seed(random_seed) species_tree_root = yuleTree(number_of_species, mutation_speciation_rate_ratio) # print '#########' @@ -479,11 +484,12 @@ def calcExpectedHeight(num_species, mu_over_s): # print ' expected height =',expected_height actual_height = calcActualHeight(species_tree_root) # print ' actual height =',actual_height -# print ' newick: ',newick[0] + print ' newick: ',newick[0] alpha = 0.5 ### gamma shape parameter rate categories # 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) +# try1 = readSeq.patterns() result2 = mcmcbrn(postorder, readSeq.patterns(), rates_list) diff --git a/readSeq.py b/readSeq.py index 6453201..e7ece02 100644 --- a/readSeq.py +++ b/readSeq.py @@ -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/simSequences.py b/simSequences.py index bf3cca1..4181da1 100644 --- a/simSequences.py +++ b/simSequences.py @@ -165,10 +165,10 @@ def readnewick(tree): output_filename = os.path.join('simulated_output.nexus') - yuletree = '(((1:0.54019,(5:0.40299,10:0.40299):0.13720):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.11120)' + yuletree = '(5:1.86010,((3:0.47109,2:0.47109):0.492,(4:0.05805,1:0.05805):0.906):0.896)' preorder = readnewick(yuletree) - ntax = 10 - num_sites = 10000 + ntax = 5 + num_sites = 100 result = simulate(preorder, ntax, num_sites, output_filename) diff --git a/tree.tre b/tree.tre index 9e8d81b..2d84738 100644 --- a/tree.tre +++ b/tree.tre @@ -1 +1 @@ -(((1:0.00039,5:0.03915):0.987,(4:0.42253,2:0.000042):0.84):0.118,3:0.1433) \ No newline at end of file +(5:1.8601,((3:0.47109,2:0.47109):0.492,(4:0.05805,1:0.05805):0.906):0.896) \ No newline at end of file diff --git a/treeLike_JCGAMMA.py b/treeLike_JCGAMMA.py index 4011773..b5bc82f 100644 --- a/treeLike_JCGAMMA.py +++ b/treeLike_JCGAMMA.py @@ -160,7 +160,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, 'tree2.tre') + path = os.path.join(script_dir, 'treecopy.tre') with open(path, 'r') as content: newick = content.read() return newick @@ -308,6 +308,8 @@ def readnewick(tree): post = pre[:] post.reverse() + for nd in post: + print nd.number, nd.edgelen return post def Makenewick(pre):