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; + } """