some modifications to avoid redundancy

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)

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()

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())
+    result =  prepareTree(postorder, readSeq.patterns(sequence_file))

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)