Cymon working

.buildconfig

@@ -0,0 +1,10 @@
+[default]
+name=Default
+runtime=host
+config-opts=
+run-opts=
+prefix=/home/cymon/.cache/gnome-builder/install/p4-phylogenetics.git.cymon/host
+app-id=
+postbuild=
+prebuild=
+default=true

.gitignore

@@ -29,4 +29,5 @@ build/
 # Sphinx documentation
 share/sphinxdoc/_build/
 
-
+make_pf.sh
+setup.py

Pf/defines.h

@@ -29,6 +29,11 @@
 #define RMATRIX_STMTREV 118
 #define RMATRIX_VT 119
 #define RMATRIX_PMB 120
+#define RMATRIX_PRASREV 121
+#define RMATRIX_GNETREV 122
+
+#define PIVEC_NONE	0
+#define PIVEC_TUPLE	1
 
 #define GAP_CODE       -1
 #define QMARK_CODE     -2

Pf/p4_model.c

@@ -384,7 +384,6 @@ void p4_newRMatrix(p4_model *aModel, int pNum, int mNum, int free, int spec)
     else if(spec == RMATRIX_RTREV) {
         rtRevRMatrix(aRMatrix->bigR);
     }
-
     else if(spec == RMATRIX_TMJTT94) {
         tmjtt94RMatrix(aRMatrix->bigR);
     }
@@ -412,6 +411,12 @@ void p4_newRMatrix(p4_model *aModel, int pNum, int mNum, int free, int spec)
     else if(spec == RMATRIX_STMTREV) {
         stmtREVRMatrix(aRMatrix->bigR);
     }
+    else if(spec == RMATRIX_PRASREV) {
+        prasREVRMatrix(aRMatrix->bigR);
+    }
+    else if(spec == RMATRIX_GNETREV) {
+        gnetREVRMatrix(aRMatrix->bigR);
+    }
     else if(spec == RMATRIX_VT) {
         vtRMatrix(aRMatrix->bigR);
     }

Pf/proteinModels.h

@@ -16,3 +16,5 @@ void gcpREVRMatrix(double **r);
 void stmtREVRMatrix(double **r);
 void vtRMatrix(double **r);
 void pmbRMatrix(double **r);
+void prasREVRMatrix(double **r);
+void gnetREVRMatrix(double **r);

bin/p4

@@ -38,9 +38,11 @@ argvAfterDoubleDash = []
 # back to being an empty string, the way it usually is when you run interactive
 # python.
 
-# Cymon seems to think this is needed to pull in the correct libraries:
-# sys.path = [os.path.dirname(sys.path[0])] + sys.path
-sys.path[0] = ''
+# If you have multiple p4 repositories you will need this path modification
+# to pull in the correct p4 libs: ensure that the p4 script is being executed
+# from the <repo>/bin directory 
+sys.path = [os.path.dirname(sys.path[0])] + sys.path
+#sys.path[0] = ''
 
 if len(sys.argv) > 1:  # Is there stuff on the command line?
     for f in sys.argv[1:]:
@@ -109,7 +111,6 @@ if doDrawTreesWithNodeNumbers:
     for t in var.trees:
         t.draw(showInternalNodeNames=1, addToBrLen=0.2, width=None, showNodeNums=1, partNum=None, model=None)
 
-
 if forceExitAtEnd:
     sys.exit()
 if exitAtEnd:
@@ -127,7 +128,7 @@ os.environ['PYTHONINSPECT'] = '1'
 # Set the prompt (default is '>>> ').  I suspect that doing this might
 # screw up using the python interpreter in python-mode in emacs, but
 # you'll want to confirm that.
-sys.ps1 = 'p4> '
+sys.ps1 = 'p4-gitc> '
 
 # When the "from p4 import *" was done above, it read in any *.py
 # files in ~/.p4.  Now we want to execfile one more file, which of
@@ -153,12 +154,10 @@ if var.interactiveHelper:
     elif var.interactiveHelper == 'ipython':
         from IPython import start_ipython
         if var.excepthookEditor:
