[galaxy-dev] [hg] galaxy 1729: Corrected a bug in branch lengths batchfile, a...

[Greg Von Kuster]

details:   http://www.bx.psu.edu/hg/galaxy/rev/743237dd0e65
changeset: 1729:743237dd0e65
user:      guru
date:      Mon Feb 02 18:57:07 2009 -0500
description:
Corrected a bug in branch lengths batchfile, as per Sergei's suggestions.

1 file(s) affected in this change:

lib/galaxy/tools/util/hyphy_util.py

diffs (306 lines):

diff -r cabfdc9dea12 -r 743237dd0e65 lib/galaxy/tools/util/hyphy_util.py
--- a/lib/galaxy/tools/util/hyphy_util.py	Fri Jan 30 11:37:51 2009 -0500
+++ b/lib/galaxy/tools/util/hyphy_util.py	Mon Feb 02 18:57:07 2009 -0500
@@ -368,143 +368,281 @@
 
 BranchLengthsMF = """
 VERBOSITY_LEVEL            = -1;
+
 fscanf                  (PROMPT_FOR_FILE, "Lines", inLines);
+
+
 
 _linesIn = Columns (inLines);
 
+
+
 /*---------------------------------------------------------*/
 
+
+
 _currentGene   = 1; 
+
 _currentState = 0;
+
 geneSeqs      = "";
+
 geneSeqs      * 128;
 
+
+
 for (l=0; l<_linesIn; l=l+1)
+
 {
+
     if (Abs(inLines[l]) == 0)
+
     {
+
         if (_currentState == 1)
+
         {
+
             geneSeqs      * 0;
+
             DataSet       ds            = ReadFromString (geneSeqs);
+
             _processAGene (_currentGene);
-            geneSeqs * 128;
+
+             geneSeqs * 128;
+
             _currentGene = _currentGene + 1;
+
         }
+
     }
+
     else
+
     {
+
         if (_currentState == 0)
+
         {
+
             _currentState = 1;
+
         }
+
         geneSeqs * inLines[l];
+
         geneSeqs * "\\n";
+
     }
+
 }
 
+
+
 if (_currentState == 1)
+
 {
+
     geneSeqs      * 0;
+
     if (Abs(geneSeqs))
+
     {
+
         DataSet       ds            = ReadFromString (geneSeqs);
+
         _processAGene (_currentGene);
+
     }
+
 }
+
+
 
 fprintf (resultFile,CLOSE_FILE);
 
+
+
 /*---------------------------------------------------------*/
 
+
+
 function _processAGene (_geneID)
+
 {
+
     DataSetFilter             filteredData = CreateFilter (ds,1);
+
     if (_currentGene == 1)
+
     {
+
         SelectTemplateModel        (filteredData);
+
         
-        SetDialogPrompt         ("Tree file");
+
+        SetDialogPrompt           ("Tree file");
+
         fscanf                    (PROMPT_FOR_FILE, "Tree",  givenTree);
+
         fscanf                    (stdin, "String", resultFile);
+
         
+
         /* do sequence to branch map */
+
         
+
         validNames = {};
+
         taxonNameMap = {};
+
         
+
         for (k=0; k<TipCount(givenTree); k=k+1)
+
         {
+
             validNames[TipName(givenTree,k)&&1] = 1;
+
         }
+
         
+
         for (k=0; k<BranchCount(givenTree); k=k+1)
+
         {
+
             thisName = BranchName(givenTree,k);
+
             taxonNameMap[thisName&&1] = thisName;
+
         }
+
         
+
         storeValidNames = validNames;
+
         fprintf         (resultFile,CLEAR_FILE,KEEP_OPEN,"Block\\tBranch\\tLength\\tLowerBound\\tUpperBound\\n"); 
+
     }
+
     else
+
     {
+
+        HarvestFrequencies (vectorOfFrequencies, filteredData, 1,1,1);
+
         validNames = storeValidNames;
+
     }
+
     
+
     for (k=0; k<ds.species; k=k+1)
+
     {
+
         GetString (thisName, ds,k);
+
         shortName = (thisName^{{"\\\\..+",""}})&&1;
+
         if (validNames[shortName])
+
         {
+
             taxonNameMap[shortName] = thisName;
+
             validNames - (shortName);
+
             SetParameter (ds,k,shortName);
+
         }
+
         else
+
         {
+
             fprintf         (resultFile,"ERROR:", thisName, " could not be matched to any of the leaves in tree ", givenTree,"\\n"); 
+
             return 0;
+
         }
+
     }
+
     
+
     /* */
+
     
+
     LikelihoodFunction lf = (filteredData,givenTree);
+
+    Optimize                 (res,lf);
+
     
-    Optimize                 (res,lf);
+
+    timer = Time(0)-timer;
+
     
-    timer = Time(0)-timer;
+
+    branchNames   = BranchName   (givenTree,-1);
+
+    branchLengths = BranchLength (givenTree,-1);
+
     
-    branchNames   = BranchName   (givenTree,-1);
-    branchLengths = BranchLength (givenTree,-1);
+
     
+
+    for (k=0; k<Columns(branchNames)-1; k=k+1)
+
+    {
+
+        COVARIANCE_PARAMETER = "givenTree."+branchNames[k]+".t";
+
+        COVARIANCE_PRECISION = 0.95;
+
+        CovarianceMatrix (cmx,lf);
+
+        if (k==0)
+
+        {
+
+            /* compute a scaling factor */
+
+            ExecuteCommands ("givenTree."+branchNames[0]+".t=1");
+
+            scaleFactor = BranchLength (givenTree,0);
+
+            ExecuteCommands ("givenTree."+branchNames[0]+".t="+cmx[0][1]);
+
+        }
+
+        fprintf (resultFile,_geneID,"\\t",taxonNameMap[branchNames[k]&&1],"\\t",branchLengths[k],"\\t",scaleFactor*cmx[0][0],"\\t",scaleFactor*cmx[0][2],"\\n");
+
+    }
+
     
-    for (k=0; k<Columns(branchNames)-1; k=k+1)
-    {
-        COVARIANCE_PARAMETER = "givenTree."+branchNames[k]+".t";
-        COVARIANCE_PRECISION = 0.95;
-        CovarianceMatrix (cmx,lf);
-        if (k==0)
-        {
-            /* compute a scaling factor */
-            ExecuteCommands ("givenTree."+branchNames[0]+".t=1");
-            scaleFactor = BranchLength (givenTree,0);
-            ExecuteCommands ("givenTree."+branchNames[0]+".t="+cmx[0][1]);
-        }
-        fprintf (resultFile,_geneID,"\\t",taxonNameMap[branchNames[k]&&1],"\\t",branchLengths[k],"\\t",scaleFactor*cmx[0][0],"\\t",scaleFactor*cmx[0][2],"\\n");
-    }
-    
+
     ttl = (branchLengths*(Transpose(branchLengths["1"])))[0];
+
     global treeScaler = 1;
+
     ReplicateConstraint ("this1.?.t:=treeScaler*this2.?.t__",givenTree,givenTree);
+
     COVARIANCE_PARAMETER = "treeScaler";
+
     COVARIANCE_PRECISION = 0.95;
+
     CovarianceMatrix (cmx,lf);
+
     fprintf (resultFile,_geneID,"\\tTotal Tree\\t",ttl,"\\t",ttl*cmx[0][0],"\\t",ttl*cmx[0][2],"\\n");
+
     ClearConstraints (givenTree);
+
     return 0;
+
 }
 """