-            sys.exit(start_ipython(argv=["-m", "p4.interactive.ipython", "-i"], 
+            sys.exit(start_ipython(argv=["-m", "p4.interactive.ipython", "-i"],
                                    user_ns=locals()))
         else:
             sys.exit(start_ipython(argv=["-i"], user_ns=locals()))
-
-
 
 # Maybe the user just types 'p4' and wants to see a splash...
 if len(sys.argv) == 1:

p4/__init__.py

@@ -85,7 +85,7 @@
 
 if 1:
     verboseStartupFiles = False  # Turn on for debugging...
-    if 0:  # If you want it, turn it on.
+    if 1:  # If you want it, turn it on.
         try:
             import __main__
             exec(open(os.environ['P4_STARTUP']).read())     # an official python anti-idiom!

p4/alignment_manip.py

@@ -1266,16 +1266,23 @@ def recodeDayhoff(self, firstLetter=False, symbols="123456"):
                     pass  # They stay as they are.
                 elif c == 'x':
                     s.sequence[i] = '-'
+                elif c in 'bzj':
+                    if ambigsBecomeGaps:
+                        s.sequence[i] = '-'
+                    else:
+                        gm.append("Ambiguity character '%s' is not handled." % c)
+                        gm.append("Perhaps set 'ambigsBecomeGaps' as True.")
+                        raise P4Error(gm)
                 else:
-                    # Maybe this should raise a P4Error?
-                    print("skipping character '%s'" % c)
-            s.sequence = ''.join(s.sequence)
+                    gm.append("Unknown character state '%s'" % c)
+                    raise P4Error(gm)
+            s.sequence = string.join(s.sequence, '')
         self.dataType = 'standard'
         self.equates = {}
         self.dim = 6
         self.symbols = symbols
 
-    def recodeProteinIntoGroups(self, groups, firstLetter=False):
+    def recodeProteinIntoGroups(self, groups, firstLetter=False, startAtZero=False):
         """Recode protein data into user-specified groups, in place.
 
         A generalization of :meth:`p4.alignment.Alignment.recodeDayhoff`
@@ -1288,7 +1295,10 @@ def recodeProteinIntoGroups(self, groups, firstLetter=False):
         It does not make a new alignment-- it does the re-coding 'in-place'.
 
         If arg *firstLetter* is set, then the character is recoded as the
-        first letter of its group rather than as a number. 
+        first letter of its group rather than as a number.
+
+        If arg *startAtZero* groups will be coded with number beginning with
+        zero, else beginning with 1.
         """
 
         gm = ['Alignment.recodeProteinIntoGroups()']
@@ -1308,7 +1318,10 @@ def recodeProteinIntoGroups(self, groups, firstLetter=False):
         for s in theseSymbols:
             assert s in self.symbols
             assert theseSymbols.count(s) == 1
-        numeralSymbols = ['%i' % (i + 1) for i in range(nGroups)]
+        if startAtZero:
+            numeralSymbols = ['%i' % i for i in range(nGroups)]
+        else:
+            numeralSymbols = ['%i' % (i + 1) for i in range(nGroups)]
         firstLetters = [gr[0] for gr in myGroups]
 
         for s in self.sequences:
@@ -2653,6 +2666,8 @@ def getMinmaxChiSqGroups(self, percent_cutoff=0.05, min_bins=2, max_bins=20,
                              n_choices=1000, seed=42, verbose=False):
         """An interface for Susko and Roger's minmax-chisq program.
 
+        This returns the maximum number of bins (groups of amino-acids) that maintains composition homogeneity.
+
         Susko, E. and Roger, A.J. (2007). On reduced amino acid alphabets for
         phylogenetic inference.  Mol. Biol. Evol. 24:2139-2150.
 
@@ -2706,31 +2721,56 @@ def getMinmaxChiSqGroups(self, percent_cutoff=0.05, min_bins=2, max_bins=20,
         # Need the binning before pvalue falls below percent_cutoff
         # They could be all <= 0.05 or all >= 0.05 because of the binning range
         pvalues = [float(bin[0].split()[1]) for bin in results]
-        if not pvalues[0] >= percent_cutoff:
-            print("No p-value <= %s" % percent_cutoff)
-            return None
         if not any(pv for pv in pvalues if pv >= percent_cutoff):
             print("No p-value >= %s" % percent_cutoff)
+            print("\nminmax-chisq output:\n%s" % stdout)
             return None
+        if not any(pv for pv in pvalues if pv <= percent_cutoff):
+            print("No p-value <= %s" % percent_cutoff)
+            print("\nminmax-chisq output:\n%s" % stdout)
+            return None
+        opt_bin = False
+        #Find the maximum number of groups that maintains homogeneity
+        #Loop through p-values and find the first value that is less than the
+        #cutoff, then pick the previous bin
+        #Find the first homogeneous bin
+        found_homogeneous = False
         for i, bin_result in enumerate(results):
             nbins, pvalue = bin_result[0].split()
-            if float(pvalue) <= percent_cutoff:
-                opt_bin = results[i - 1]
-                break
-        nbins, pvalue = opt_bin[0].split()
-        scores = opt_bin[1].split()
-        amino_acid_order = "A R N D C Q E G H I L K M F P S T W Y V".lower().split()
-        c = zip(scores, amino_acid_order)
-        groups = {}
-        for bin, amino in c:
-            sbin = str(bin)
-            groups[sbin] = groups.get(sbin, "") + amino
-        if verbose:
+            if verbose:
+                print("Doing bin %s, p-value = %s" % (bin_result, pvalue))
+            if not found_homogeneous:
+                if float(pvalue) >= percent_cutoff:
+                    found_homogeneous = True
+                    if verbose:
+                        print("\tFound Homogen %s, p-value = %s" % (bin_result, pvalue))
+                    continue
+            else:
+                #Find the next result < cuttoff and pick the previous bin
+                if float(pvalue) <= percent_cutoff:
+                    opt_bin = results[i-1]
+                    break
+        if not opt_bin:
+            print("Error: unable to find an optimal bin")
+            print("- try increasing number of bins, or lowering the %% cutoff")
             print("\nminmax-chisq output:\n%s" % stdout)
-            print("\nMaximum number of bins that maintains homogeneity: %s" % nbins)
-            print("\nGroups: %s" % ", ".join(groups.values()))
-
-        return groups.values()
+            return
+        else:
+            nbins, pvalue = opt_bin[0].split()
+            scores = opt_bin[1].split()
+            amino_acid_order = "A R N D C Q E G H I L K M F P S T W Y V".lower().split()
+            c = zip(scores, amino_acid_order)
+            groups = {}
+            for bin, amino in c:
+                sbin = str(bin)
+                groups[sbin] = groups.get(sbin, "") + amino
+            if verbose:
+                print("\nminmax-chisq output:\n%s" % stdout)
+                print("\nMaximum number of bins that maintains homogeneity: %s" % nbins)
+                print("\nGroups: %s" % ", ".join(groups.values()))
+            r = groups.values()
+            r.sort()
+            return r
 
     def getKosiolAISGroups(self, tree, n_bins, remove_files=False, verbose=True):
         """An interface for Kosiol's program AIS, for grouping amino acids.

p4/pnumbers.py

@@ -41,7 +41,7 @@ class Numbers(object):
 
     """
 
-    def __init__(self, inThing, col=0, skip=0):
+    def __init__(self, inThing, col=0, stop=0, skip=0):
         self.data = []
         self.bins = None
         # self.binSize is a property
@@ -55,7 +55,7 @@ def __init__(self, inThing, col=0, skip=0):
         if isinstance(inThing, numpy.ndarray):
             inThing = list(inThing)
         if inThing:
-            self.read(inThing, col, skip)
+            self.read(inThing, col, stop, skip)
 
     @property
     def binSize(self):
@@ -76,7 +76,7 @@ def binSize(self, binSize):
     def binSize(self):
         self._binSize = None
 
-    def read(self, inThing, col=0, skip=0):
+    def read(self, inThing, col=0, stop=0, skip=0):
         """Slurp in some more numbers.
 
         You can use this repeatedly, eg if your numbers are in more
@@ -88,6 +88,7 @@ def read(self, inThing, col=0, skip=0):
             try:
                 self.col = int(col)
                 self.skip = int(skip)
+                self.stop = int(stop)
             except (ValueError, TypeError):
                 gm.append("Args col and skip must be ints")
                 raise P4Error(gm)
@@ -100,7 +101,7 @@ def read(self, inThing, col=0, skip=0):
             #        gm.append("File '%s' has %i lines, " % (inThing, len(theLines)))
             #        gm.append("but skip is set to %i." % self.skip)
             #        raise P4Error(gm)
-            skipsDone = 0
+            count = 0
             digitsPlusMinus = '0123456789+-'
             for aLine in theLines:
                 ll = aLine.lstrip()
@@ -116,8 +117,11 @@ def read(self, inThing, col=0, skip=0):
                 elif ll[0] not in digitsPlusMinus:
                     pass
                 else:
-                    if skipsDone < self.skip:
-                        skipsDone += 1
+                    if count < self.skip:
+                        count += 1
+                        continue
+                    elif self.stop > 0 and count >= self.stop:
+                        break
                     else:
                         splitLine = aLine.split()
                         try:
@@ -134,6 +138,8 @@ def read(self, inThing, col=0, skip=0):
                             gm.append("Line '%s'.  " % aLine.rstrip())
                             gm.append("Can't make sense of '%s'" % theOne)
                             raise P4Error(gm)
+                        count += 1
+
         elif isinstance(inThing, list):
             for thing in inThing:
                 try:

p4/tree.py

@@ -1539,6 +1539,23 @@ def getSeqNumsAbove(self, nodeSpecifier):
                 seqNums.append(n.seqNum)
         return seqNums
 
+    def treeHeight(self):
+        """Mean, (min, max) of all leaf to root branch lengths"""
+        theights = self.getLeafToRootHeights()
+        return (p4.func.mean(theights), (min(theights), max(theights)))
+
+    def getLeafToRootHeights(self):
+        """Get all leaf to root branch lengths"""
+        theights = []
+        for nodeNum in self.iterLeavesNoRoot():
+            nodeHeight = 0
+            while self.node(nodeNum).parent:
+                if nodeNum != self.root.nodeNum:
+                    nodeHeight += self.node(nodeNum).br.len
+                nodeNum = self.node(nodeNum).parent.nodeNum
+            theights.append(nodeHeight)
+        return theights
+
     # def getAllChildrenNums(self, specifier):
     #     """Returns a list of the nodeNums of all children of the specified node
     #     Ambiguous, unused.

setup.py

@@ -44,8 +44,8 @@
 #my_include_dirs = ["/my/weird/include"]
 #my_lib_dirs = ["/my/weird/lib"]
 
-my_include_dirs = []
-my_lib_dirs = []
+my_include_dirs = ["/share/apps/include"]
+my_lib_dirs = ["/usr/local/lib64", "/share/apps/lib64"]
 
 likelyDirs = [ "/usr",
                "/usr/local",