diff --git a/res/TemplateBatchFiles/GARD.bf b/res/TemplateBatchFiles/GARD.bf index 8bdd67f74..176ebaf09 100644 --- a/res/TemplateBatchFiles/GARD.bf +++ b/res/TemplateBatchFiles/GARD.bf @@ -1,115 +1,148 @@ -RequireVersion("0.9920061128"); +RequireVersion("2.3"); + +LoadFunctionLibrary("libv3/UtilityFunctions.bf"); +LoadFunctionLibrary("libv3/IOFunctions.bf"); +LoadFunctionLibrary("libv3/all-terms.bf"); if (MPI_NODE_COUNT <= 1) { - fprintf (stdout, "[ERROR] This analysis requires an MPI environment to run\n"); - return 0; + fprintf (stdout, "[ERROR] This analysis requires an MPI environment to run\n"); + return 0; } -RequireVersion (".9920060901"); -partCount = 2; - -produceOffspring = MPI_NODE_COUNT-1; -populationSize = 2*produceOffspring; -incestDistance = 0; -generationCount = 5000; -maxSampleTries = populationSize*10; -mutationThreshhold = 0.0001; -mutationProb = 0.15; -mutationProbDecrease = 0.95; -annealingPhase = 100; -SHORT_MPI_RETURN = 1; -totalSampleCounter = 0; -localMutationRate = 0.05; -localMutationInterval = 20; +partCount = 2; +produceOffspring = MPI_NODE_COUNT-1; +populationSize = 2*produceOffspring; +incestDistance = 0; +generationCount = 5000; +maxSampleTries = populationSize*10; +mutationThreshhold = 0.0001; +mutationProb = 0.15; +mutationProbDecrease = 0.95; +annealingPhase = 100; +SHORT_MPI_RETURN = 1; +totalSampleCounter = 0; +localMutationRate = 0.05; +localMutationInterval = 20; shortestAllowedSegment = 0; -stoppingCriterion = 50; -sampleCount = 0; -familyControlSize = produceOffspring$6; +stoppingCriterion = 50; +sampleCount = 0; +familyControlSize = produceOffspring$6; -verboseFlag = 0; -rateClassesCount = 2; -MESSAGE_LOGGING = 0; -cAICPenaltyPerSite = 100; -adjustAICScores = 1; -maxTimeAllowed = 3600*1000; -startAtBreakpoint = 1; +verboseFlag = 0; +rateClassesCount = 2; +MESSAGE_LOGGING = 0; +cAICPenaltyPerSite = 100; +adjustAICScores = 1; +maxTimeAllowed = 3600*1000; +startAtBreakpoint = 1; /* ________________________________________________________________________________________________*/ -global AC = 1; -global AT = 1; -global CG = 1; -global CT = 1; -global GT = 1; +global AC = 1; +global AT = 1; +global CG = 1; +global CT = 1; +global GT = 1; /* ________________________________________________________________________________________________*/ -MasterList = {}; -REPLACE_TREE_STRUCTURE = 1; -bppMap = {}; -SHORT_MPI_RETURN = 1; -totalBitSize = 0; -LIKELIHOOD_FUNCTION_OUTPUT = 0; -FILE_SEPARATOR = "__FILE_SEPARATOR__"; +MasterList = {}; +REPLACE_TREE_STRUCTURE = 1; +bppMap = {}; +SHORT_MPI_RETURN = 1; +totalBitSize = 0; +LIKELIHOOD_FUNCTION_OUTPUT = 0; +FILE_SEPARATOR = "__FILE_SEPARATOR__"; ExecuteAFile (HYPHY_LIB_DIRECTORY+"TemplateBatchFiles"+DIRECTORY_SEPARATOR+"Utility"+DIRECTORY_SEPARATOR+"NJ.bf"); /* ________________________________________________________________________________________________*/ +function _reportSubstitutionMatrix (treeName, filterName) +{ + + ExecuteCommands ("bl = BranchLength (`treeName`,-1);"); + ExecuteCommands ("bn = BranchName (`treeName`,-1);"); + ExecuteCommands ("_size = BranchCount(`treeName`) + TipCount (`treeName`) - 1;"); + + for (k = _size; k >= 0 ; k = k - 1) + { + if (bl[k] > 0 || k == 0) + { + c = 1; + ExecuteCommands ("GetInformation(branchMx,`treeName`." + bn[k] +");"); + if (bl[k]) + { + normalizer = bl[k]; + } + else + { + normalizer = 1; + } + + rateMatrix = branchMx * (1/ normalizer); + break; + } + } + + return rateMatrix; + +} + function ComputeDistanceFormula (s1,s2) { - GetDataInfo (siteDifferenceCount, filteredData, s1, s2, DIST); - - totalSitesCompared = Transpose(summingVector)*(siteDifferenceCount*summingVector); - totalSitesCompared = totalSitesCompared[0]; - - if (_useK2P) - { - _dTransitionCounts = siteDifferenceCount[0][2]+siteDifferenceCount[2][0] /* A-G and G-A */ - +siteDifferenceCount[1][3]+siteDifferenceCount[3][1]; /* C-T and T-C */ - - _dTransversionCounts = (siteDifferenceCount[0][0]+siteDifferenceCount[1][1]+siteDifferenceCount[2][2]+siteDifferenceCount[3][3])+_dTransitionCounts; - - _dTransitionCounts = _dTransitionCounts/totalSitesCompared; - _dTransversionCounts = 1-_dTransversionCounts/totalSitesCompared; - - _d1C = 1-2*_dTransitionCounts-_dTransversionCounts; - _d2C = 1-2*_dTransversionCounts; - - if (_d1C>0 && _d2C>0) - { - return -(0.5*Log(_d1C)+.25*Log(_d2C)); - } - } - else - { - _dAGCounts = siteDifferenceCount[0][2]+siteDifferenceCount[2][0] /* A-G and G-A */; - _dCTCounts = siteDifferenceCount[1][3]+siteDifferenceCount[3][1]; /* C-T and T-C */ - - _dTransversionCounts = (siteDifferenceCount[0][0]+siteDifferenceCount[1][1]+siteDifferenceCount[2][2]+ - siteDifferenceCount[3][3])+_dAGCounts+_dCTCounts; - - _dAGCounts = _dAGCounts/totalSitesCompared; - _dCTCounts = _dCTCounts/totalSitesCompared; - - _dTransversionCounts = 1-_dTransversionCounts/totalSitesCompared; - - _d1C = 1-_dAGCounts/_d_TN_K1-0.5*_dTransversionCounts/_d_fR; - _d2C = 1-_dCTCounts/_d_TN_K2-0.5*_dTransversionCounts/_d_fY; - _d3C = 1-0.5*_dTransversionCounts/_d_fY/_d_fR; - - if ((_d1C>0)&&(_d2C>0)&&(_d3C>0)) - { - return -_d_TN_K1*Log(_d1C)-_d_TN_K2*Log(_d2C)-_d_TN_K3*Log(_d3C); - } - } - - return 1000; + GetDataInfo (siteDifferenceCount, filteredData, s1, s2, DIST); + + totalSitesCompared = Transpose(summingVector)*(siteDifferenceCount*summingVector); + totalSitesCompared = totalSitesCompared[0]; + + if (_useK2P) + { + _dTransitionCounts = siteDifferenceCount[0][2]+siteDifferenceCount[2][0] /* A-G and G-A */ + +siteDifferenceCount[1][3]+siteDifferenceCount[3][1]; /* C-T and T-C */ + + _dTransversionCounts = (siteDifferenceCount[0][0]+siteDifferenceCount[1][1]+siteDifferenceCount[2][2]+siteDifferenceCount[3][3])+_dTransitionCounts; + + _dTransitionCounts = _dTransitionCounts/totalSitesCompared; + _dTransversionCounts = 1-_dTransversionCounts/totalSitesCompared; + + _d1C = 1-2*_dTransitionCounts-_dTransversionCounts; + _d2C = 1-2*_dTransversionCounts; + + if (_d1C>0 && _d2C>0) + { + return -(0.5*Log(_d1C)+.25*Log(_d2C)); + } + } + else + { + _dAGCounts = siteDifferenceCount[0][2]+siteDifferenceCount[2][0] /* A-G and G-A */; + _dCTCounts = siteDifferenceCount[1][3]+siteDifferenceCount[3][1]; /* C-T and T-C */ + + _dTransversionCounts = (siteDifferenceCount[0][0]+siteDifferenceCount[1][1]+siteDifferenceCount[2][2]+ + siteDifferenceCount[3][3])+_dAGCounts+_dCTCounts; + + _dAGCounts = _dAGCounts/totalSitesCompared; + _dCTCounts = _dCTCounts/totalSitesCompared; + + _dTransversionCounts = 1-_dTransversionCounts/totalSitesCompared; + + _d1C = 1-_dAGCounts/_d_TN_K1-0.5*_dTransversionCounts/_d_fR; + _d2C = 1-_dCTCounts/_d_TN_K2-0.5*_dTransversionCounts/_d_fY; + _d3C = 1-0.5*_dTransversionCounts/_d_fY/_d_fR; + + if ((_d1C>0)&&(_d2C>0)&&(_d3C>0)) + { + return -_d_TN_K1*Log(_d1C)-_d_TN_K2*Log(_d2C)-_d_TN_K3*Log(_d3C); + } + } + + return 1000; } @@ -117,179 +150,176 @@ function ComputeDistanceFormula (s1,s2) function StringToMatrix (zz) { - return zz; + return zz; } /*---------------------------------------------------------------------------------------------------------------------------------------------*/ function ExportAMatrix (fileName, rateMatrix,dummy) { - if (Abs(rateMatrix)) - { - sortedBP = ConvertToPart (rateMatrix); - } - else - { - v = 0; - } - theAVL = {}; - - theAVL ["BP"] = {}; - theAVL ["Trees"] = {}; - - if (dataType) - { - bpF = 0; - bpF2 = 0; - } - else - { - bpF = -1; - bpF2 = -1; - } - - USE_DISTANCES = 0; - - lfDef = ""; - lfDef * 128; - lfDef * "LikelihoodFunction multiPart = ("; - - - partSpecs = {}; - - for (h=0; h1) - { - while (1) - { - for (nodeCounter = 0; nodeCounter < MPI_NODE_COUNT-1; nodeCounter = nodeCounter+1) - { - if (MPINodeState[nodeCounter][0]==1) - { - fromNode = ReceiveJobs (0,0); - break; - } - } - if (nodeCounter == MPI_NODE_COUNT-1) - { - break; - } - } - } - return 0; + if (MPI_NODE_COUNT>1) + { + while (1) + { + for (nodeCounter = 0; nodeCounter < MPI_NODE_COUNT-1; nodeCounter = nodeCounter+1) + { + if (MPINodeState[nodeCounter][0]==1) + { + fromNode = ReceiveJobs (0,0); + break; + } + } + if (nodeCounter == MPI_NODE_COUNT-1) + { + break; + } + } + } + return 0; } /*---------------------------------------------------------------------------------------------------------------------------------------------*/ function adjustAICScore (theLL,bpMatrix) { - daAICScore = 2*(baseParams*(baseSites/(baseSites-baseParams-1)) - theLL) ; - lastBpos = 0; - - for (_pid = 0; _pid < Rows(bpMatrix); _pid = _pid+1) - { - thisSpan = bppMap[bpMatrix[_pid]] - lastBpos+1; - lastBpos = bppMap[bpMatrix[_pid]]; - if (thisSpan > perPartParameterCount + 1) - { - daAICScore = daAICScore + 2*(perPartParameterCount*(thisSpan/(thisSpan-perPartParameterCount-1))); - } - else - { - daAICScore = daAICScore + 2*perPartParameterCount*cAICPenaltyPerSite; - } - } - - thisSpan = baseSites-lastBpos; - if (thisSpan > perPartParameterCount + 1) - { - daAICScore = daAICScore + 2*(perPartParameterCount*(thisSpan/(thisSpan-perPartParameterCount-1))); - } - else - { - daAICScore = daAICScore + 2*perPartParameterCount*cAICPenaltyPerSite; - } - - return -daAICScore; + daAICScore = 2*(baseParams*(baseSites/(baseSites-baseParams-1)) - theLL) ; + lastBpos = 0; + + for (_pid = 0; _pid < Rows(bpMatrix); _pid = _pid+1) + { + thisSpan = bppMap[bpMatrix[_pid]] - lastBpos+1; + lastBpos = bppMap[bpMatrix[_pid]]; + if (thisSpan > perPartParameterCount + 1) + { + daAICScore = daAICScore + 2*(perPartParameterCount*(thisSpan/(thisSpan-perPartParameterCount-1))); + } + else + { + daAICScore = daAICScore + 2*perPartParameterCount*cAICPenaltyPerSite; + } + } + + thisSpan = baseSites-lastBpos; + if (thisSpan > perPartParameterCount + 1) + { + daAICScore = daAICScore + 2*(perPartParameterCount*(thisSpan/(thisSpan-perPartParameterCount-1))); + } + else + { + daAICScore = daAICScore + 2*perPartParameterCount*cAICPenaltyPerSite; + } + + return -daAICScore; } @@ -297,120 +327,120 @@ function adjustAICScore (theLL,bpMatrix) function ReceiveJobs (sendOrNot, ji) { - if (MPI_NODE_COUNT>1) - { - MPIReceive (-1, fromNode, result_String); - mji = MPINodeState[fromNode-1][1]; - - if (sendOrNot) - { - MPISend (fromNode,mpiStringToSend); - MPINodeState[fromNode-1][1] = ji; - } - else - { - MPINodeState[fromNode-1][0] = 0; - MPINodeState[fromNode-1][1] = -1; - } - ExecuteCommands (result_String); - myDF = lf_MLES[1][1]+baseParams; - myAIC = 2*(lf_MLES[1][0]-myDF*(baseSites/(baseSites-myDF-1))); - myLL = lf_MLES[1][0]; - ji = mji; - } - else - { - myDF = res[1][1]+baseParams; - myAIC = 2*(res[1][0]-myDF*(baseSites/(baseSites-myDF-1))); - } - - sortedBP = {{-1}}; - - if (resultProcessingContext==0) - { - sortedScores[ji][0] = myAIC; - if (ji>=0) - { - jobPrint = ConvertToPartString (currentPopulation[ji]); - sortedBP = ConvertToPart (currentPopulation[ji]); - if (adjustAICScores) - { - myAIC = adjustAICScore (myLL, sortedBP); - } - v = Rows (sortedBP); - sortedScores[ji][0] = myAIC; - } - if (verboseFlag) - { - fprintf (stdout, "Individual ",ji," AIC-c = ",-myAIC," "); - } - } - else - { - intermediateProbs[ji][0] = myAIC; - if (ji>=0) - { - jobPrint = ConvertToPartString (children[ji-populationSize]); - sortedBP = ConvertToPart (children[ji-populationSize]); - if (adjustAICScores) - { - myAIC = adjustAICScore (myLL, sortedBP); - } - v = Rows (sortedBP); - intermediateProbs[ji][0] = myAIC; - } - if (verboseFlag) - { - fprintf (stdout, "Offspring ",ji," AIC-c = ",-myAIC," "); - } - } - - if (sortedBP[0]>=0) - { - if (dataType) - { - bpF = 0; - } - else - { - bpF = -1; - } - filterString = ""; - for (h=0; h1) + { + MPIReceive (-1, fromNode, result_String); + mji = MPINodeState[fromNode-1][1]; + + if (sendOrNot) + { + MPISend (fromNode,mpiStringToSend); + MPINodeState[fromNode-1][1] = ji; + } + else + { + MPINodeState[fromNode-1][0] = 0; + MPINodeState[fromNode-1][1] = -1; + } + ExecuteCommands (result_String); + myDF = lf_MLES[1][1]+baseParams; + myAIC = 2*(lf_MLES[1][0]-myDF*(baseSites/(baseSites-myDF-1))); + myLL = lf_MLES[1][0]; + ji = mji; + } + else + { + myDF = res[1][1]+baseParams; + myAIC = 2*(res[1][0]-myDF*(baseSites/(baseSites-myDF-1))); + } + + sortedBP = {{-1}}; + + if (resultProcessingContext==0) + { + sortedScores[ji][0] = myAIC; + if (ji>=0) + { + jobPrint = ConvertToPartString (currentPopulation[ji]); + sortedBP = ConvertToPart (currentPopulation[ji]); + if (adjustAICScores) + { + myAIC = adjustAICScore (myLL, sortedBP); + } + v = Rows (sortedBP); + sortedScores[ji][0] = myAIC; + } + if (verboseFlag) + { + fprintf (stdout, "Individual ",ji," AIC-c = ",-myAIC," "); + } + } + else + { + intermediateProbs[ji][0] = myAIC; + if (ji>=0) + { + jobPrint = ConvertToPartString (children[ji-populationSize]); + sortedBP = ConvertToPart (children[ji-populationSize]); + if (adjustAICScores) + { + myAIC = adjustAICScore (myLL, sortedBP); + } + v = Rows (sortedBP); + intermediateProbs[ji][0] = myAIC; + } + if (verboseFlag) + { + fprintf (stdout, "Offspring ",ji," AIC-c = ",-myAIC," "); + } + } + + if (sortedBP[0]>=0) + { + if (dataType) + { + bpF = 0; + } + else + { + bpF = -1; + } + filterString = ""; + for (h=0; h0) - { - minPartLength = curSegLength; - } - } - - if (bpF20) - { - minPartLength = curSegLength; - } - } - - _ConstraintString * 0; - return _ConstraintString; + sortedBP = ConvertToPart (pString); + if (dataType) + { + bpF = 0; + } + else + { + bpF = -1; + } + + minPartLength = 1e100; + + _ConstraintString = ""; + _ConstraintString * 256; + + + for (h=0; h0) + { + minPartLength = curSegLength; + } + } + + if (bpF20) + { + minPartLength = curSegLength; + } + } + + _ConstraintString * 0; + return _ConstraintString; } /*---------------------------------------------------------------------------------------------------------------------------------------------*/ function ConvertToPart (pString) { - partitionHits = {}; - h = 0; - v = bppSize; - - for (mpiNode=0; mpiNode=Abs(bppMap)-1) - { - aBP = aBP - 2^(bppSize-1); - } - - v = v + bppSize; - partitionHits[aBP] = 1; - } - - meKeys = Rows(partitionHits); - v = Columns(meKeys); - sortedBP = {v,1}; - - for (h=0; h=Abs(bppMap)-1) + { + aBP = aBP - 2^(bppSize-1); + } + + v = v + bppSize; + partitionHits[aBP] = 1; + } + + meKeys = Rows(partitionHits); + v = Columns(meKeys); + sortedBP = {v,1}; + + for (h=0; h1) && (jobIndex>=0)) - { - mpiStringToSend * ("Optimize(res,lf);return \"lf_MLES=\"+res+\";\";"); - mpiStringToSend * 0; - for (mpiNode = 0; mpiNode < MPI_NODE_COUNT-1; mpiNode = mpiNode+1) - { - if (MPINodeState[mpiNode][0]==0) - { - break; - } - } - if (mpiNode==MPI_NODE_COUNT-1) - { - mpiNode = ReceiveJobs (1,jobIndex); - } - else - { - MPISend (mpiNode+1,mpiStringToSend); - MPINodeState[mpiNode][0] = 1; - MPINodeState[mpiNode][1] = jobIndex; - } - } - else - { - mpiStringToSend * 0; - ExecuteCommands (mpiStringToSend); - Optimize (res,lf); - if (jobIndex>=0) - { - mpiNode = ReceiveJobs (1, jobIndex); - } - else - { - myAIC = 2*(res[1][0]-res[1][1]-baseParams); - } - } - return 0; +{ + myAIC = MasterList[ConvertToPartString (cString)]; + + if (myAIC<0) + { + if (resultProcessingContext==0) + { + sortedScores[jobIndex][0] = myAIC; + if (verboseFlag) + { + fprintf (stdout, "Individual ",jobIndex," AIC-c = ",-myAIC, "\n"); + } + } + else + { + intermediateProbs[jobIndex][0] = myAIC; + if (verboseFlag) + { + fprintf (stdout, "Offspring ",jobIndex," AIC-c = ",-myAIC,"\n"); + } + } + return 0; + } + + if (dataType) + { + bpF = 0; + } + else + { + bpF = -1; + } + + mpiStringToSend = ""; + mpiStringToSend * 8192; + + ConstraintString = ""; + LikelihoodFunctionString = ""; + ConstraintString * 8192; + LikelihoodFunctionString * 256; + LikelihoodFunctionString * "LikelihoodFunction lf=("; + + + for (h=0; h1) && (jobIndex>=0)) + { + mpiStringToSend * ("Optimize(res,lf);return \"lf_MLES=\"+res+\";\";"); + mpiStringToSend * 0; + for (mpiNode = 0; mpiNode < MPI_NODE_COUNT-1; mpiNode = mpiNode+1) + { + if (MPINodeState[mpiNode][0]==0) + { + break; + } + } + if (mpiNode==MPI_NODE_COUNT-1) + { + mpiNode = ReceiveJobs (1,jobIndex); + } + else + { + MPISend (mpiNode+1,mpiStringToSend); + MPINodeState[mpiNode][0] = 1; + MPINodeState[mpiNode][1] = jobIndex; + } + } + else + { + mpiStringToSend * 0; + ExecuteCommands (mpiStringToSend); + Optimize (res,lf); + if (jobIndex>=0) + { + mpiNode = ReceiveJobs (1, jobIndex); + } + else + { + myAIC = 2*(res[1][0]-res[1][1]-baseParams); + } + } + return 0; } @@ -687,81 +717,81 @@ function RunASample (dummy,jobIndex) function SpawnRandomString (clsCnt) { - rModel = {totalBitSize,1}; - for (h=0; h (-0.1); _lc = _lc + 1) - { - if (Rows (currentPopulation[_lc])) - { - myAIC = -(sampleString == ConvertToPartString (currentPopulation[_lc])); - } - } - for (_lc = 0; _lc < Abs(children) && myAIC > (-0.1); _lc = _lc + 1) - { - if (Rows (children[_lc])) - { - myAIC = -(sampleString == ConvertToPartString (children[_lc])); - } - } - - while ((myAIC<(-0.1) || minPartLength 1) - { - fprintf (stdout,"Adjusting the child to avoid a duplicate. Min(fragment) = ",minPartLength,". Pass ", mutPassCount, "\n"); - } - - mutPassCount = mutPassCount + 1; - sampleString = Min(Random(0,stateVectorDimension)$1,stateVectorDimension-1); - myAIC = testChild[sampleString]; - newValue = Random (0,rateClassesCount-0.0000001)$1; - - while (newValue == myAIC) - { - newValue = Random (0,rateClassesCount-0.0000001)$1; - } - - testChild [sampleString] = newValue; - sampleString = ConvertToPartString (testChild); - myAIC = MasterList [sampleString]; - for (_lc = 0; _lc < populationSize && myAIC > (-0.1); _lc = _lc + 1) - { - if (Rows (currentPopulation[_lc])) - { - myAIC = -(sampleString == ConvertToPartString (currentPopulation[_lc])); - } - } - for (_lc = 0; _lc < Abs(children) && myAIC > (-0.1); _lc = _lc + 1) - { - if (Rows (children[_lc])) - { - myAIC = -(sampleString == ConvertToPartString (children[_lc])); - } - } - } - return testChild; + sampleString = ConvertToPartString (putativeChild); + myAIC = MasterList[sampleString]; + testChild = putativeChild; + mutPassCount = 1; + + for (_lc = 0; _lc < populationSize && myAIC > (-0.1); _lc = _lc + 1) + { + if (Rows (currentPopulation[_lc])) + { + myAIC = -(sampleString == ConvertToPartString (currentPopulation[_lc])); + } + } + for (_lc = 0; _lc < Abs(children) && myAIC > (-0.1); _lc = _lc + 1) + { + if (Rows (children[_lc])) + { + myAIC = -(sampleString == ConvertToPartString (children[_lc])); + } + } + + while ((myAIC<(-0.1) || minPartLength 1) + { + fprintf (stdout,"Adjusting the child to avoid a duplicate. Min(fragment) = ",minPartLength,". Pass ", mutPassCount, "\n"); + } + + mutPassCount = mutPassCount + 1; + sampleString = Min(Random(0,stateVectorDimension)$1,stateVectorDimension-1); + myAIC = testChild[sampleString]; + newValue = Random (0,rateClassesCount-0.0000001)$1; + + while (newValue == myAIC) + { + newValue = Random (0,rateClassesCount-0.0000001)$1; + } + + testChild [sampleString] = newValue; + sampleString = ConvertToPartString (testChild); + myAIC = MasterList [sampleString]; + for (_lc = 0; _lc < populationSize && myAIC > (-0.1); _lc = _lc + 1) + { + if (Rows (currentPopulation[_lc])) + { + myAIC = -(sampleString == ConvertToPartString (currentPopulation[_lc])); + } + } + for (_lc = 0; _lc < Abs(children) && myAIC > (-0.1); _lc = _lc + 1) + { + if (Rows (children[_lc])) + { + myAIC = -(sampleString == ConvertToPartString (children[_lc])); + } + } + } + return testChild; } /*---------------------------------------------------------------------------------------------------------------------------------------------*/ function UpdateBL (dummy) { - return 0; + return 0; } @@ -770,51 +800,51 @@ function UpdateBL (dummy) dataType = 0; -fprintf (stdout, "Initialized GARD on ", MPI_NODE_COUNT, " MPI nodes.\nPopulation size is ", populationSize, " models\n"); -SetDialogPrompt ("Nucleotide file to screen:"); -DataSet ds = ReadDataFile (PROMPT_FOR_FILE); +fprintf (stdout, "Initialized GARD on ", MPI_NODE_COUNT, " MPI nodes.\nPopulation size is ", populationSize, " models\n"); +SetDialogPrompt ("Nucleotide file to screen:"); +DataSet ds = ReadDataFile (PROMPT_FOR_FILE); -baseFilePath = LAST_FILE_PATH; +baseFilePath = LAST_FILE_PATH; done = 0; while (!done) { - fprintf (stdout,"\nPlease enter a 6 character model designation (e.g:010010 defines HKY85):"); - fscanf (stdin,"String", modelDesc); - if (Abs(modelDesc)==6) - { - done = 1; - } + fprintf (stdout,"\nPlease enter a 6 character model designation (e.g:010010 defines HKY85):"); + fscanf (stdin,"String", modelDesc); + if (Abs(modelDesc)==6) + { + done = 1; + } } ChoiceList (rvChoice, "Rate variation options", 1, SKIP_NONE, - "None", "Homogeneous rates across sites (fastest)", - "General Discrete", "General discrete distribution on N-bins", - "Beta-Gamma", "Adaptively discretized gamma N-bins. Try is GDD doesn't converge well, or if you suspect a lot of rate classes and do not want to overparameterize the model"); - + "None", "Homogeneous rates across sites (fastest)", + "General Discrete", "General discrete distribution on N-bins", + "Beta-Gamma", "Adaptively discretized gamma N-bins. Try is GDD doesn't converge well, or if you suspect a lot of rate classes and do not want to overparameterize the model"); + if (rvChoice<0) { - return 0; + return 0; } if (rvChoice) { - rateClasses = 0; - while (rateClasses < 2 || rateClasses > 32) - { - fprintf (stdout, "How many distribution bins [2-32]?:"); - fscanf (stdin, "Number", rateClasses); - rateClasses = rateClasses $ 1; - } - fprintf (stdout, "\nUsing ", rateClasses, " distribution bins\n"); + rateClasses = 0; + while (rateClasses < 2 || rateClasses > 32) + { + fprintf (stdout, "How many distribution bins [2-32]?:"); + fscanf (stdin, "Number", rateClasses); + rateClasses = rateClasses $ 1; + } + fprintf (stdout, "\nUsing ", rateClasses, " distribution bins\n"); } -DEFAULT_FILE_SAVE_NAME = baseFilePath + ".html"; +DEFAULT_FILE_SAVE_NAME = baseFilePath + ".html"; -SetDialogPrompt ("Save results to:"); -fprintf (PROMPT_FOR_FILE, CLEAR_FILE, KEEP_OPEN); -outputFilePath = LAST_FILE_PATH; +SetDialogPrompt ("Save results to:"); +fprintf (PROMPT_FOR_FILE, CLEAR_FILE, KEEP_OPEN); +outputFilePath = LAST_FILE_PATH; @@ -822,30 +852,31 @@ DataSetFilter filteredData = CreateFilter (ds,1); InferTreeTopology (0); treeString = TreeMatrix2TreeString(0); +baseTreeString = treeString; /* find "informative sites" */ if (dataType==0) { - for (h=0; h0) - { - m1=m1+1; - } - } - if (m1>1) - { - bppMap[Abs(bppMap)] = h; - } - } + for (h=0; h0) + { + m1=m1+1; + } + } + if (m1>1) + { + bppMap[Abs(bppMap)] = h; + } + } } @@ -858,274 +889,274 @@ h = Abs(bppMap); if (h <= partCount) { - fprintf (outputFilePath, "ERROR: There are too few potential break points to support ", partCount-1, " recombination events.\n"); - return 0; + fprintf (outputFilePath, "ERROR: There are too few potential break points to support ", partCount-1, " recombination events.\n"); + return 0; } -maxBPP = h-1; +maxBPP = h-1; ModelTitle = ""+modelDesc[0]; - + rateBiasTerms = {{"AC","1","AT","CG","CT","GT"}}; -paramCount = 0; +paramCount = 0; modelConstraintString = ""; for (customLoopCounter2=1; customLoopCounter2<6; customLoopCounter2=customLoopCounter2+1) { - for (customLoopCounter=0; customLoopCounter 8) - { - rateClasses = 8; - } - } - - if (ds.species < 5) - { - rateClasses = Min (3,rateClasses); - } - else - { - if (ds.species < 10) - { - rateClasses = Min (4,rateClasses); - } - else - { - if (ds.species < 30) - { - rateClasses = Min (5,rateClasses); - } - } - } - - if (rvChoice == 1) - { - gdDefString = ""; - gdDefString * 1024; - for (mi=1; mi0.000000001;\n"); - - for (mi=3; mi<=rateClasses; mi=mi+1) - { - gdDefString*("global RS_"+mi+" = 1.5;"+"\nRS_"+mi+":>1;RS_"+mi+":<100000;\n"); - } - - rateStrMx = {rateClasses,1}; - rateStrMx[0] = "RS_1"; - rateStrMx[1] = "1"; - - for (mi=3; mi<=rateClasses; mi=mi+1) - { - rateStrMx[mi-1] = rateStrMx[mi-2]+"*RS_"+mi; - } - - freqStrMx = {rateClasses,1}; - freqStrMx[0] = "PS_1"; - - for (mi=1; mi0.05;betaP:<85; - betaQ:>0.05;betaQ:<85; - category pc = (rateClasses-1, EQUAL, MEAN, - _x_^(betaP-1)*(1-_x_)^(betaQ-1)/Beta(betaP,betaQ), /* density */ - IBeta(_x_,betaP,betaQ), /*CDF*/ - 0, /*left bound*/ - 1, /*right bound*/ - IBeta(_x_,betaP+1,betaQ)*betaP/(betaP+betaQ) - ); - - global alpha = .5; - alpha:>0.01;alpha:<100; - category c = (rateClasses, pc, MEAN, - GammaDist(_x_,alpha,alpha), - CGammaDist(_x_,alpha,alpha), - 0 , - 1e25, - CGammaDist(_x_,alpha+1,alpha) - ); - - } - NucleotideMatrix = {{*,c*AC*t,c*t,c*AT*t} - {c*AC*t,*,c*CG*t,c*CT*t} - {c*t,c*CG*t,*,c*GT*t} - {c*AT*t,c*CT*t,c*GT*t,*}}; + rateClasses = rateClasses $ 1; + if (rateClasses < 2) + { + rateClasses = 2; + } + else + { + if (rateClasses > 8) + { + rateClasses = 8; + } + } + + if (ds.species < 5) + { + rateClasses = Min (3,rateClasses); + } + else + { + if (ds.species < 10) + { + rateClasses = Min (4,rateClasses); + } + else + { + if (ds.species < 30) + { + rateClasses = Min (5,rateClasses); + } + } + } + + if (rvChoice == 1) + { + gdDefString = ""; + gdDefString * 1024; + for (mi=1; mi0.000000001;\n"); + + for (mi=3; mi<=rateClasses; mi=mi+1) + { + gdDefString*("global RS_"+mi+" = 1.5;"+"\nRS_"+mi+":>1;RS_"+mi+":<100000;\n"); + } + + rateStrMx = {rateClasses,1}; + rateStrMx[0] = "RS_1"; + rateStrMx[1] = "1"; + + for (mi=3; mi<=rateClasses; mi=mi+1) + { + rateStrMx[mi-1] = rateStrMx[mi-2]+"*RS_"+mi; + } + + freqStrMx = {rateClasses,1}; + freqStrMx[0] = "PS_1"; + + for (mi=1; mi0.05;betaP:<85; + betaQ:>0.05;betaQ:<85; + category pc = (rateClasses-1, EQUAL, MEAN, + _x_^(betaP-1)*(1-_x_)^(betaQ-1)/Beta(betaP,betaQ), /* density */ + IBeta(_x_,betaP,betaQ), /*CDF*/ + 0, /*left bound*/ + 1, /*right bound*/ + IBeta(_x_,betaP+1,betaQ)*betaP/(betaP+betaQ) + ); + + global alpha = .5; + alpha:>0.01;alpha:<100; + category c = (rateClasses, pc, MEAN, + GammaDist(_x_,alpha,alpha), + CGammaDist(_x_,alpha,alpha), + 0 , + 1e25, + CGammaDist(_x_,alpha+1,alpha) + ); + + } + NucleotideMatrix = {{*,c*AC*t,c*t,c*AT*t} + {c*AC*t,*,c*CG*t,c*CT*t} + {c*t,c*CG*t,*,c*GT*t} + {c*AT*t,c*CT*t,c*GT*t,*}}; } else { - NucleotideMatrix = {{*,AC*t,t,AT*t}{AC*t,*,CG*t,CT*t}{t,CG*t,*,GT*t}{AT*t,CT*t,GT*t,*}}; + NucleotideMatrix = {{*,AC*t,t,AT*t}{AC*t,*,CG*t,CT*t}{t,CG*t,*,GT*t}{AT*t,CT*t,GT*t,*}}; } -HarvestFrequencies (nucEFV, filteredData, 1, 1, 1); -Model nucModel = (NucleotideMatrix, nucEFV, 1); -Tree givenTree = treeString; -branchNames = BranchName (givenTree,-1); +HarvestFrequencies (nucEFV, filteredData, 1, 1, 1); +Model nucModel = (NucleotideMatrix, nucEFV, 1); +Tree givenTree = treeString; +branchNames = BranchName (givenTree,-1); LikelihoodFunction lf = (filteredData,givenTree); /* check parameter counts */ -GetString (varList, lf, -1); -baseParams = Columns(varList["Global Independent"])+3; -perPartParameterCount = Columns(varList["Local Independent"]); -baseSites = filteredData.sites; +GetString (varList, lf, -1); +baseParams = Columns(varList["Global Independent"])+3; +perPartParameterCount = Columns(varList["Local Independent"]); +baseSites = filteredData.sites; if (baseParams + (partCount+1) * perPartParameterCount >= baseSites - 1) { - fprintf (stdout, "ERROR: Too few sites for c-AIC inference.\n"); - return 0; + fprintf (stdout, "ERROR: Too few sites for c-AIC inference.\n"); + return 0; } if (baseParams + 2 * perPartParameterCount >= baseSites - 1) { - fprintf (stdout, "ERROR: Too few sites for reliable c-AIC inference.\n"); - return 0; + fprintf (stdout, "ERROR: Too few sites for reliable c-AIC inference.\n"); + return 0; } -fprintf (stdout, "\nFitting a baseline nucleotide model...\n"); -Optimize (res,lf); -baseLL = res[1][0]; +fprintf (stdout, "\nFitting a baseline nucleotide model...\n"); +Optimize (res,lf); +baseLL = res[1][0]; currentPopulation = {}; -sortedScores = {populationSize,2}; -/*baseParams = res[1][2];*/ +sortedScores = {populationSize,2}; +/*baseParams = res[1][2];*/ -myDF = baseParams+perPartParameterCount; +myDF = baseParams+perPartParameterCount; sortedScores[0][0] = 2*(res[1][0]-myDF*(baseSites/(baseSites-myDF-1))); sortedScores[0][1] = 0; fprintf (stdout, CLEAR_FILE, "Done with single partition analysis. ", - "Log(L) = ",lf, - ", c-AIC = ",-sortedScores[0][0], "\n"); + "Log(L) = ",lf, + ", c-AIC = ",-sortedScores[0][0], "\n"); if (baseParams>3) { - ConstraintString = ""; - ConstraintString*256; - for (h=0; h1) { - MPINodeState = {MPI_NODE_COUNT-1,3}; - MPINodeBreakpoints = {}; + MPINodeState = {MPI_NODE_COUNT-1,3}; + MPINodeBreakpoints = {}; } -currentBEST_IC = crapAIC; -ibfPath = "GARD_GA_CHC.ibf"; +currentBEST_IC = crapAIC; +ibfPath = "GARD_GA_CHC.ibf"; current_BPP_done_counter = 0; lastImprovedBPC = 0; byBPImprovement = {}; -byBPSplits = {}; +byBPSplits = {}; DataSetFilter allData = CreateFilter (ds, 1); @@ -1136,250 +1167,408 @@ fprintf (stdout, "Starting the GA...\n"); for (currentBPC = startAtBreakpoint; currentBPC < maxBPP; currentBPC = currentBPC + 1) { - partCount = currentBPC; - totalBitSize = bppSize * partCount; - stateVectorDimension = totalBitSize; - - if (currentBPC == startAtBreakpoint) - { - currentPopulation [0] = {totalBitSize,1}; - } - else - { - for (k=0; k kf || currentBPC == 1) - { - byBPSplits [currentBPC] = ConvertToPart (currentPopulation[populationSize-1]); - byBPImprovement [currentBPC] = currentBEST_IC-kf; - if (currentBEST_IC > kf) - { - lastImprovedBPC = currentBPC; - bestIndividualOverall = currentPopulation[populationSize-1]; - } - currentBEST_IC = Min(kf,currentBEST_IC); - } - else - { - break; - } + partCount = currentBPC; + totalBitSize = bppSize * partCount; + stateVectorDimension = totalBitSize; + + if (currentBPC == startAtBreakpoint) + { + currentPopulation [0] = {totalBitSize,1}; + } + else + { + for (k=0; k kf || currentBPC == 1) + { + byBPSplits [currentBPC] = ConvertToPart (currentPopulation[populationSize-1]); + byBPImprovement [currentBPC] = currentBEST_IC-kf; + if (currentBEST_IC > kf) + { + lastImprovedBPC = currentBPC; + bestIndividualOverall = currentPopulation[populationSize-1]; + } + currentBEST_IC = Min(kf,currentBEST_IC); + } + else + { + break; + } } fprintf (outputFilePath, "GARD Results"); fprintf (outputFilePath, "
Processed alignment ", baseFilePath," with ", ds.species, " sequences and ", ds.sites, " sites.", - "

Processor time taken: ", Time(1)-startTimer, " seconds.

"); - + "

Processor time taken: ", Time(1)-startTimer, " seconds."); + if (timeLeft < 0) { - fprintf (outputFilePath, "

This analysis was stopped before convergence, because the CPU time limit per job was reached. The reported results may, therefore be incomplete, as there may be more breakpoints, for example. ", - "Consider reducing the size of the alignment, either by removing sequences or selecting shorted sequence regions, or choosing a simpler model (no rate variation etc).
\n"); + fprintf (outputFilePath, "
This analysis was stopped before convergence, because the CPU time limit per job was reached. The reported results may, therefore be incomplete, as there may be more breakpoints, for example. ", + "Consider reducing the size of the alignment, either by removing sequences or selecting shorted sequence regions, or choosing a simpler model (no rate variation etc).
\n"); } fprintf (outputFilePath, "
Results

", reportImprovementHTML (0), "

"); - + fprintf (outputFilePath, "
Substitution process parameters

Nucleotide Model (",modelDesc, - ") Fit Results

Log likelihood : ",outputSpan (150,baseLL),"
", - "c-AIC : ", outputSpan(150, crapAIC),"
"); + ") Fit Results

Log likelihood : ",outputSpan (150,baseLL),"
", + "c-AIC : ", outputSpan(150, crapAIC),"
"); if (rvChoice) { - if (rvChoice == 1) - { - fprintf (outputFilePath, "

Using general discrete distribution of rates across sites"); - } - else - { - fprintf (outputFilePath, "

Using the beta-gamma distribution of rates across sites"); - } - GetInformation (cI,c); - for (k = 0; k < Columns (cI); k=k+1) - { - fprintf (outputFilePath, "
Rate : ", Format (cI[0][k],7,3), outputSpan (150, "Weight : " + Format (cI[1][k],7,3))); - } + if (rvChoice == 1) + { + fprintf (outputFilePath, "

Using general discrete distribution of rates across sites"); + } + else + { + fprintf (outputFilePath, "

Using the beta-gamma distribution of rates across sites"); + } + GetInformation (cI,c); + for (k = 0; k < Columns (cI); k=k+1) + { + fprintf (outputFilePath, "
Rate : ", Format (cI[0][k],7,3), outputSpan (150, "Weight : " + Format (cI[1][k],7,3))); + } } else { - fprintf (outputFilePath, "

Homogeneous rate distribution across sites"); + fprintf (outputFilePath, "

Homogeneous rate distribution across sites"); } fprintf (outputFilePath, "

Nucleotide substitution rate matrix

", - "", - "", - "", - "", - "
 ACGT
A*",AC,"1",AT,"
C-*",CG,"",CT,"
G--*",GT,"
T---*
", - "

"); + "", + "", + "", + "", + "
 ACGT
A*",AC,"1",AT,"
C-*",CG,"",CT,"
G--*",GT,"
T---*
", + ""); fprintf (outputFilePath, "", CLOSE_FILE); -fprintf (stdout, "Performing the final optimization...\n"); +fprintf (stdout, "Performing the final optimization...\n"); rawOutFile = outputFilePath + "_splits"; fprintf (rawOutFile,CLEAR_FILE); +// print rate matrix to new JSON file + +rateMatrix = _reportSubstitutionMatrix("givenTree", "filteredData"); +jsonOutFile = outputFilePath + ".json"; + +gard.json = { + "rateMatrix" : rateMatrix +}; +gard.json[utility.getGlobalValue("terms.json.input")] = {}; +(gard.json[utility.getGlobalValue("terms.json.input")])[utility.getGlobalValue("terms.json.file")] = baseFilePath; +(gard.json[utility.getGlobalValue("terms.json.input")])[utility.getGlobalValue("terms.json.sequences")] = Eval("filteredData.species"); +(gard.json[utility.getGlobalValue("terms.json.input")])[utility.getGlobalValue("terms.json.sites")] = filteredData.sites; +(gard.json[utility.getGlobalValue("terms.json.input")])[utility.getGlobalValue("terms.json.partition_count")] = partCount; +gard.json["timeElapsed"] = Time(1) - startTimer; +gard.json["rvChoice"] = rvChoice; +gard.json["rvRates"] = rateClasses; +gard.json["potentialBreakpoints"] = Abs(bppMap); if (Rows(bestIndividualOverall)) { - totalBitSize = Rows (bestIndividualOverall)*Columns(bestIndividualOverall); - partCount = totalBitSize/bppSize; - outAVL = ExportAMatrix (rawOutFile,bestIndividualOverall,0); + totalBitSize = Rows (bestIndividualOverall)*Columns(bestIndividualOverall); + partCount = totalBitSize/bppSize; + outAVL = ExportAMatrix (rawOutFile,bestIndividualOverall,0); + + filterStrings = outAVL["Filters"]; + treeStrings = outAVL["Trees"]; + + results = KHTest(outAVL); + gard.json["likelihoodScores"] = results["likelihoodScores"]; + gard.json["pairwiseP"] = results["pairwiseP"]; + gard.json["singleTreeAICc"] = results["singleTreeAICc"]; + } else { - outAVL = ExportAMatrix (rawOutFile,0,0); + outAVL = ExportAMatrix (rawOutFile,0,0); } +io.SpoolJSON (gard.json, jsonOutFile); + rawOutFile = outputFilePath + "_ga_details"; fprintf (rawOutFile,CLEAR_FILE, KEEP_OPEN); masterKeys = Rows(MasterList); for (h=Rows(masterKeys)*Columns(masterKeys)-1; h>=0; h=h-1) { - aKey = masterKeys[h]; - fprintf (rawOutFile,(-MasterList[aKey]),"\n",aKey,"\n"); + aKey = masterKeys[h]; + fprintf (rawOutFile,(-MasterList[aKey]),"\n",aKey,"\n"); } fprintf (rawOutFile,CLOSE_FILE); IS_TREE_PRESENT_IN_DATA = 0; -DATA_FILE_PRINT_FORMAT = 6; +DATA_FILE_PRINT_FORMAT = 6; -tempFile = outputFilePath + "_finalout"; +tempFile = outputFilePath + "_finalout"; -fprintf (tempFile,CLEAR_FILE,KEEP_OPEN,allData,"\nBEGIN TREES;\n"); -bpL = outAVL ["BP"]; -bpT = outAVL ["Trees"]; -for (rvChoice = 0; rvChoice < Abs(bpT); rvChoice = rvChoice + 1) +fprintf (tempFile,CLEAR_FILE,KEEP_OPEN,allData,"\nBEGIN TREES;\n"); +bpL = outAVL ["BP"]; +bpT = outAVL ["Trees"]; +for (rvChoice = 0; rvChoice < Abs(bpT); rvChoice = rvChoice + 1) { - fprintf (tempFile, "\tTREE part_", rvChoice+1, " = ", bpT[rvChoice], ";\n"); + fprintf (tempFile, "\tTREE part_", rvChoice+1, " = ", bpT[rvChoice], ";\n"); } -fprintf (tempFile,"\nEND;\n\nBEGIN ASSUMPTIONS;\n"); -for (rvChoice = 0; rvChoice < Abs(bpL); rvChoice = rvChoice + 1) +fprintf (tempFile,"\nEND;\n\nBEGIN ASSUMPTIONS;\n"); +for (rvChoice = 0; rvChoice < Abs(bpL); rvChoice = rvChoice + 1) { - fprintf (tempFile, "\tCHARSET span_", rvChoice+1, " = ", bpL[rvChoice], ";\n"); + fprintf (tempFile, "\tCHARSET span_", rvChoice+1, " = ", bpL[rvChoice], ";\n"); } -fprintf (tempFile,"\nEND;\n",CLOSE_FILE); - +fprintf (tempFile,"\nEND;\n",CLOSE_FILE); /*---------------------------------------------------------------------------------------------------------------------------------------------*/ -function reportImprovementHTML (_IC) +function testLRT (vec1, vec2, itCount) { - if (lastImprovedBPC) + size1 = Columns(vec1); + + sumVec1 = {size1,1}; + jvec = {2,size1}; + resMx1 = {itCount,1}; + resMx2 = {itCount,1}; + + for (k=0; kGARD found evidence of "+lastImprovedBPC+" breakpoints

"+spoolAICTable(0) + "\n"; - + resampled = Random(jvec,1); + resampled = resampled*sumVec1; + resMx1[k] = resampled[0]; + resMx2[k] = resampled[1]; } - return "

GARD found no evidence of recombination
"; + + return (resMx1-resMx2)*2; } /*---------------------------------------------------------------------------------------------------------------------------------------------*/ +function KHTest(outAVL) { + + filterStrings = outAVL["Filters"]; + treeStrings = outAVL["Trees"]; + + USE_DISTANCES = 0; + lfDef = ""; + lfDef * 128; + lfDef * "LikelihoodFunction multiPart = ("; + + readPCount = Abs (filterStrings); + + for (pccounter = 0; pccounter < readPCount; pccounter = pccounter + 1) + { + ExecuteCommands ("DataSetFilter part_" + pccounter + " = CreateFilter (ds, 1, \"" + filterStrings[pccounter] + "\");"); + ExecuteCommands ("Tree tree_" + pccounter + " = baseTreeString;"); + if (pccounter) + { + lfDef * ","; + } + lfDef * ("part_"+pccounter+",tree_"+pccounter); + } + lfDef * ");"; + lfDef * 0; + + ExecuteCommands (lfDef); + + Optimize (res2, multiPart); + myDF = res2[1][1]+baseParams; + myAIC = -2*(res2[1][0]-myDF*(baseSites/(baseSites-myDF-1))); + + bpLocations = {readPCount, 1}; + for (pccounter = 0; pccounter < readPCount; pccounter = pccounter + 1) + { + if (pccounter) + { + bpLocations[pccounter] = siteCount + bpLocations[pccounter-1]; + } + ExecuteCommands ("siteCount = part_" + pccounter + ".sites;"); + } + + conditionals = {}; + likelihoodScores = {readPCount,readPCount}; + pairwiseP = {readPCount,readPCount}; + + treeSplitMatches = 0; + khIterations = 10000; + + LIKELIHOOD_FUNCTION_OUTPUT = 7; + + for (pccounter = 0; pccounter < readPCount; pccounter += 1) + { + for (pc2 = 0; pc2 < readPCount; pc2 += 1) + { + if (Abs(pc2-pccounter) <= 1) + { + DataSetFilter aPart = CreateFilter (ds,1,filterStrings[pccounter]); + Tree aTree = treeStrings[pc2]; + LikelihoodFunction aLF = (aPart,aTree); + Optimize (aRes, aLF); + ConstructCategoryMatrix (conds, aLF,SITE_LOG_LIKELIHOODS); + conditionals [pccounter*readPCount+pc2] = conds; + likelihoodScores [pccounter][pc2] = aRes[1][0]; + } + } + + partTreeConds = conditionals[pccounter*readPCount+pccounter]; + + for (pc2 = 0; pc2 < readPCount; pc2 += 1) + { + if (Abs(pc2-pccounter) == 1) + { + otherPartTree = conditionals[pccounter*readPCount+pc2]; + baseLRT = 2*(likelihoodScores[pccounter][pccounter]-likelihoodScores[pccounter][pc2]); + textMx = testLRT(partTreeConds,otherPartTree,khIterations) % 0; + for (kk=0; kk=2*OPTIMIZATION_PRECISION) + { + break; + } + } + pval = Max(1,kk)/khIterations; + pairwiseP[pccounter][pc2] = pval; + } + } + } + + results = { + "likelihoodScores" : likelihoodScores, + "singleTreeAICc" : myAIC, + "pairwiseP": pairwiseP + }; + + return results; -function outputSpan (_offset, _text) -{ - return "" + _text + ""; } /*---------------------------------------------------------------------------------------------------------------------------------------------*/ -function spoolAICTable (dummy) -{ - colorList = {{"red","black","blue","green","purple","orange"}}; - fcolorList = {{"white","white","white","white","white","black"}}; - - htmlAICTable = ""; - htmlAICTable * 128; - +/*---------------------------------------------------------------------------------------------------------------------------------------------*/ - htmlAICTable * "
"; - htmlAICTable * ""; +function reportImprovementHTML (_IC) +{ + if (lastImprovedBPC) + { + return "
GARD found evidence of "+lastImprovedBPC+" breakpoints

"+spoolAICTable(0) + "

\n"; + + } + return "
GARD found no evidence of recombination
"; +} +/*---------------------------------------------------------------------------------------------------------------------------------------------*/ - currentAIC = byBPImprovement [0]; +function outputSpan (_offset, _text) +{ + return "" + _text + ""; +} - for (_partCount = 0; _partCount "; - - if (k"; - prv = bpLocs[k]; - } - sp = sp + "
BPsc-AICΔ c-AICSegments
"; - } - else - { - ci = ""; - sp = "
1-"+allData.sites+"
"; - } - htmlAICTable * ("\n"+ _partCount+ - "
"+currentAIC+ - "
"+ ci+ - ""+ sp+ - ""); - - } +/*---------------------------------------------------------------------------------------------------------------------------------------------*/ - htmlAICTable * "\n
"; - htmlAICTable * 0; - return htmlAICTable; +function spoolAICTable (dummy) +{ + colorList = {{"red","black","blue","green","purple","orange"}}; + fcolorList = {{"white","white","white","white","white","black"}}; + + htmlAICTable = ""; + htmlAICTable * 128; + + + htmlAICTable * "
"; + htmlAICTable * ""; + + + currentAIC = byBPImprovement [0]; + + for (_partCount = 0; _partCount "; + + if (k"; + prv = bpLocs[k]; + } + sp = sp + "
BPsc-AICΔ c-AICSegments
"; + } + else + { + ci = ""; + sp = "
1-"+allData.sites+"
"; + } + htmlAICTable * ("\n"+ _partCount+ + "
"+currentAIC+ + "
"+ ci+ + ""+ sp+ + ""); + + } + + htmlAICTable * "\n
"; + htmlAICTable * 0; + return htmlAICTable; } diff --git a/res/TemplateBatchFiles/LEISR.bf b/res/TemplateBatchFiles/LEISR.bf index e6a3df4e4..f6712148c 100644 --- a/res/TemplateBatchFiles/LEISR.bf +++ b/res/TemplateBatchFiles/LEISR.bf @@ -20,6 +20,8 @@ LoadFunctionLibrary("libv3/models/DNA/JC69.bf"); LoadFunctionLibrary("libv3/models/protein.bf"); LoadFunctionLibrary("libv3/models/protein/empirical.bf"); +// for JSON storage compatibility +LoadFunctionLibrary("SelectionAnalyses/modules/io_functions.ibf"); /*------------------------------------------------------------------------------*/ @@ -28,13 +30,38 @@ utility.ToggleEnvVariable ("NORMALIZE_SEQUENCE_NAMES", 1); leisr.analysis_description = { terms.io.info: "LEISR (Likelihood Estimation of Individual Site Rates) infer relative amino-acid or nucleotide rates from a fixed nucleotide or amino-acid alignment and tree. Relative site-specific substitution rates are inferred by first optimizing alignment-wide branch lengths, and then inferring a site-specific uniform tree scaler", - terms.io.version: "0.1alpha", - terms.io.reference: "@TBD. Analysis based on Rate4Site method, : Pupko, T., Bell, R. E., Mayrose, I., Glaser, F. & Ben-Tal, N. Rate4Site: an algorithmic tool for the identification of functional regions in proteins by surface mapping of evolutionary determinants within their homologues. Bioinformatics 18, S71–S77 (2002).", + terms.io.version: "0.2", + terms.io.reference: "Spielman, S.J. and Kosakovsky Pond, S.L. Relative evolutionary rate inference in HyPhy with LEISR. bioRxiv. https://doi.org/10.1101/206011. (2017); Pupko, T., Bell, R. E., Mayrose, I., Glaser, F. & Ben-Tal, N. Rate4Site: an algorithmic tool for the identification of functional regions in proteins by surface mapping of evolutionary determinants within their homologues. Bioinformatics 18, S71–S77 (2002).", terms.io.authors: "Sergei L Kosakovsky Pond and Stephanie J Spielman", terms.io.contact: "{spond,stephanie.spielman}@temple.edu" }; io.DisplayAnalysisBanner(leisr.analysis_description); +/*******************************************************************************************************************/ + + +/***************************************** LOAD DATASET **********************************************************/ +SetDialogPrompt ("Specify a multiple sequence alignment file"); +leisr.alignment_info = alignments.ReadNucleotideDataSet ("leisr.dataset", None); + +leisr.name_mapping = leisr.alignment_info[utility.getGlobalValue("terms.data.name_mapping")]; +if (None == leisr.name_mapping) { + leisr.name_mapping = {}; + utility.ForEach (alignments.GetSequenceNames ("leisr.dataset"), "_value_", "`&leisr.name_mapping`[_value_] = _value_"); +} +leisr.partitions_and_trees = trees.LoadAnnotatedTreeTopology.match_partitions (leisr.alignment_info[utility.getGlobalValue("terms.data.partitions")], leisr.name_mapping); +leisr.partition_count = Abs (leisr.partitions_and_trees); + +io.CheckAssertion ("leisr.partition_count==1", "This analysis can only handle a single partition"); + + +io.ReportProgressMessageMD ("relative_rates", "Data", "Input alignment description"); +io.ReportProgressMessageMD ("relative_rates", "Data", "Loaded **" + + leisr.alignment_info [terms.data.sequences] + "** sequences, **" + + leisr.alignment_info [terms.data.sites] + "** sites, and **" + leisr.partition_count + "** partitions from \`" + leisr.alignment_info [terms.data.file] + "\`"); +leisr.filter_specification = alignments.DefineFiltersForPartitions (leisr.partitions_and_trees, "leisr.dataset" , "leisr.filter.", leisr.alignment_info); +/*******************************************************************************************************************/ + /***************************************** MODEL SELECTION **********************************************************/ @@ -45,40 +72,20 @@ leisr.analysis_type = io.SelectAnOption ({{leisr.protein_type , "Infer relative if (leisr.analysis_type == leisr.protein_type) { - leisr.baseline_model = io.SelectAnOption (models.protein.empirical_models, - "Select a protein model:"); - - leisr.use_rate_variation = io.SelectAnOption( {{"Gamma", "Use a four-category discrete gamma distribution when optimizing branch lengths."}, - {"GDD", "Use a four-category general discrete distribution when optimizing branch lengths."}, - {"No", "Do not consider rate variation when optimizing branch lengths."} - }, - "Optimize branch lengths with rate variation?"); - // "Yes", "No" - leisr.plusF = io.SelectAnOption ({{"Yes", "Use empirical (+F) amino-acid frequencies ."}, {"No", "Use default amino-acid frequencies."}}, - "Use a +F model for initial branch length optimization?"); - if (leisr.plusF == "Yes"){ - leisr.generators = models.protein.empirical.plusF_generators; - } - else { - leisr.generators = models.protein.empirical.default_generators; - } -} + leisr.baseline_model = io.SelectAnOption (models.protein.empirical_models, "Select a protein model:"); + leisr.generators = models.protein.empirical.default_generators; +} else { - leisr.baseline_model = io.SelectAnOption (models.DNA.models, - "Select a nucleotide model:"); - leisr.use_rate_variation = io.SelectAnOption( {{"Gamma", "Use a four-category discrete gamma distribution when optimizing branch lengths."}, - {"GDD", "Use a four-category general discrete distribution when optimizing branch lengths."}, - {"No", "Do not consider rate variation when optimizing branch lengths."} - }, - "Optimize branch lengths with rate variation?"); - + leisr.baseline_model = io.SelectAnOption (models.DNA.models, "Select a nucleotide model:"); leisr.generators = models.DNA.generators; - leisr.plusF = "No"; } - - +leisr.use_rate_variation = io.SelectAnOption( {{"Gamma", "Use a four-category discrete gamma distribution when optimizing branch lengths."}, + {"GDD", "Use a four-category general discrete distribution when optimizing branch lengths."}, + {"No", "Do not consider rate variation when optimizing branch lengths."} + }, + "Optimize branch lengths with rate variation?"); function leisr.Baseline.ModelDescription(type){ @@ -99,52 +106,23 @@ function leisr.Baseline.ModelDescription.withGDD4(type){ leisr.baseline_model_name = leisr.baseline_model; -if (leisr.plusF == "Yes"){ - leisr.baseline_model_name = leisr.baseline_model_name + "+F"; +if (leisr.analysis_type == leisr.protein_type) { + leisr.baseline_model_name = leisr.baseline_model_name + "F"; } - if (leisr.use_rate_variation == "Gamma"){ - leisr.baseline_model_name = leisr.baseline_model_name + " with 4 category Gamma rates"; + leisr.baseline_model_name = leisr.baseline_model_name + "+4Gamma"; leisr.baseline_model_desc = "leisr.Baseline.ModelDescription.withGamma"; -} +} else { if (leisr.use_rate_variation == "GDD"){ - leisr.baseline_model_name = leisr.baseline_model_name + " with 4 category GDD rates"; + leisr.baseline_model_name = leisr.baseline_model_name + "+4GDD"; leisr.baseline_model_desc = "leisr.Baseline.ModelDescription.withGDD4"; - } + } else { leisr.baseline_model_name = leisr.baseline_model_name; leisr.baseline_model_desc = "leisr.Baseline.ModelDescription"; } } -/**************************************************************/ - - - - -/*******************************************************************************************************************/ - - -/***************************************** LOAD DATASET **********************************************************/ -SetDialogPrompt ("Specify a multiple sequence alignment file"); -leisr.alignment_info = alignments.ReadNucleotideDataSet ("leisr.dataset", NOne); - -name_mapping = leisr.alignment_info[utility.getGlobalValue("terms.data.name_mapping")]; -if (None == name_mapping) { - name_mapping = {}; - utility.ForEach (alignments.GetSequenceNames ("leisr.dataset"), "_value_", "`&name_mapping`[_value_] = _value_"); -} -leisr.partitions_and_trees = trees.LoadAnnotatedTreeTopology.match_partitions (leisr.alignment_info[utility.getGlobalValue("terms.data.partitions")], name_mapping); -leisr.partition_count = Abs (leisr.partitions_and_trees); - -io.CheckAssertion ("leisr.partition_count==1", "This analysis can only handle a single partition"); - - -io.ReportProgressMessageMD ("relative_rates", "Data", "Input alignment description"); -io.ReportProgressMessageMD ("relative_rates", "Data", "Loaded **" + - leisr.alignment_info [terms.data.sequences] + "** sequences, **" + - leisr.alignment_info [terms.data.sites] + "** sites, and **" + leisr.partition_count + "** partitions from \`" + leisr.alignment_info [terms.data.file] + "\`"); -leisr.filter_specification = alignments.DefineFiltersForPartitions (leisr.partitions_and_trees, "leisr.dataset" , "leisr.filter.", leisr.alignment_info); /*******************************************************************************************************************/ @@ -158,24 +136,28 @@ leisr.trees = utility.Map (leisr.partitions_and_trees, "_value_", "_value_[terms leisr.filter_names = utility.Map (leisr.filter_specification, "_value_", "_value_[terms.data.name]"); // value => value['name'] leisr.alignment_wide_MLES = estimators.FitSingleModel_Ext ( leisr.filter_names, - leisr.trees, + leisr.trees, leisr.baseline_model_desc, None, None); - - - + + + estimators.fixSubsetOfEstimates(leisr.alignment_wide_MLES, leisr.alignment_wide_MLES[terms.global]); io.ReportProgressMessageMD ("relative_rates", "overall", ">Fitted an alignment-wide model. **Log-L = " + leisr.alignment_wide_MLES [terms.fit.log_likelihood] + "**."); -/** - Set up the table to display to the screen +/** + Set up the table to display to the screen */ leisr.table_screen_output = {{"Site", "Rel. rate (MLE)", "95% profile likelihood CI"}}; leisr.table_output_options = {terms.table_options.header : TRUE, terms.table_options.minimum_column_width : 16, terms.table_options.align : "center"}; +leisr.table_headers = {{"MLE", "Relative rate estimate at a site"} + {"Lower", "Lower bound of 95% profile likelihood CI"} + {"Upper", "Upper bound of 95% profile likelihood CI"}}; +leisr.site_results = {leisr.alignment_info [terms.data.sites], Rows (leisr.table_headers)}; leisr.site_patterns = alignments.Extract_site_patterns (leisr.filter_names[0]); @@ -186,17 +168,17 @@ leisr.site_model = model.generic.DefineModel("leisr.Baseline.ModelDescription", }, leisr.filter_names[0], None); - - + + leisr.site_model_mapping = {"relative_rates_site_model_instance" : leisr.site_model}; - + // leisr.site_tree is created from the information in leisr.trees[0] -// and populated with (the default) model +// and populated with (the default) model model.ApplyModelToTree( "leisr.site_tree", leisr.trees[0], {terms.default : leisr.site_model}, None); // create a site filter; this is an ugly hack for the time being -// alignments.serialize_site_filter returns HBL code as string in +// alignments.serialize_site_filter returns HBL code as string in // which the function `__make_filter` is defined. ExecuteCommands (alignments.serialize_site_filter ( leisr.filter_names[0], @@ -208,7 +190,7 @@ LikelihoodFunction leisr.site_likelihood = (leisr.site_filter, leisr.site_tree); leisr.site_model_scaler_name = "leisr.site_rate_estimate"; -leisr.rate_estimates = {}; +leisr.rate_estimates = {}; /** this will store site estimates, which will then be dumped to JSON @@ -219,7 +201,7 @@ parameters.DeclareGlobal (leisr.site_model_scaler_name, None); estimators.ApplyExistingEstimates ("leisr.site_likelihood", leisr.site_model_mapping, leisr.alignment_wide_MLES, {"0" : leisr.site_model_scaler_name} // proportional scaler ); - + leisr.queue = mpi.CreateQueue ({terms.mpi.LikelihoodFunctions: {{"leisr.site_likelihood"}}, terms.mpi.Models : {{"leisr.site_model"}}, @@ -252,30 +234,51 @@ io.ReportProgressMessageMD ("relative_rates", "Stats", "* **Median**: " + Forma io.ReportProgressMessageMD ("relative_rates", "Stats", "* **Std.Dev**: " + Format (leisr.stats[terms.math.stddev], 6, 2)); io.ReportProgressMessageMD ("relative_rates", "Stats", "* **95% Range**: [" + Format (leisr.stats[terms.math._2.5], 5,2) + "," + Format (leisr.stats[terms.math._97.5], 5,2) + "]"); - -if (leisr.use_rate_variation == "No"){ - leisr.storerv = "None"; -} else { - leisr.storerv = utility.Map( - utility.Map (leisr.alignment_wide_MLES[utility.getGlobalValue("terms.global")], "_value_", ' {terms.fit.MLE : _value_[terms.fit.MLE]}'), - "_value_", - "_value_[terms.fit.MLE]"); -} -tree_definition = utility.Map (leisr.partitions_and_trees, "_partition_", '_partition_[terms.data.tree]'); -io.SpoolJSON ({ terms.json.input : {terms.json.file: leisr.alignment_info[terms.data.file], - terms.json.sequences: leisr.alignment_info[terms.data.sequences], - terms.json.sites: leisr.alignment_info[terms.data.sites], - terms.json.tree_string: (tree_definition[0])[terms.trees.newick_with_lengths]}, - terms.json.analysis : leisr.analysis_description, - terms.json.relative_site_rates : leisr.rate_estimates, - terms.json.global: {terms.json.model: leisr.baseline_model_name, - terms.model.rate_variation: leisr.storerv, - terms.efv_estimate: (leisr.alignment_wide_MLES[utility.getGlobalValue("terms.efv_estimate")])["VALUEINDEXORDER"][0], - terms.json.tree_string: (leisr.alignment_wide_MLES[terms.fit.trees])[0], - terms.json.log_likelihood: leisr.alignment_wide_MLES[terms.fit.log_likelihood]} - }, - leisr.alignment_info[terms.data.file] + ".site-rates.json"); + +/************************************* JSON STORAGE ***********************************/ + +leisr.store_global = utility.Map( + utility.Map (leisr.alignment_wide_MLES[utility.getGlobalValue("terms.global")], "_value_", ' {terms.fit.MLE : _value_[terms.fit.MLE]}'), + "_value_", + "_value_[terms.fit.MLE]"); + +if (Abs(((leisr.partitions_and_trees["0"])[terms.data.tree])[terms.branch_length]) == 0){ + store_tree = utility.Map (leisr.partitions_and_trees, "_pt_", '(_pt_[terms.data.tree])[terms.trees.newick]'); +} else { + store_tree = utility.Map (leisr.partitions_and_trees, "_pt_", '(_pt_[terms.data.tree])[terms.trees.newick_with_lengths]'); +} + +leisr.aicc = leisr.getIC(leisr.alignment_wide_MLES[terms.fit.log_likelihood], leisr.alignment_wide_MLES[terms.parameters], leisr.alignment_info[utility.getGlobalValue("terms.data.sites")] * leisr.alignment_info[utility.getGlobalValue("terms.data.sequences")]); +leisr.json_content = { terms.json.input : + {terms.json.file: leisr.alignment_info[terms.data.file], + terms.json.sequences: leisr.alignment_info[terms.data.sequences], + terms.json.sites: leisr.alignment_info[terms.data.sites], + terms.json.partition_count: leisr.partition_count, + terms.json.trees: store_tree + }, + terms.json.analysis : leisr.analysis_description, + terms.json.fits: + { + leisr.baseline_model_name: + { + terms.json.log_likelihood: leisr.alignment_wide_MLES[terms.fit.log_likelihood], + terms.json.parameters: leisr.alignment_wide_MLES[terms.parameters], + terms.json.AICc: leisr.aicc, + terms.json.rate_distribution: leisr.store_global, + terms.efv_estimate: (leisr.alignment_wide_MLES[utility.getGlobalValue("terms.efv_estimate")])["VALUEINDEXORDER"][0] + } + }, + terms.json.MLE : {terms.json.headers : leisr.table_headers, + terms.json.content : {"0":leisr.site_results} + } + }; +selection.io.json_store_branch_attribute(leisr.json_content, utility.getGlobalValue ("terms.original_name"), utility.getGlobalValue ("terms.json.node_label"), -1, 0, utility.getGlobalValue ("leisr.name_mapping")); +selection.io.json_store_branch_attribute(leisr.json_content, utility.getGlobalValue ("leisr.baseline_model_name"), utility.getGlobalValue ("terms.branch_length"), 0, 0, + selection.io.extract_branch_info((leisr.alignment_wide_MLES[terms.branch_length])[0], "selection.io.branch.length") + ); + +io.SpoolJSON (leisr.json_content, leisr.alignment_info[terms.data.file] + ".LEISR.json"); //---------------------------------------------------------------------------------------- @@ -288,25 +291,23 @@ io.SpoolJSON ({ terms.json.input : {terms.json.file: leisr.alignment_info[terms. //---------------------------------------------------------------------------------------- lfunction leisr.handle_a_site (lf, filter_data, pattern_info, model_mapping) { - GetString (lfInfo, ^lf,-1); ExecuteCommands (filter_data); - + __make_filter ((lfInfo["Datafilters"])[0]); utility.SetEnvVariable ("USE_LAST_RESULTS", TRUE); - + if (pattern_info [utility.getGlobalValue("terms.data.is_constant")]) { - // the MLE for a constant site is 0; + // the MLE for a constant site is 0; // only the CI is non-trivial ^(utility.getGlobalValue("leisr.site_model_scaler_name")) = 0; - + } else { - + ^(utility.getGlobalValue("leisr.site_model_scaler_name")) = 1; Optimize (results, ^lf); } - return parameters.GetProfileCI (utility.getGlobalValue("leisr.site_model_scaler_name"), lf, 0.95); } @@ -322,15 +323,14 @@ lfunction leisr.store_results (node, result, arguments) { if ((^'leisr.table_output_options')[utility.getGlobalValue("terms.table_options.header")]) { - - io.ReportProgressMessageMD ("relative_rates", "sites", "Site rate estimates and associated 95% profile likelihood estimates\n"); + + io.ReportProgressMessageMD ("relative_rates", "sites", "Site rate estimates and associated 95% profile likelihood estimates\n"); fprintf (stdout, io.FormatTableRow (^'leisr.table_screen_output',^'leisr.table_output_options')); (^'leisr.table_output_options')[utility.getGlobalValue("terms.table_options.header")] = FALSE; } - - + utility.ForEach (pattern_info[utility.getGlobalValue("terms.data.sites")], "_site_index_", " leisr.rate_estimates [_site_index_+1] = `&result`; @@ -340,6 +340,12 @@ lfunction leisr.store_results (node, result, arguments) { result_row [2] = Format((`&result`)[terms.lower_bound],6,3) + ' :' +Format((`&result`)[terms.upper_bound],6,3); fprintf (stdout, io.FormatTableRow (result_row,leisr.table_output_options)); + + // JSON-related + //console.log(_site_index_); + leisr.site_results[_site_index_][0] = (leisr.rate_estimates[_site_index_+1])[terms.fit.MLE]; + leisr.site_results[_site_index_][1] = (leisr.rate_estimates[_site_index_+1])[terms.lower_bound]; + leisr.site_results[_site_index_][2] = (leisr.rate_estimates[_site_index_+1])[terms.upper_bound]; " ); @@ -349,3 +355,6 @@ lfunction leisr.store_results (node, result, arguments) { +lfunction leisr.getIC(logl, params, samples) { + return -2 * logl + 2 * samples / (samples - params - 1) * params; +} \ No newline at end of file diff --git a/res/TemplateBatchFiles/SelectionAnalyses/FUBAR.bf b/res/TemplateBatchFiles/SelectionAnalyses/FUBAR.bf index 93f2a4b48..5daca9303 100644 --- a/res/TemplateBatchFiles/SelectionAnalyses/FUBAR.bf +++ b/res/TemplateBatchFiles/SelectionAnalyses/FUBAR.bf @@ -39,6 +39,7 @@ namespace terms.fubar { }; fubar.json = { + terms.json.analysis: {}, terms.json.input: {}, terms.json.fits : {}, terms.json.timers : {}, @@ -97,6 +98,8 @@ fubar.analysis_description = {terms.io.info : io.DisplayAnalysisBanner ( fubar.analysis_description ); +fubar.json[terms.json.analysis] = fubar.analysis_description; + namespace fubar { LoadFunctionLibrary ("modules/shared-load-file.bf"); load_file ({utility.getGlobalValue("terms.prefix"): "fubar", utility.getGlobalValue("terms.settings") : {utility.getGlobalValue("terms.settings.branch_selector") : "selection.io.SelectAllBranches"}}); diff --git a/res/TemplateBatchFiles/SelectionAnalyses/SLAC.bf b/res/TemplateBatchFiles/SelectionAnalyses/SLAC.bf index 59382775c..dbc4b6d46 100644 --- a/res/TemplateBatchFiles/SelectionAnalyses/SLAC.bf +++ b/res/TemplateBatchFiles/SelectionAnalyses/SLAC.bf @@ -75,12 +75,12 @@ slac.json = { /** The dictionary of results to be written to JSON at the end of the run */ - + slac.display_orders = {terms.original_name: -1, terms.json.nucleotide_gtr: 0, terms.json.global_mg94xrev: 1 }; - + selection.io.startTimer (slac.json [terms.json.timers], "Total time", 0); @@ -264,13 +264,13 @@ for (slac.i = 0; slac.i < Abs (slac.filter_specification); slac.i += 1) { } slac.branch_attributes = selection.substitution_mapper (slac.ancestors["MATRIX"], slac.ancestors["TREE_AVL"], slac.ancestors["AMBIGS"], slac.counts, slac.ancestors ["MAPPING"], slac.codon_data_info[terms.code]); - + /* selection.io.json_store_branch_attribute(slac.json, terms.original_name, terms.json.node_label, 0, slac.i, slac.name_mapping); */ - + selection.io.json_store_branch_attribute(slac.json, terms.codon, terms.json.node_label, 0, slac.i, slac.branch_attributes[terms.codon]); @@ -329,7 +329,11 @@ if (slac.samples > 0) { slac.table_output_options_samplers[terms.table_options.header] = TRUE; slac.sample.results = {}; - slac.queue = mpi.CreateQueue ({"LikelihoodFunctions": {{slac.partitioned_mg_results[terms.likelihood_function]}}}); + slac.queue = mpi.CreateQueue ({"LikelihoodFunctions": {{slac.partitioned_mg_results[terms.likelihood_function]}}, + "Headers" : {{"libv3/all-terms.bf"}}, + "Variables" : {{"slac.by_site","slac.AVERAGED","slac.RESOLVED","slac.by_branch"}} + }); + for (slac.s = 0; slac.s < slac.samples; slac.s+=1) { diff --git a/res/TemplateBatchFiles/libv3/all-terms.bf b/res/TemplateBatchFiles/libv3/all-terms.bf index 5586e4420..20551d23b 100644 --- a/res/TemplateBatchFiles/libv3/all-terms.bf +++ b/res/TemplateBatchFiles/libv3/all-terms.bf @@ -84,6 +84,9 @@ namespace terms{ /* Term functions */ + function characterFrequency (character) { + return "Equilibrium frequency for " + character; + } function nucleotideRate(fromC, toC) { return "Substitution rate from nucleotide " + fromC + " to nucleotide " + toC; } @@ -172,6 +175,7 @@ namespace terms{ equal = "Equal frequencies"; CF3x4 = "Corrected 3x4 frequency estimator"; _20x1 = "Protein 20x1 estimator"; + MLE = "Maximum likelihood frequency estimator"; predefined = "Based on a training set"; binary = "Binary character frequency estimator"; } diff --git a/res/TemplateBatchFiles/libv3/models/frequencies.bf b/res/TemplateBatchFiles/libv3/models/frequencies.bf index 458b33012..8b9174eb4 100644 --- a/res/TemplateBatchFiles/libv3/models/frequencies.bf +++ b/res/TemplateBatchFiles/libv3/models/frequencies.bf @@ -70,6 +70,44 @@ function frequencies.empirical.protein (model, namespace, datafilter) { return model; } +/** + * Sets model's equilibrium frequency estimator to ML for protein data + * @name frequencies.empirical.protein + * @param {Dictionary} model + * @param {String} namespace + * @param {DataSetFilter} datafilter + * @returns {Dictionary} updated model + */ +function frequencies.ML.protein (model, namespace, datafilter) { + model = frequencies._aux.empirical.singlechar(model, namespace, datafilter); + // define 20 frequency parameters + // initialize to empirical freqs + // add to model parameter manifest + frequencies.ML.protein.emp = model[terms.efv_estimate]; + model[terms.efv_estimate] = {20,1}; + frequencies.ML.protein.variables = {20,1}; + frequencies.ML.protein.scaler = namespace + ".frequency_scaler"; + + utility.ForEachPair (model[terms.alphabet], "_index", "_letter", + ' + _idx = _index[1]; + _freq_parameter = namespace + ".equilibrium_frequency_of." + _letter; + frequencies.ML.protein.variables [_idx] = _freq_parameter; + (model[terms.efv_estimate]) [_idx] = _freq_parameter + "/" + frequencies.ML.protein.scaler; + parameters.DeclareGlobalWithRanges (_freq_parameter, frequencies.ML.protein.emp[_idx], 0, 1); + model.generic.AddGlobal (model, _freq_parameter, terms.characterFrequency (_letter)); + ' + ); + + // constrain pi_A + ... + pi_A = 1, + + parameters.SetConstraint ( frequencies.ML.protein.scaler, Join ("+", frequencies.ML.protein.variables), "global"); + + model[terms.model.efv_estimate_name] = terms.frequencies.ml; + (model[terms.parameters])[terms.model.empirical] = -1; // correct for the restricted sum + return model; +} + /** * @name frequencies.empirical.corrected.CF3x4 * @param {Dictionary} model diff --git a/res/TemplateBatchFiles/libv3/models/model_functions.bf b/res/TemplateBatchFiles/libv3/models/model_functions.bf index 8ef51e878..7a3e2a872 100644 --- a/res/TemplateBatchFiles/libv3/models/model_functions.bf +++ b/res/TemplateBatchFiles/libv3/models/model_functions.bf @@ -228,7 +228,7 @@ function model.generic.DefineModel (model_spec, id, arguments, data_filter, esti // Set EFV field - Call (model.generic.DefineModel.model [terms.model.frequency_estimator], model.generic.DefineModel.model, + model.generic.DefineModel.model = Call (model.generic.DefineModel.model [terms.model.frequency_estimator], model.generic.DefineModel.model, id, data_filter); @@ -240,10 +240,14 @@ function model.generic.DefineModel (model_spec, id, arguments, data_filter, esti model.generic.DefineModel.model [terms.id] = id; - parameters.StringMatrixToFormulas (model.generic.DefineModel.model [terms.model.matrix_id],model.generic.DefineModel.model[terms.model.rate_matrix]); - utility.SetEnvVariable (model.generic.DefineModel.model [terms.model.efv_id], model.generic.DefineModel.model[terms.efv_estimate]); - + + if (Type ((model.generic.DefineModel.model[terms.efv_estimate])[0]) == "String") { + parameters.StringMatrixToFormulas (model.generic.DefineModel.model [terms.model.efv_id],model.generic.DefineModel.model[terms.efv_estimate]); + } else { + utility.SetEnvVariable (model.generic.DefineModel.model [terms.model.efv_id], model.generic.DefineModel.model[terms.efv_estimate]); + } + model.define_from_components (id, model.generic.DefineModel.model [terms.model.matrix_id], model.generic.DefineModel.model [terms.model.efv_id], model.generic.DefineModel.model [terms.model.canonical]); diff --git a/res/TemplateBatchFiles/libv3/models/parameters.bf b/res/TemplateBatchFiles/libv3/models/parameters.bf index f835d635b..382c01f13 100644 --- a/res/TemplateBatchFiles/libv3/models/parameters.bf +++ b/res/TemplateBatchFiles/libv3/models/parameters.bf @@ -209,7 +209,16 @@ function parameters.SetValues(set) { * @returns nothing */ lfunction parameters.ConstrainMeanOfSet (set, mean, namespace) { - unscaled = utility.Map (utility.Values (set), "_name_", "_name_ + '_scaler_variable'"); + if (Type (set) == "AssociativeList") { + unscaled = utility.Map (utility.Values (set), "_name_", "_name_ + '_scaler_variable'"); + } else { + if (Type (set) == "Matrix") { + unscaled = utility.Map (set, "_name_", "_name_ + '_scaler_variable'"); + } + else { + return; + } + } global_scaler = namespace + ".scaler_variable"; parameters.SetConstraint (global_scaler, Join ("+", unscaled), "global"); utility.ForEach (set, "_name_", ' @@ -287,18 +296,18 @@ lfunction parameters.AddMultiplicativeTerm(matrix, term, do_empties) { */ function parameters.StringMatrixToFormulas(id, matrix) { __N = Rows(matrix); + __M = Columns(matrix); - ExecuteCommands("`id` = {__N,__N}"); + ExecuteCommands("`id` = {__N,__M}"); for (__r = 0; __r < __N; __r += 1) { - for (__c = 0; __c < __N; __c += 1) { + for (__c = 0; __c < __M; __c += 1) { - if (__r != __c && Abs(matrix[__r][__c])) { + if (Abs(matrix[__r][__c])) { ExecuteCommands("`id`[__r][__c] := " + matrix[__r][__c]); } } } - } /** @@ -399,6 +408,7 @@ lfunction parameters.GetConstraint(parameter) { function parameters.SetConstraint(id, value, global_tag) { if (Type(id) == "String") { if (Abs(id)) { + //console.log ("`global_tag` `id` := " + value); ExecuteCommands("`global_tag` `id` := " + value); } } else { diff --git a/res/TemplateBatchFiles/libv3/models/protein.bf b/res/TemplateBatchFiles/libv3/models/protein.bf index dd00da4cd..170c2be49 100644 --- a/res/TemplateBatchFiles/libv3/models/protein.bf +++ b/res/TemplateBatchFiles/libv3/models/protein.bf @@ -14,7 +14,12 @@ models.protein.empirical_models = {{"LG", "Empirical model of protein evolution {"WAG", "Empirical model of protein evolution from Whelan and Goldman (2001). Ref: https://doi.org/10.1093/oxfordjournals.molbev.a003851"}, {"JTT", "Empirical model of protein evolution from Jones, Taylor, and Thornton (1996). Ref: https://doi.org/10.1093/bioinformatics/8.3.275"}, {"JC69", "Empirical model of protein evolution with equal exchangeability rates among all amino acids, also known as JC69."}, - {"mtMAM", "Empirical model of protein evolution for mammalian mitochondrial genomes from Yang N, Nielsen R, and Hasegawa M. (1998). Ref: http://dx.doi.org/10.1093/oxfordjournals.molbev.a025888"}}; + {"mtMAM", "Empirical model of protein evolution for mammalian mitochondrial genomes from Yang N, Nielsen R, and Hasegawa M. (1998). Ref: http://dx.doi.org/10.1093/oxfordjournals.molbev.a025888"}, + {"cpREV", "Empirical model of protein evolution for chloroplast genomes from from Adachi et al. (2000). Ref: https://www.ncbi.nlm.nih.gov/pubmed/10795826"}, + {"HIVBm", "Empirical model of protein evolution for between-host HIV sequences from Nickle et al. (2007). Ref: https://doi.org/10.1371/journal.pone.0000503"}, + {"HIVWm", "Empirical model of protein evolution for within-host HIV sequences from Nickle et al. (2007). Ref: https://doi.org/10.1371/journal.pone.0000503"}, + {"AB", "Empirical model of protein evolution for antibody sequences from Mirsky et al. (2015). Ref: https://doi.org/10.1093/molbev/msu340"} + }; models.protein.dimensions = 20; diff --git a/res/TemplateBatchFiles/libv3/models/protein/empirical.bf b/res/TemplateBatchFiles/libv3/models/protein/empirical.bf index 3957c79f9..43777c43d 100644 --- a/res/TemplateBatchFiles/libv3/models/protein/empirical.bf +++ b/res/TemplateBatchFiles/libv3/models/protein/empirical.bf @@ -10,14 +10,31 @@ models.protein.empirical.default_generators = {"LG": "models.protein.LG.ModelDes "WAG": "models.protein.WAG.ModelDescription", "JTT": "models.protein.JTT.ModelDescription", "JC69": "models.protein.JC69.ModelDescription", - "mtMAM": "models.protein.mtMAM.ModelDescription"}; + "mtMAM": "models.protein.mtMAM.ModelDescription", + "cpREV": "models.protein.cpREV.ModelDescription", + "HIVBm": "models.protein.HIVBm.ModelDescription", + "HIVWm": "models.protein.HIVWm.ModelDescription", + "AB" : "models.protein.AB.ModelDescription"}; models.protein.empirical.plusF_generators = {"LG": "models.protein.LGF.ModelDescription", "WAG": "models.protein.WAGF.ModelDescription", "JTT": "models.protein.JTTF.ModelDescription", "JC69": "models.protein.JC69F.ModelDescription", - "mtMAM": "models.protein.mtMAMF.ModelDescription"}; - + "mtMAM": "models.protein.mtMAMF.ModelDescription", + "cpREV": "models.protein.cpREVF.ModelDescription", + "HIVBm": "models.protein.HIVBmF.ModelDescription", + "HIVWm": "models.protein.HIVWmF.ModelDescription", + "AB" : "models.protein.ABF.ModelDescription"}; + +models.protein.empirical.mleF_generators = {"LG": "models.protein.LGML.ModelDescription", + "WAG": "models.protein.WAGML.ModelDescription", + "JTT": "models.protein.JTTML.ModelDescription", + "JC69": "models.protein.JC69ML.ModelDescription", + "mtMAM": "models.protein.mtMAMML.ModelDescription", + "cpREV": "models.protein.cpREVML.ModelDescription", + "HIVBm": "models.protein.HIVBmML.ModelDescription", + "HIVWm": "models.protein.HIVWmML.ModelDescription", + "AB" : "models.protein.ABML.ModelDescription"}; /** @module models.protein.empirical */ /** @@ -295,6 +312,19 @@ function models.protein.WAGF.ModelDescription(type) { return models.protein.WAGF.ModelDescription.model_definition; } +/** + * @name models.protein.WAGML.ModelDescription + * @description Create the baseline schema (dictionary) for the WAG+ML model of protein evolution + * @returns {Dictionary} model description + * @param {String} type + */ +function models.protein.WAGML.ModelDescription(type) { + models.protein.WAGML.ModelDescription.model_definition = models.protein.WAG.ModelDescription(type); + models.protein.WAGML.ModelDescription.model_definition [terms.model.frequency_estimator] = "frequencies.ML.protein"; + models.protein.WAGML.ModelDescription.model_definition [terms.model.efv_estimate_name] = utility.getGlobalValue("terms.frequencies.MLE"); + return models.protein.WAGML.ModelDescription.model_definition; +} + @@ -328,6 +358,19 @@ function models.protein.LGF.ModelDescription(type) { } +/** + * @name models.protein.LGML.ModelDescription + * @description Create the baseline schema (dictionary) for the LG+ML model of protein evolution + * @returns {Dictionary} model description + * @param {String} type + */ +function models.protein.LGML.ModelDescription(type) { + models.protein.LGML.ModelDescription.model_definition = models.protein.LG.ModelDescription(type); + models.protein.LGML.ModelDescription.model_definition [terms.model.frequency_estimator] = "frequencies.ML.protein"; + models.protein.LGML.ModelDescription.model_definition [terms.model.efv_estimate_name] = utility.getGlobalValue("terms.frequencies.MLE"); + return models.protein.LGML.ModelDescription.model_definition; +} + /**************************************** JTT functions *************************************/ @@ -360,6 +403,19 @@ function models.protein.JTTF.ModelDescription(type) { } +/** + * @name models.protein.JTTML.ModelDescription + * @description Create the baseline schema (dictionary) for the JTT+ML model of protein evolution + * @returns {Dictionary} model description + * @param {String} type + */ +function models.protein.JTTML.ModelDescription(type) { + models.protein.JTTML.ModelDescription.model_definition = models.protein.JTT.ModelDescription(type); + models.protein.JTTML.ModelDescription.model_definition [terms.model.frequency_estimator] = "frequencies.ML.protein"; + models.protein.JTTML.ModelDescription.model_definition [terms.model.efv_estimate_name] = utility.getGlobalValue("terms.frequencies.MLE"); + return models.protein.JTTML.ModelDescription.model_definition; +} + /**************************************** JC69 functions *************************************/ @@ -390,6 +446,20 @@ function models.protein.JC69F.ModelDescription(type) { } +/** + * @name models.protein.JC69ML.ModelDescription + * @description Create the baseline schema (dictionary) for the JC69+ML model of protein evolution + * @returns {Dictionary} model description + * @param {String} type + */ +function models.protein.JC69ML.ModelDescription(type) { + models.protein.JC69ML.ModelDescription.model_definition = models.protein.JC69.ModelDescription(type); + models.protein.JC69ML.ModelDescription.model_definition [terms.model.frequency_estimator] = "frequencies.ML.protein"; + models.protein.JC69ML.ModelDescription.model_definition [terms.model.efv_estimate_name] = utility.getGlobalValue("terms.frequencies.MLE"); + return models.protein.JC69ML.ModelDescription.model_definition; +} + + /**************************************** mtMAM functions *************************************/ @@ -421,6 +491,195 @@ function models.protein.mtMAMF.ModelDescription(type) { } +/** + * @name models.protein.mtMAMML.ModelDescription + * @description Create the baseline schema (dictionary) for the mtMAM+ML model of protein evolution + * @returns {Dictionary} model description + * @param {String} type + */ +function models.protein.mtMAMML.ModelDescription(type) { + models.protein.mtMAMML.ModelDescription.model_definition = models.protein.mtMAM.ModelDescription(type); + models.protein.mtMAMML.ModelDescription.model_definition [terms.model.frequency_estimator] = "frequencies.ML.protein"; + models.protein.mtMAMML.ModelDescription.model_definition [terms.model.efv_estimate_name] = utility.getGlobalValue("terms.frequencies.MLE"); + return models.protein.mtMAMML.ModelDescription.model_definition; +} + + + +/**************************************** cpREV functions *************************************/ + +/** + * @name models.protein.cpREV.ModelDescription + * @description Create the baseline schema (dictionary) for the cpREV model of protein evolution + * @returns {Dictionary} model description + * @param {String} type + */ +function models.protein.cpREV.ModelDescription(type) { + models.protein.cpREV.ModelDescription.model_definition = models.protein.empirical.ModelDescription(type); + models.protein.cpREV.ModelDescription.model_definition [terms.model.empirical_rates] = models.protein.cpREV.Rij; + models.protein.cpREV.ModelDescription.model_definition [terms.model.frequency_estimator] = "models.protein.cpREV.frequencies"; + return models.protein.cpREV.ModelDescription.model_definition; +} + +/** + * @name models.protein.cpREVF.ModelDescription + * @description Create the baseline schema (dictionary) for the cpREV+F model of protein evolution + * @returns {Dictionary} model description + * @param {String} type + */ +function models.protein.cpREVF.ModelDescription(type) { + models.protein.cpREVF.ModelDescription.model_definition = models.protein.cpREV.ModelDescription(type); + models.protein.cpREVF.ModelDescription.model_definition [terms.model.frequency_estimator] = "frequencies.empirical.protein"; + models.protein.cpREVF.ModelDescription.model_definition [terms.model.efv_estimate_name] = utility.getGlobalValue("terms.frequencies._20x1"); + return models.protein.cpREVF.ModelDescription.model_definition; +} + +/** + * @name models.protein.cpREVML.ModelDescription + * @description Create the baseline schema (dictionary) for the cpREV+ML model of protein evolution + * @returns {Dictionary} model description + * @param {String} type + */ +function models.protein.cpREVML.ModelDescription(type) { + models.protein.cpREVML.ModelDescription.model_definition = models.protein.cpREV.ModelDescription(type); + models.protein.cpREVML.ModelDescription.model_definition [terms.model.frequency_estimator] = "frequencies.ML.protein"; + models.protein.cpREVML.ModelDescription.model_definition [terms.model.efv_estimate_name] = utility.getGlobalValue("terms.frequencies.MLE"); + return models.protein.cpREVML.ModelDescription.model_definition; +} + + +/**************************************** HIVBm functions *************************************/ + + +/** + * @name models.protein.HIVBm.ModelDescription + * @description Create the baseline schema (dictionary) for the HIVBm model of protein evolution + * @returns {Dictionary} model description + * @param {String} type + */ +function models.protein.HIVBm.ModelDescription(type) { + models.protein.HIVBm.ModelDescription.model_definition = models.protein.empirical.ModelDescription(type); + models.protein.HIVBm.ModelDescription.model_definition [terms.model.empirical_rates] = models.protein.HIVBm.Rij; + models.protein.HIVBm.ModelDescription.model_definition [terms.model.frequency_estimator] = "models.protein.HIVBm.frequencies"; + return models.protein.HIVBm.ModelDescription.model_definition; +} + +/** + * @name models.protein.HIVBmF.ModelDescription + * @description Create the baseline schema (dictionary) for the HIVBm+F model of protein evolution + * @returns {Dictionary} model description + * @param {String} type + */ +function models.protein.HIVBmF.ModelDescription(type) { + models.protein.HIVBmF.ModelDescription.model_definition = models.protein.HIVBm.ModelDescription(type); + models.protein.HIVBmF.ModelDescription.model_definition [terms.model.frequency_estimator] = "frequencies.empirical.protein"; + models.protein.HIVBmF.ModelDescription.model_definition [terms.model.efv_estimate_name] = utility.getGlobalValue("terms.frequencies._20x1"); + return models.protein.HIVBmF.ModelDescription.model_definition; +} + +/** + * @name models.protein.HIVBmML.ModelDescription + * @description Create the baseline schema (dictionary) for the HIVBm+ML model of protein evolution + * @returns {Dictionary} model description + * @param {String} type + */ +function models.protein.HIVBmML.ModelDescription(type) { + models.protein.HIVBmML.ModelDescription.model_definition = models.protein.HIVBm.ModelDescription(type); + models.protein.HIVBmML.ModelDescription.model_definition [terms.model.frequency_estimator] = "frequencies.ML.protein"; + models.protein.HIVBmML.ModelDescription.model_definition [terms.model.efv_estimate_name] = utility.getGlobalValue("terms.frequencies.MLE"); + return models.protein.HIVBmML.ModelDescription.model_definition; +} + + + +/**************************************** HIVWm functions *************************************/ + + +/** + * @name models.protein.HIVWm.ModelDescription + * @description Create the baseline schema (dictionary) for the HIVWm model of protein evolution + * @returns {Dictionary} model description + * @param {String} type + */ +function models.protein.HIVWm.ModelDescription(type) { + models.protein.HIVWm.ModelDescription.model_definition = models.protein.empirical.ModelDescription(type); + models.protein.HIVWm.ModelDescription.model_definition [terms.model.empirical_rates] = models.protein.HIVWm.Rij; + models.protein.HIVWm.ModelDescription.model_definition [terms.model.frequency_estimator] = "models.protein.HIVWm.frequencies"; + return models.protein.HIVWm.ModelDescription.model_definition; +} + +/** + * @name models.protein.HIVWmF.ModelDescription + * @description Create the baseline schema (dictionary) for the HIVWm+F model of protein evolution + * @returns {Dictionary} model description + * @param {String} type + */ +function models.protein.HIVWmF.ModelDescription(type) { + models.protein.HIVWmF.ModelDescription.model_definition = models.protein.HIVWm.ModelDescription(type); + models.protein.HIVWmF.ModelDescription.model_definition [terms.model.frequency_estimator] = "frequencies.empirical.protein"; + models.protein.HIVWmF.ModelDescription.model_definition [terms.model.efv_estimate_name] = utility.getGlobalValue("terms.frequencies._20x1"); + return models.protein.HIVWmF.ModelDescription.model_definition; +} + +/** + * @name models.protein.HIVWmML.ModelDescription + * @description Create the baseline schema (dictionary) for the HIVWm+ML model of protein evolution + * @returns {Dictionary} model description + * @param {String} type + */ +function models.protein.HIVWmML.ModelDescription(type) { + models.protein.HIVWmML.ModelDescription.model_definition = models.protein.HIVWm.ModelDescription(type); + models.protein.HIVWmML.ModelDescription.model_definition [terms.model.frequency_estimator] = "frequencies.ML.protein"; + models.protein.HIVWmML.ModelDescription.model_definition [terms.model.efv_estimate_name] = utility.getGlobalValue("terms.frequencies.MLE"); + return models.protein.HIVWmML.ModelDescription.model_definition; +} + + +/**************************************** AB functions *************************************/ + + + +/** + * @name models.protein.AB.ModelDescription + * @description Create the baseline schema (dictionary) for the AB model of protein evolution + * @returns {Dictionary} model description + * @param {String} type + */ +function models.protein.AB.ModelDescription(type) { + models.protein.AB.ModelDescription.model_definition = models.protein.empirical.ModelDescription(type); + models.protein.AB.ModelDescription.model_definition [terms.model.empirical_rates] = models.protein.AB.Rij; + models.protein.AB.ModelDescription.model_definition [terms.model.frequency_estimator] = "models.protein.AB.frequencies"; + return models.protein.AB.ModelDescription.model_definition; +} + +/** + * @name models.protein.ABF.ModelDescription + * @description Create the baseline schema (dictionary) for the AB+F model of protein evolution + * @returns {Dictionary} model description + * @param {String} type + */ +function models.protein.ABF.ModelDescription(type) { + models.protein.ABF.ModelDescription.model_definition = models.protein.AB.ModelDescription(type); + models.protein.ABF.ModelDescription.model_definition [terms.model.frequency_estimator] = "frequencies.empirical.protein"; + models.protein.ABF.ModelDescription.model_definition [terms.model.efv_estimate_name] = utility.getGlobalValue("terms.frequencies._20x1"); + return models.protein.ABF.ModelDescription.model_definition; +} + +/** + * @name models.protein.ABML.ModelDescription + * @description Create the baseline schema (dictionary) for the AB+ML model of protein evolution + * @returns {Dictionary} model description + * @param {String} type + */ +function models.protein.ABML.ModelDescription(type) { + models.protein.ABML.ModelDescription.model_definition = models.protein.AB.ModelDescription(type); + models.protein.ABML.ModelDescription.model_definition [terms.model.frequency_estimator] = "frequencies.ML.protein"; + models.protein.ABML.ModelDescription.model_definition [terms.model.efv_estimate_name] = utility.getGlobalValue("terms.frequencies.MLE"); + return models.protein.ABML.ModelDescription.model_definition; +} + + + /*=============================================================================================*/ /** Below this section are all of the empirical matrices and frequency vectors, including Rij **/ @@ -1007,6 +1266,8 @@ function models.protein.JTT.frequencies (model, namespace, datafilter) { return model; } + + /* Define a dictionary of amino-acid exchangeability rates for the JTT model of protein evolution. */ models.protein.JTT.Rij = { "A":{ @@ -1745,8 +2006,1038 @@ models.protein.mtMAM.Rij = { }; +lfunction models.protein.HIVWm.frequencies (model, namespace, datafilter) { + model[utility.getGlobalValue("terms.efv_estimate")] = + {{0.0377494} + {0.0240105} + {0.0342034} + {0.0618606} + {0.0422741} + {0.0838496} + {0.0156076} + {0.0983641} + {0.0641682} + {0.0577867} + {0.0158419} + {0.0891129} + {0.0458601} + {0.0437824} + {0.057321} + {0.0550846} + {0.0813774} + {0.0515639} + {0.019597} + {0.0205847}}; + model[utility.getGlobalValue("terms.model.efv_estimate_name")] = utility.getGlobalValue("terms.frequencies.predefined"); + (model[utility.getGlobalValue("terms.parameters")])[utility.getGlobalValue("terms.model.empirical")] = 0; + return model; +} + +models.protein.HIVWm.Rij = { + "A":{ + "C":0.167653, + "D":4.43521, + "E":5.56325, + "F":0.597923, + "G":1.8685, + "H":0.005, + "I":0.005, + "K":0.592784, + "L":0.16024, + "M":0.005, + "N":0.617509, + "P":1.00981, + "Q":0.005, + "R":0.0744808, + "S":8.594200000000001, + "T":24.1422, + "V":24.8094, + "W":0.005, + "Y":0.005 + }, + "C":{ + "D":0.005, + "E":0.005, + "F":0.362959, + "G":0.0489798, + "H":0.005, + "I":0.005, + "K":0.005, + "L":0.005, + "M":0.005, + "N":0.0604932, + "P":0.005, + "Q":0.005, + "R":2.86364, + "S":1.12195, + "T":0.005, + "V":0.005, + "W":5.49894, + "Y":8.34835 + }, + "D":{ + "E":12.1233, + "F":0.005, + "G":10.3969, + "H":2.31779, + "I":0.145124, + "K":0.894313, + "L":0.005, + "M":0.005, + "N":29.4087, + "P":0.005, + "Q":0.005, + "R":0.0674539, + "S":0.427881, + "T":0.630395, + "V":2.91786, + "W":0.005, + "Y":2.28154 + }, + "E":{ + "F":0.005, + "G":14.7801, + "H":0.005, + "I":0.0390512, + "K":23.9626, + "L":0.129839, + "M":0.005, + "N":0.201526, + "P":0.005, + "Q":3.20656, + "R":0.0251632, + "S":0.005, + "T":0.458743, + "V":2.19952, + "W":0.005, + "Y":0.005 + }, + "F":{ + "G":0.005, + "H":0.005, + "I":1.48288, + "K":0.005, + "L":7.48781, + "M":0.005, + "N":0.005, + "P":0.0342252, + "Q":0.005, + "R":0.005, + "S":4.27939, + "T":0.114512, + "V":2.28, + "W":0.005, + "Y":4.12728 + }, + "G":{ + "H":0.005, + "I":0.005, + "K":0.279425, + "L":0.0489798, + "M":0.0489798, + "N":0.0604932, + "P":0.005, + "Q":0.0604932, + "R":13.4379, + "S":6.27966, + "T":0.0489798, + "V":2.79622, + "W":2.8258, + "Y":0.005 + }, + "H":{ + "I":0.005, + "K":0.22406, + "L":1.76382, + "M":0.005, + "N":8.59876, + "P":13.9444, + "Q":18.5465, + "R":6.84405, + "S":0.7251570000000001, + "T":0.95956, + "V":0.827479, + "W":0.005, + "Y":47.4889 + }, + "I":{ + "K":0.817481, + "L":9.102460000000001, + "M":17.3064, + "N":0.987028, + "P":0.005, + "Q":0.0342252, + "R":1.34069, + "S":0.7400910000000001, + "T":9.36345, + "V":24.8231, + "W":0.005, + "Y":0.114512 + }, + "K":{ + "L":0.005, + "M":4.09564, + "N":10.6655, + "P":0.111928, + "Q":13.0705, + "R":39.8897, + "S":0.005, + "T":4.04802, + "V":0.128065, + "W":0.005, + "Y":0.005 + }, + "L":{ + "M":11.3839, + "N":0.005, + "P":9.83095, + "Q":2.89048, + "R":0.586757, + "S":6.14396, + "T":0.005, + "V":2.95344, + "W":1.37031, + "Y":0.005 + }, + "M":{ + "N":0.201526, + "P":0.005, + "Q":0.005, + "R":3.28652, + "S":0.392575, + "T":7.41313, + "V":14.7683, + "W":0.005, + "Y":0.579198 + }, + "N":{ + "P":0.344848, + "Q":0.342068, + "R":0.16024, + "S":14.5699, + "T":4.54206, + "V":0.0744808, + "W":0.005, + "Y":5.06475 + }, + "P":{ + "Q":3.04502, + "R":0.404723, + "S":14.249, + "T":4.33701, + "V":0.005, + "W":0.005, + "Y":0.005 + }, + "Q":{ + "R":10.6746, + "S":0.16024, + "T":0.203091, + "V":0.005, + "W":0.0443298, + "Y":0.005 + }, + "R":{ + "S":8.350239999999999, + "T":0.928203, + "V":0.279425, + "W":5.96564, + "Y":0.005 + }, + "S":{ + "T":6.34079, + "V":0.862637, + "W":1.10156, + "Y":0.933142 + }, + "T":{ + "V":0.005, + "W":0.005, + "Y":0.490608 + }, + "V":{ + "W":0.005, + "Y":1.35482 + }, + "W":{ + "Y":0.005 + } +}; + + +lfunction models.protein.HIVBm.frequencies (model, namespace, datafilter) { + model[utility.getGlobalValue("terms.efv_estimate")] = + {{ 0.060490222} + { 0.020075899} + { 0.042109048} + { 0.071567447} + { 0.028809447} + { 0.072308239} + { 0.022293943} + { 0.069730629} + { 0.056968211} + { 0.098851122} + { 0.019768318} + { 0.044127815} + { 0.046025282} + { 0.053606488} + { 0.066039665} + { 0.05060433} + { 0.053636813} + { 0.061625237} + { 0.033011601} + { 0.028350243}}; + model[utility.getGlobalValue("terms.model.efv_estimate_name")] = utility.getGlobalValue("terms.frequencies.predefined"); + (model[utility.getGlobalValue("terms.parameters")])[utility.getGlobalValue("terms.model.empirical")] = 0; + return model; +} + +models.protein.HIVBm.Rij = { + "A":{ + "C":0.123758, + "D":1.45504, + "E":1.48135, + "F":0.0141269, + "G":2.13536, + "H":0.0847613, + "I":0.005, + "K":0.005, + "L":0.215256, + "M":0.0186643, + "N":0.005, + "P":2.12217, + "Q":0.0551128, + "R":0.307507, + "S":2.46633, + "T":15.9183, + "V":7.61428, + "W":0.005, + "Y":0.005 + }, + "C":{ + "D":0.005, + "E":0.005, + "F":9.29815, + "G":0.897871, + "H":0.240073, + "I":0.005, + "K":0.005, + "L":0.129777, + "M":0.005, + "N":0.08606419999999999, + "P":0.005, + "Q":0.005, + "R":0.351721, + "S":4.69314, + "T":0.739969, + "V":0.420027, + "W":2.63277, + "Y":7.57932 + }, + "D":{ + "E":10.5872, + "F":0.005, + "G":2.83806, + "H":1.9169, + "I":0.0176792, + "K":0.005, + "L":0.008760479999999999, + "M":0.005, + "N":17.6612, + "P":0.0342658, + "Q":0.005, + "R":0.005, + "S":0.52823, + "T":0.274724, + "V":1.04793, + "W":0.005, + "Y":0.6746529999999999 + }, + "E":{ + "F":0.005, + "G":3.92775, + "H":0.11974, + "I":0.00609079, + "K":4.61482, + "L":0.005, + "M":0.175789, + "N":0.07926329999999999, + "P":0.0120226, + "Q":2.5602, + "R":0.0749218, + "S":0.005, + "T":0.289774, + "V":1.02847, + "W":0.005, + "Y":0.07926329999999999 + }, + "F":{ + "G":0.291561, + "H":0.145558, + "I":3.39836, + "K":0.0342658, + "L":8.524839999999999, + "M":0.188025, + "N":0.005, + "P":0.005, + "Q":0.005, + "R":0.005, + "S":0.956472, + "T":0.0141269, + "V":0.723274, + "W":0.8293430000000001, + "Y":15.34 + }, + "G":{ + "H":0.005, + "I":0.005, + "K":0.521705, + "L":0.005, + "M":0.005, + "N":0.323401, + "P":0.005, + "Q":0.0619137, + "R":3.65345, + "S":4.38041, + "T":0.369615, + "V":0.953155, + "W":1.21674, + "Y":0.005 + }, + "H":{ + "I":0.103111, + "K":0.005, + "L":1.74171, + "M":0.005, + "N":7.64585, + "P":2.45318, + "Q":7.05545, + "R":9.04044, + "S":0.382747, + "T":0.7115939999999999, + "V":0.005, + "W":0.06951789999999999, + "Y":18.6943 + }, + "I":{ + "K":0.322319, + "L":5.95879, + "M":11.2065, + "N":0.680565, + "P":0.0410593, + "Q":0.005, + "R":0.677289, + "S":1.21803, + "T":8.612170000000001, + "V":17.7389, + "W":0.005, + "Y":0.148168 + }, + "K":{ + "L":0.0814995, + "M":1.28246, + "N":7.90443, + "P":0.0313862, + "Q":6.54737, + "R":20.45, + "S":0.504111, + "T":4.67142, + "V":0.265829, + "W":0.005, + "Y":0.005 + }, + "L":{ + "M":5.31961, + "N":0.005, + "P":2.07757, + "Q":1.49456, + "R":0.701427, + "S":0.927656, + "T":0.0437673, + "V":1.41036, + "W":0.748843, + "Y":0.111986 + }, + "M":{ + "N":0.005, + "P":0.005, + "Q":0.303676, + "R":2.51394, + "S":0.005, + "T":4.94026, + "V":6.8532, + "W":0.089078, + "Y":0.005 + }, + "N":{ + "P":0.00739578, + "Q":0.672052, + "R":0.295543, + "S":13.1447, + "T":6.88667, + "V":0.026656, + "W":0.005, + "Y":1.76417 + }, + "P":{ + "Q":4.47211, + "R":1.28355, + "S":5.37762, + "T":2.01417, + "V":0.005, + "W":0.0444506, + "Y":0.0304381 + }, + "Q":{ + "R":3.4215, + "S":0.116311, + "T":0.243589, + "V":0.0209153, + "W":0.026656, + "Y":0.113033 + }, + "R":{ + "S":3.4791, + "T":2.86868, + "V":0.0812454, + "W":0.9913380000000001, + "Y":0.00991826 + }, + "S":{ + "T":8.93107, + "V":0.0749218, + "W":0.0248728, + "Y":0.648024 + }, + "T":{ + "V":0.709226, + "W":0.005, + "Y":0.105652 + }, + "V":{ + "W":0.005, + "Y":0.0410593 + }, + "W":{ + "Y":1.28022 + } +}; +lfunction models.protein.cpREV.frequencies (model, namespace, datafilter) { + model[utility.getGlobalValue("terms.efv_estimate")] = + {{ 0.076} + { 0.009} + { 0.037} + { 0.049} + { 0.051} + { 0.084} + { 0.025} + { 0.081} + { 0.05} + { 0.101} + { 0.022} + { 0.041} + { 0.043} + { 0.038} + { 0.062} + { 0.062} + { 0.054} + { 0.066} + { 0.018} + { 0.031}}; + model[utility.getGlobalValue("terms.model.efv_estimate_name")] = utility.getGlobalValue("terms.frequencies.predefined"); + (model[utility.getGlobalValue("terms.parameters")])[utility.getGlobalValue("terms.model.empirical")] = 0; + return model; +} + +models.protein.cpREV.Rij = { + "A":{ + "C":1.2816975, + "D":0.33527215, + "E":0.95600458, + "F":0.13027718, + "G":1.2740342, + "H":0.1264455, + "I":0.27779692, + "K":0.45213844, + "L":0.37742065, + "M":0.35443056, + "N":0.43489587, + "P":0.93876201, + "Q":0.30078701, + "R":0.20116329, + "S":4.6746517, + "T":2.5672267, + "V":1.8545339, + "W":0.026821772, + "Y":0.10728709 + }, + "C":{ + "D":0.019158408, + "E":0.019158408, + "F":1.3909005, + "G":0.58049978, + "H":0.84488581, + "I":0.53643544, + "K":0.091960361, + "L":0.75867297, + "M":0.30461869, + "N":1.0307224, + "P":0.54601464, + "Q":0.019158408, + "R":1.576737, + "S":4.465825, + "T":1.1035243, + "V":1.1341778, + "W":0.8333907699999999, + "Y":2.8086227 + }, + "D":{ + "E":7.0713686, + "F":0.042148499, + "G":0.8257274, + "H":0.63414332, + "I":0.019158408, + "K":0.78932643, + "L":0.019158408, + "M":0.09004452, + "N":8.4967541, + "P":0.32569294, + "Q":0.7663363399999999, + "R":0.082381156, + "S":1.1303461, + "T":0.50961366, + "V":0.14368806, + "W":0.034485135, + "Y":0.53835128 + }, + "E":{ + "F":0.27779692, + "G":0.72610368, + "H":0.31036622, + "I":0.28354444, + "K":5.0367456, + "L":0.15709895, + "M":0.21649002, + "N":2.0212121, + "P":0.35443056, + "Q":5.9812551, + "R":0.29120781, + "S":1.0881976, + "T":0.70694527, + "V":0.38316817, + "W":0.12069797, + "Y":0.2720494 + }, + "F":{ + "G":0.047896021, + "H":0.24331179, + "I":0.86979174, + "K":0.13794054, + "L":2.4292862, + "M":0.6264799599999999, + "N":0.18583656, + "P":0.082381156, + "Q":0.019158408, + "R":0.10153956, + "S":0.93301449, + "T":0.28354444, + "V":0.60732155, + "W":0.89661351, + "Y":4.5405428 + }, + "G":{ + "H":0.036400976, + "I":0.07663363400000001, + "K":0.50386614, + "L":0.038316817, + "M":0.040232658, + "N":1.2510441, + "P":0.053643544, + "Q":0.25480683, + "R":0.46554933, + "S":1.323846, + "T":0.17625736, + "V":0.17434152, + "W":0.15709895, + "Y":0.019158408 + }, + "H":{ + "I":0.055559384, + "K":0.58433146, + "L":0.1264455, + "M":0.019158408, + "N":2.6917564, + "P":0.29120781, + "Q":2.431202, + "R":1.3698262, + "S":0.58049978, + "T":0.061306907, + "V":0.047896021, + "W":0.13219302, + "Y":3.7761223 + }, + "I":{ + "K":0.66096509, + "L":3.3431423, + "M":3.39487, + "N":0.32186126, + "P":0.22415338, + "Q":0.17625736, + "R":0.26055435, + "S":0.41382162, + "T":1.9924745, + "V":9.190288499999999, + "W":0.080465315, + "Y":0.17050984 + }, + "K":{ + "L":0.4176533, + "M":0.36975728, + "N":4.6554933, + "P":0.57858393, + "Q":6.3471807, + "R":8.586798699999999, + "S":1.6629499, + "T":1.7587419, + "V":0.47704437, + "W":0.019158408, + "Y":0.47321269 + }, + "L":{ + "M":2.588301, + "N":0.21649002, + "P":0.41956914, + "Q":0.54793048, + "R":0.38891569, + "S":0.98857388, + "T":0.29887117, + "V":1.6572023, + "W":0.30461869, + "Y":0.36209392 + }, + "M":{ + "N":0.11686629, + "P":0.19158408, + "Q":0.38699985, + "R":0.23948011, + "S":0.1781732, + "T":1.2357173, + "V":0.9100244, + "W":0.16476231, + "Y":0.41190578 + }, + "N":{ + "P":0.33144047, + "Q":1.4713658, + "R":0.68395518, + "S":3.9945282, + "T":2.6687663, + "V":0.15901479, + "W":0.07663363400000001, + "Y":1.444544 + }, + "P":{ + "Q":0.61881659, + "R":0.16667815, + "S":2.3028407, + "T":0.49811862, + "V":0.23373258, + "W":0.09387620100000001, + "Y":0.18583656 + }, + "Q":{ + "R":3.3431423, + "S":0.75867297, + "T":0.46171764, + "V":0.10345541, + "W":0.10153956, + "Y":0.74909377 + }, + "R":{ + "S":0.73759872, + "T":0.60157402, + "V":0.17625736, + "W":0.44064339, + "Y":0.61881659 + }, + "S":{ + "T":4.1209737, + "V":0.31994542, + "W":0.13985638, + "Y":1.0000689 + }, + "T":{ + "V":1.456039, + "W":0.055559384, + "Y":0.1360247 + }, + "V":{ + "W":0.019158408, + "Y":0.22798506 + }, + "W":{ + "Y":0.66288093 + } +}; +lfunction models.protein.AB.frequencies (model, namespace, datafilter) { + model[utility.getGlobalValue("terms.efv_estimate")] = + {{6.541704e-02} + {4.708366e-02} + {3.168984e-02} + {4.688141e-02} + {2.150693e-02} + {4.240711e-02} + {2.842211e-02} + {1.005278e-01} + {9.812606e-03} + {3.424424e-02} + {6.222565e-02} + {4.844488e-02} + {1.760370e-02} + {3.478555e-02} + {3.962469e-02} + {1.280566e-01} + {8.199314e-02} + {3.393045e-02} + {7.586119e-02} + {4.948141e-02}}; + model[utility.getGlobalValue("terms.model.efv_estimate_name")] = utility.getGlobalValue("terms.frequencies.predefined"); + (model[utility.getGlobalValue("terms.parameters")])[utility.getGlobalValue("terms.model.empirical")] = 0; + return model; +} + +models.protein.AB.Rij = { + "A":{ + "C":0.1784266, + "D":0.09291290000000001, + "E":1.241095, + "F":0.008929181, + "G":0.1992269, + "H":0.9521821, + "I":1.851951, + "K":5.241316, + "L":0.1140412, + "M":0.06969101, + "N":0.07388355000000001, + "P":0.06299270999999999, + "Q":0.1130146, + "R":1.800713, + "S":0.9988358000000001, + "T":2.912317, + "V":0.07939549, + "W":0.1433528, + "Y":3.774477 + }, + "C":{ + "D":0.782913, + "E":0.05795374, + "F":0.1821885, + "G":1.923901, + "H":0.06273863, + "I":1.0894, + "K":10.4955, + "L":0.3245175, + "M":0.3932002, + "N":7.54924, + "P":0.3362326, + "Q":0.08208677, + "R":0.3498713, + "S":1.926435, + "T":1.135258, + "V":0.5724286, + "W":0.1711315, + "Y":0.1366145 + }, + "D":{ + "E":7.185182, + "F":1.374268e-06, + "G":0.08705989, + "H":0.5038373, + "I":0.4868901, + "K":14.05444, + "L":1.762721, + "M":0.02769442, + "N":6.190318, + "P":0.03972173, + "Q":0.1955446, + "R":0.007342554, + "S":7.348346, + "T":2.147175, + "V":0.0007310937, + "W":2.622763, + "Y":0.04931206 + }, + "E":{ + "F":0.02340019, + "G":0.1843856, + "H":7.426619, + "I":2.1124, + "K":11.26995, + "L":0.03916999, + "M":0.03020502, + "N":0.06622772, + "P":0.03357577, + "Q":0.1031734, + "R":0.1509482, + "S":0.5822988, + "T":0.1516881, + "V":0.01423897, + "W":0.9078338, + "Y":0.4076074 + }, + "F":{ + "G":1.046446e-08, + "H":7.519215e-11, + "I":0.05891123, + "K":3.963388, + "L":0.0006594967, + "M":0.006079219, + "N":3.722878e-16, + "P":0.007213178, + "Q":0.1993818, + "R":0.004878395, + "S":0.2639482, + "T":3.225214e-06, + "V":0.4440833, + "W":0.7741612, + "Y":0.02243512 + }, + "G":{ + "H":3.691671, + "I":0.0551634, + "K":8.908434, + "L":6.712736e-06, + "M":0.6802781, + "N":3.030805, + "P":0.001233336, + "Q":0.00149661, + "R":0.7426909, + "S":0.0005906405, + "T":0.1202094, + "V":4.332983e-05, + "W":0.02737091, + "Y":0.009047737 + }, + "H":{ + "I":1.38937, + "K":7.29808, + "L":0.0001029959, + "M":0.001283121, + "N":3.608816, + "P":0.07659566, + "Q":0.05288625, + "R":0.02889815, + "S":0.06776709, + "T":0.06016624, + "V":0.02252612, + "W":0.1240642, + "Y":0.5795409 + }, + "I":{ + "K":9.139518000000001, + "L":0.03560482, + "M":0.02157936, + "N":0.055044, + "P":0.02187264, + "Q":0.1984772, + "R":0.07915055999999999, + "S":0.9984215, + "T":0.07862767, + "V":0.1386853, + "W":0.2295842, + "Y":0.42282 + }, + "K":{ + "L":4.706586, + "M":5.879103, + "N":1.455741, + "P":2.298295, + "Q":5.642309, + "R":10.49496, + "S":5.439116, + "T":3.443285, + "V":7.01389, + "W":20.55414, + "Y":6.890244 + }, + "L":{ + "M":1.601123, + "N":0.5059793, + "P":10.96748, + "Q":2.714705, + "R":0.05016568, + "S":0.6007607, + "T":3.087152, + "V":0.06318748, + "W":0.2903165, + "Y":7.926675 + }, + "M":{ + "N":0.02158451, + "P":5.647985, + "Q":3.390618, + "R":1.149931, + "S":0.1580539, + "T":0.5702792, + "V":0.3378544, + "W":0.152132, + "Y":3.59531 + }, + "N":{ + "P":1.238634, + "Q":0.004649035, + "R":0.009948993999999999, + "S":0.08688405, + "T":1.039298, + "V":0.008024263, + "W":0.07109973, + "Y":0.0349344 + }, + "P":{ + "Q":3.94794, + "R":0.07417279, + "S":0.01861354, + "T":1.415612, + "V":0.1011149, + "W":0.002246759, + "Y":4.39672 + }, + "Q":{ + "R":0.3556198, + "S":0.9813064, + "T":0.03674486, + "V":0.2199856, + "W":7.074464, + "Y":1.643946 + }, + "R":{ + "S":1.284651, + "T":0.9057112, + "V":0.005516074, + "W":0.1992133, + "Y":0.2217442 + }, + "S":{ + "T":3.058575, + "V":0.1385142, + "W":0.8104751, + "Y":0.07477041 + }, + "T":{ + "V":0.01412361, + "W":0.09984255, + "Y":0.2166054 + }, + "V":{ + "W":0.6121284, + "Y":0.09663569 + }, + "W":{ + "Y":0.5010635 + } +}; diff --git a/res/TemplateBatchFiles/libv3/models/protein/matrix2dict_aux.bf b/res/TemplateBatchFiles/libv3/models/protein/matrix2dict_aux.bf index 4671bf672..b5507c575 100644 --- a/res/TemplateBatchFiles/libv3/models/protein/matrix2dict_aux.bf +++ b/res/TemplateBatchFiles/libv3/models/protein/matrix2dict_aux.bf @@ -18,7 +18,7 @@ function mapMatrixToDict (letters, matrix) { return dict; } -models.protein.WAG.empirical_Q = { +WAG = { {0, 3.26324, 2.34804, 5.02923, 0.668808, 4.50138, 1.00707, 0.6142879999999999, 2.8795, 1.26431, 2.83893, 1.61995, 4.57074, 2.88691, 1.75252, 10.7101, 6.73946, 6.37375, 0.35946, 0.764894} {3.26324, 0, 0.0962568, 0.06784229999999999, 1.26464, 0.974403, 0.791065, 0.540574, 0.23523, 1.22101, 1.24069, 0.842805, 0.347612, 0.313977, 1.67824, 4.4726, 1.62992, 3.18413, 2.27837, 1.72794} {2.34804, 0.0962568, 0, 19.6173, 0.148478, 2.75024, 2.95706, 0.125304, 1.52466, 0.269452, 0.32966, 17.251, 1.34714, 1.95972, 0.468033, 3.40533, 1.19107, 0.484018, 0.412312, 1.03489} @@ -41,12 +41,12 @@ models.protein.WAG.empirical_Q = { {0.764894, 1.72794, 1.03489, 0.623719, 20.5074, 0.329184, 12.3072, 1.33502, 0.423423, 1.26654, 1.36128, 3.45058, 0.686449, 0.723509, 1.21225, 2.50053, 0.925072, 1, 7.8969, 0} }; -//fprintf (stdout, mapMatrixToDict (models.protein.alphabet, models.protein.WAG.empirical_Q), "\n"): +//fprintf (stdout, mapMatrixToDict (models.protein.alphabet, WAG), "\n"): -models.protein.LG.empirical_Q = { +LG = { {0, 7.908630, 1.255501, 3.299796, 0.806091, 6.564482, 1.140209, 0.476059, 1.704692, 1.256114, 3.571425, 0.879541, 3.741781, 3.081669, 1.350659, 15.019796, 6.797891, 8.095413, 0.574197, 0.695704} {7.908630, 0, 0.198761, 0.011117, 3.511742, 1.808740, 2.035214, 1.018736, 0.042150, 1.887354, 2.839512, 1.680068, 0.239513, 0.269463, 1.698443, 8.847193, 3.633208, 6.225305, 2.129215, 3.703275} {1.255501, 0.198761, 0, 16.661482, 0.055336, 2.684605, 2.945743, 0.033966, 0.899053, 0.047901, 0.081174, 16.128578, 1.253315, 1.662968, 0.393842, 3.940757, 1.353096, 0.120634, 0.094970, 0.429279} @@ -68,10 +68,10 @@ models.protein.LG.empirical_Q = { {0.574197, 2.129215, 0.094970, 0.247361, 7.807073, 0.853083, 1.897035, 0.354780, 0.158568, 1.968773, 2.211975, 0.144174, 0.302262, 0.750481, 1.886083, 0.790716, 0.447447, 0.602135, 0, 10.014342} {0.695704, 3.703275, 0.429279, 0.381397, 24.795538, 0.173733, 16.861539, 0.738801, 0.419191, 0.952079, 1.529266, 1.944603, 0.284730, 0.817640, 0.999078, 1.272668, 0.781117, 0.792149, 10.014342, 0} }; -//fprintf (stdout, mapMatrixToDict (models.protein.alphabet, models.protein.LG.empirical_Q), "\n"): +//fprintf (stdout, mapMatrixToDict (models.protein.alphabet, LG), "\n"): -models.protein.JTT.empirical_Q = { +JTT = { {0, 1.825304, 2.629062, 3.389193, 0.439402, 5.529052, 0.698916, 1.149188, 1.173822, 0.984993, 1.491421, 1.772843, 6.226284, 1.768897, 1.689314, 12.350566, 14.560301, 9.291012, 0.267941, 0.443212} {1.825304, 0, 0.335605, 0.171280, 2.155291, 1.736056, 2.303555, 0.478375, 0.155718, 0.522966, 1.300168, 0.995491, 0.392886, 0.290103, 3.240373, 6.848188, 1.492781, 1.974145, 3.508343, 6.719573} {2.629062, 0.335605, 0, 24.676880, 0.103333, 4.042937, 3.280087, 0.368468, 0.897486, 0.195361, 0.603695, 17.632664, 0.404041, 1.657439, 0.492165, 1.872296, 1.350869, 1.001687, 0.182588, 1.442353} @@ -94,89 +94,83 @@ models.protein.JTT.empirical_Q = { {0.443212, 6.719573, 1.442353, 0.201608, 17.425203, 0.166810, 18.582271, 0.964144, 0.279300, 0.766034, 0.603279, 2.226330, 0.361739, 0.809408, 0.748581, 1.997292, 0.638941, 0.525762, 2.376287, 0} }; -//fprintf (stdout, mapMatrixToDict (models.protein.alphabet, models.protein.JTT.empirical_Q), "\n"): - - - -/*****************************************************************************************************************************************************************************************************************************/ - -models.protein.WAG.empirical_Q_normalized = { - {0., 2.081175e-02, 4.424356e-02, 9.644885e-02, 8.490242e-03, 1.237844e-01, 8.127019e-03, 9.834135e-03, 5.899778e-02, 3.600239e-02, 1.828840e-02, 2.091646e-02, 6.909218e-02, 3.502343e-02, 2.545459e-02, 2.459330e-01, 1.358225e-01, 1.492591e-01, 1.708106e-03, 8.912199e-03} - {9.337576e-02, 0., 1.813745e-03, 1.301055e-03, 1.605407e-02, 2.679532e-02, 6.383892e-03, 8.654049e-03, 4.819601e-03, 3.476936e-02, 7.992529e-03, 1.088210e-02, 5.254569e-03, 3.809101e-03, 2.437562e-02, 1.027029e-01, 3.284827e-02, 7.456520e-02, 1.082647e-02, 2.013312e-02} - {6.718774e-02, 6.138903e-04, 0., 3.762142e-01, 1.884863e-03, 7.562952e-02, 2.386347e-02, 2.005993e-03, 3.123852e-02, 7.672926e-03, 2.123675e-03, 2.227414e-01, 2.036354e-02, 2.377493e-02, 6.797970e-03, 7.819566e-02, 2.400406e-02, 1.133463e-02, 1.959249e-03, 1.205807e-02} - {1.439085e-01, 4.326730e-04, 3.696449e-01, 0., 3.272523e-03, 4.960369e-02, 1.461601e-02, 6.480045e-03, 1.682461e-01, 1.395734e-02, 6.450074e-03, 3.885870e-02, 3.277285e-02, 2.108299e-01, 2.026676e-02, 5.143238e-02, 5.268470e-02, 4.380510e-02, 2.363730e-03, 7.267292e-03} - {1.913755e-02, 8.065405e-03, 2.797733e-03, 4.943786e-03, 0., 4.362670e-03, 1.741975e-02, 5.389076e-02, 5.783216e-03, 1.913755e-01, 2.437025e-02, 3.945040e-03, 7.754001e-03, 3.851615e-03, 4.740036e-03, 3.983115e-02, 1.100758e-02, 4.835584e-02, 2.309483e-02, 2.389425e-01} - {1.288044e-01, 6.214386e-03, 5.182222e-02, 3.459307e-02, 2.013959e-03, 0., 6.395123e-03, 1.548868e-03, 2.431859e-02, 5.546612e-03, 3.563543e-03, 4.617598e-02, 1.169843e-02, 1.272240e-02, 2.698185e-02, 9.789926e-02, 1.446092e-02, 1.393229e-02, 5.087841e-03, 3.835501e-03} - {2.881661e-02, 5.045123e-03, 5.571926e-02, 3.473371e-02, 2.740232e-02, 2.179194e-02, 0., 7.029141e-03, 5.796705e-02, 4.519012e-02, 8.272107e-03, 1.623063e-01, 3.343772e-02, 1.655236e-01, 9.862788e-02, 5.400277e-02, 3.030760e-02, 8.806541e-03, 3.964322e-03, 1.433978e-01} - {1.757750e-02, 3.447584e-03, 2.361080e-03, 7.762644e-03, 4.273355e-02, 2.660546e-03, 3.543330e-03, 0., 2.108143e-02, 2.869017e-01, 8.714326e-02, 2.273747e-02, 4.799484e-03, 4.391122e-03, 8.628942e-03, 2.330635e-02, 9.337141e-02, 5.819514e-01, 3.208113e-03, 1.555502e-02} - {8.239521e-02, 1.500209e-03, 2.872875e-02, 1.574787e-01, 3.583186e-03, 3.263925e-02, 2.283157e-02, 1.647196e-02, 0., 2.330296e-02, 1.912310e-02, 1.235674e-01, 2.674718e-02, 1.501354e-01, 2.469640e-01, 7.056185e-02, 8.881349e-02, 2.272611e-02, 2.076079e-03, 4.933538e-03} - {3.617737e-02, 7.787122e-03, 5.077229e-03, 9.399810e-03, 8.531504e-02, 5.356348e-03, 1.280670e-02, 1.612938e-01, 1.676681e-02, 0., 9.935387e-02, 5.395922e-03, 1.997258e-02, 3.351591e-02, 2.296714e-02, 2.515217e-02, 2.091482e-02, 1.339560e-01, 1.004497e-02, 1.475715e-02} - {8.123429e-02, 7.912656e-03, 6.211717e-03, 1.920166e-02, 4.802387e-02, 1.521181e-02, 1.036258e-02, 2.165590e-01, 6.082129e-02, 4.391801e-01, 0., 8.131996e-03, 8.228768e-03, 5.956464e-02, 3.152741e-02, 3.603533e-02, 9.708280e-02, 1.531609e-01, 7.786238e-03, 1.586107e-02} - {4.635403e-02, 5.375099e-03, 3.250575e-01, 5.771624e-02, 3.878683e-03, 9.834466e-02, 1.014432e-01, 2.819164e-02, 1.960816e-01, 1.190033e-02, 4.057260e-03, 0., 9.369553e-03, 5.950219e-02, 2.932074e-02, 2.899600e-01, 1.299922e-01, 1.460187e-02, 1.085813e-03, 4.020457e-02} - {1.307891e-01, 2.216942e-03, 2.538378e-02, 4.157839e-02, 6.511818e-03, 2.128168e-02, 1.785121e-02, 5.082956e-03, 3.625388e-02, 3.762457e-02, 3.506825e-03, 8.003178e-03, 0., 3.597839e-02, 3.135791e-02, 1.177049e-01, 5.093139e-02, 2.342947e-02, 2.104766e-03, 7.998193e-03} - {8.260732e-02, 2.002428e-03, 3.692659e-02, 3.332748e-01, 4.030289e-03, 2.883795e-02, 1.101053e-01, 5.794476e-03, 2.535576e-01, 7.866928e-02, 3.162895e-02, 6.332767e-02, 4.482896e-02, 0., 1.400860e-01, 7.506641e-02, 5.493632e-02, 2.241710e-02, 3.257243e-03, 8.430004e-03} - {5.014737e-02, 1.070316e-02, 8.819041e-03, 2.675944e-02, 4.142831e-03, 5.108451e-02, 5.479867e-02, 9.510832e-03, 3.483770e-01, 4.502808e-02, 1.398321e-02, 2.606500e-02, 3.263520e-02, 1.170084e-01, 0., 8.931696e-02, 3.550112e-02, 1.873906e-02, 1.757311e-02, 1.412465e-02} - {3.064633e-01, 2.852456e-02, 6.416591e-02, 4.295451e-02, 2.202004e-02, 1.172401e-01, 1.897867e-02, 1.624856e-02, 6.296009e-02, 3.119116e-02, 1.010943e-02, 1.630423e-01, 7.748429e-02, 3.965952e-02, 5.649546e-02, 0., 2.803409e-01, 1.731717e-02, 7.907568e-03, 2.913510e-02} - {1.928457e-01, 1.039501e-02, 2.244310e-02, 5.013408e-02, 6.933680e-03, 1.973192e-02, 1.213606e-02, 7.417044e-02, 9.029229e-02, 2.955197e-02, 3.103250e-02, 8.328300e-02, 3.820153e-02, 3.307027e-02, 2.558574e-02, 3.194205e-01, 0., 1.032926e-01, 1.673848e-03, 1.077852e-02} - {1.823811e-01, 2.030718e-02, 9.120245e-03, 3.587353e-02, 2.621329e-02, 1.636051e-02, 3.034818e-03, 3.978365e-01, 1.988373e-02, 1.628904e-01, 4.213311e-02, 8.050972e-03, 1.512372e-02, 1.161338e-02, 1.162264e-02, 1.698065e-02, 8.889353e-02, 0., 5.516418e-03, 1.165155e-02} - {1.028575e-02, 1.453057e-02, 7.769107e-03, 9.539592e-03, 6.169778e-02, 2.944354e-02, 6.732533e-03, 1.080811e-02, 8.951565e-03, 6.019556e-02, 1.055566e-02, 2.950375e-03, 6.695488e-03, 8.315944e-03, 5.371404e-02, 3.821224e-02, 7.099034e-03, 2.718563e-02, 0., 9.201110e-02} - {2.188698e-02, 1.102013e-02, 1.950020e-02, 1.196146e-02, 2.603323e-01, 9.052294e-03, 9.931890e-02, 2.137227e-02, 8.675476e-03, 3.606600e-02, 8.769406e-03, 4.455303e-02, 1.037648e-02, 8.777464e-03, 1.760746e-02, 5.741905e-02, 1.864329e-02, 2.341779e-02, 3.752494e-02, 0.} +//fprintf (stdout, mapMatrixToDict (models.protein.alphabet, JTT), "\n"): + + + +HIVBm = { +{0,0.123758,1.45504,1.48135,0.0141269,2.13536,0.0847613,0.005,0.005,0.215256,0.0186643,0.005,2.12217,0.0551128,0.307507,2.46633,15.9183,7.61428,0.005,0.005} +{0.123758,0,0.005,0.005,9.29815,0.897871,0.240073,0.005,0.005,0.129777,0.005,0.0860642,0.005,0.005,0.351721,4.69314,0.739969,0.420027,2.63277,7.57932} +{1.45504,0.005,0,10.5872,0.005,2.83806,1.9169,0.0176792,0.005,0.00876048,0.005,17.6612,0.0342658,0.005,0.005,0.52823,0.274724,1.04793,0.005,0.674653} +{1.48135,0.005,10.5872,0,0.005,3.92775,0.11974,0.00609079,4.61482,0.005,0.175789,0.0792633,0.0120226,2.5602,0.0749218,0.005,0.289774,1.02847,0.005,0.0792633} +{0.0141269,9.29815,0.005,0.005,0,0.291561,0.145558,3.39836,0.0342658,8.52484,0.188025,0.005,0.005,0.005,0.005,0.956472,0.0141269,0.723274,0.829343,15.34} +{2.13536,0.897871,2.83806,3.92775,0.291561,0,0.005,0.005,0.521705,0.005,0.005,0.323401,0.005,0.0619137,3.65345,4.38041,0.369615,0.953155,1.21674,0.005} +{0.0847613,0.240073,1.9169,0.11974,0.145558,0.005,0,0.103111,0.005,1.74171,0.005,7.64585,2.45318,7.05545,9.04044,0.382747,0.711594,0.005,0.0695179,18.6943} +{0.005,0.005,0.0176792,0.00609079,3.39836,0.005,0.103111,0,0.322319,5.95879,11.2065,0.680565,0.0410593,0.005,0.677289,1.21803,8.61217,17.7389,0.005,0.148168} +{0.005,0.005,0.005,4.61482,0.0342658,0.521705,0.005,0.322319,0,0.0814995,1.28246,7.90443,0.0313862,6.54737,20.45,0.504111,4.67142,0.265829,0.005,0.005} +{0.215256,0.129777,0.00876048,0.005,8.52484,0.005,1.74171,5.95879,0.0814995,0,5.31961,0.005,2.07757,1.49456,0.701427,0.927656,0.0437673,1.41036,0.748843,0.111986} +{0.0186643,0.005,0.005,0.175789,0.188025,0.005,0.005,11.2065,1.28246,5.31961,0,0.005,0.005,0.303676,2.51394,0.005,4.94026,6.8532,0.089078,0.005} +{0.005,0.0860642,17.6612,0.0792633,0.005,0.323401,7.64585,0.680565,7.90443,0.005,0.005,0,0.00739578,0.672052,0.295543,13.1447,6.88667,0.026656,0.005,1.76417} +{2.12217,0.005,0.0342658,0.0120226,0.005,0.005,2.45318,0.0410593,0.0313862,2.07757,0.005,0.00739578,0,4.47211,1.28355,5.37762,2.01417,0.005,0.0444506,0.0304381} +{0.0551128,0.005,0.005,2.5602,0.005,0.0619137,7.05545,0.005,6.54737,1.49456,0.303676,0.672052,4.47211,0,3.4215,0.116311,0.243589,0.0209153,0.026656,0.113033} +{0.307507,0.351721,0.005,0.0749218,0.005,3.65345,9.04044,0.677289,20.45,0.701427,2.51394,0.295543,1.28355,3.4215,0,3.4791,2.86868,0.0812454,0.991338,0.00991826} +{2.46633,4.69314,0.52823,0.005,0.956472,4.38041,0.382747,1.21803,0.504111,0.927656,0.005,13.1447,5.37762,0.116311,3.4791,0,8.93107,0.0749218,0.0248728,0.648024} +{15.9183,0.739969,0.274724,0.289774,0.0141269,0.369615,0.711594,8.61217,4.67142,0.0437673,4.94026,6.88667,2.01417,0.243589,2.86868,8.93107,0,0.709226,0.005,0.105652} +{7.61428,0.420027,1.04793,1.02847,0.723274,0.953155,0.005,17.7389,0.265829,1.41036,6.8532,0.026656,0.005,0.0209153,0.0812454,0.0749218,0.709226,0,0.005,0.0410593} +{0.005,2.63277,0.005,0.005,0.829343,1.21674,0.0695179,0.005,0.005,0.748843,0.089078,0.005,0.0444506,0.026656,0.991338,0.0248728,0.005,0.005,0,1.28022} +{0.005,7.57932,0.674653,0.0792633,15.34,0.005,18.6943,0.148168,0.005,0.111986,0.005,1.76417,0.0304381,0.113033,0.00991826,0.648024,0.105652,0.0410593,1.28022,0} }; -fprintf(stdout, "\n\n\nWAG"); -fprintf (stdout, mapMatrixToDict (models.protein.alphabet, models.protein.WAG.empirical_Q_normalized), "\n"); - - -models.protein.LG.empirical_Q_normalized = { - {0., 3.522956e-02, 2.293460e-02, 8.133689e-02, 1.174132e-02, 1.296008e-01, 8.776705e-03, 1.018879e-02, 3.791848e-02, 4.285406e-02, 2.822381e-02, 1.271276e-02, 5.674114e-02, 4.325807e-02, 2.601647e-02, 3.164948e-01, 1.247291e-01, 1.927457e-01, 2.385593e-03, 8.181845e-03} - {2.153069e-01, 0., 3.630821e-03, 2.740351e-04, 5.115122e-02, 3.570948e-02, 1.566596e-02, 2.180340e-02, 9.375764e-04, 6.438966e-02, 2.243974e-02, 2.428347e-02, 3.632027e-03, 3.782507e-03, 3.271550e-02, 1.864267e-01, 6.666287e-02, 1.482198e-01, 8.846170e-03, 4.355245e-02} - {3.418014e-02, 8.853943e-04, 0., 4.106900e-01, 8.060157e-04, 5.300146e-02, 2.267472e-02, 7.269456e-04, 1.999816e-02, 1.634220e-03, 6.414941e-04, 2.331202e-01, 1.900553e-02, 2.334345e-02, 7.586211e-03, 8.303904e-02, 2.482688e-02, 2.872194e-03, 3.945694e-04, 5.048546e-03} - {8.983463e-02, 4.952354e-05, 3.043602e-01, 0., 8.705765e-04, 2.188286e-02, 1.036699e-02, 3.010126e-03, 1.277224e-01, 7.552471e-03, 4.362376e-03, 2.487791e-02, 2.020781e-02, 1.841385e-01, 2.227563e-02, 4.097288e-02, 3.524391e-02, 1.853676e-02, 1.027702e-03, 4.485425e-03} - {2.194525e-02, 1.564331e-02, 1.010844e-03, 1.473242e-03, 0., 5.619650e-03, 1.668329e-02, 7.566810e-02, 1.690408e-03, 2.810447e-01, 4.516806e-02, 4.111401e-03, 4.551429e-03, 1.599162e-03, 3.226682e-03, 2.422454e-02, 9.619268e-03, 4.952661e-02, 3.243575e-02, 2.916085e-01} - {1.787134e-01, 8.057165e-03, 4.904047e-02, 2.732104e-02, 4.146056e-03, 0., 7.618063e-03, 5.919608e-04, 2.096479e-02, 4.797840e-03, 3.503711e-03, 6.602330e-02, 9.489902e-03, 1.195119e-02, 2.388046e-02, 1.164960e-01, 7.569210e-03, 5.802412e-03, 3.544273e-03, 2.043191e-03} - {3.104138e-02, 9.066007e-03, 5.381075e-02, 3.319756e-02, 3.156951e-02, 1.953911e-02, 0., 7.404237e-03, 4.927923e-02, 3.970833e-02, 1.111019e-02, 2.070850e-01, 2.451727e-02, 2.146863e-01, 1.485124e-01, 6.628340e-02, 3.406144e-02, 9.003304e-03, 7.881540e-03, 1.983005e-01} - {1.296036e-02, 4.538035e-03, 6.204598e-04, 3.466751e-03, 5.149721e-02, 5.460569e-04, 2.662961e-03, 0., 1.124222e-02, 4.493203e-01, 1.073075e-01, 8.794702e-03, 3.771706e-03, 3.249348e-03, 7.772081e-03, 4.291965e-03, 6.026516e-02, 4.910478e-02, 1.473992e-03, 8.688692e-03} - {4.640905e-02, 1.877620e-04, 1.642326e-02, 1.415347e-01, 1.106929e-03, 1.860771e-02, 1.705320e-02, 1.081707e-02, 0., 1.490483e-02, 1.648691e-02, 9.851189e-02, 1.880635e-02, 1.442522e-01, 3.871668e-01, 5.012591e-02, 6.627711e-02, 1.401048e-02, 6.587949e-04, 4.929905e-03} - {3.419683e-02, 8.407354e-03, 8.750282e-04, 5.456658e-03, 1.199902e-01, 2.776453e-03, 8.959132e-03, 2.818745e-01, 9.717832e-03, 0., 1.584993e-01, 3.142484e-03, 1.200011e-02, 2.597806e-02, 1.847365e-02, 1.220450e-02, 1.766063e-02, 1.288122e-01, 8.179586e-03, 1.119695e-02} - {9.722955e-02, 1.264881e-02, 1.482835e-03, 1.360659e-02, 8.325125e-02, 8.753095e-03, 1.082167e-02, 2.906155e-01, 4.640558e-02, 6.842521e-01, 0., 1.703821e-02, 4.810887e-03, 7.459796e-02, 2.962982e-02, 2.322970e-02, 1.177837e-01, 1.436375e-01, 9.190009e-03, 1.798497e-02} - {2.394488e-02, 7.483985e-03, 2.946254e-01, 4.242586e-02, 4.143233e-03, 9.018219e-02, 1.102838e-01, 1.302266e-02, 1.516037e-01, 7.417406e-03, 9.315676e-03, 0., 7.795161e-03, 7.563194e-02, 4.601631e-02, 2.683680e-01, 1.166360e-01, 6.330976e-03, 5.989957e-04, 2.286955e-02} - {1.018673e-01, 1.066928e-03, 2.289467e-02, 3.284732e-02, 4.371811e-03, 1.235518e-02, 1.244513e-02, 5.323295e-03, 2.758606e-02, 2.699781e-02, 2.507145e-03, 7.430006e-03, 0., 2.784403e-02, 2.035162e-02, 8.959077e-02, 3.331558e-02, 2.243022e-02, 1.255797e-03, 3.348570e-03} - {8.389628e-02, 1.200340e-03, 3.037792e-02, 3.233435e-01, 1.659376e-03, 1.680883e-02, 1.177254e-01, 4.954246e-03, 2.285842e-01, 6.313765e-02, 4.199715e-02, 7.787677e-02, 3.007950e-02, 0., 1.718491e-01, 8.193787e-02, 6.297003e-02, 1.591156e-02, 3.117996e-03, 9.615879e-03} - {3.677074e-02, 7.565835e-03, 7.194431e-03, 2.850544e-02, 2.439983e-03, 2.447639e-02, 5.934815e-02, 8.635692e-03, 4.470957e-01, 3.271996e-02, 1.215627e-02, 3.452971e-02, 1.602198e-02, 1.252350e-01, 0., 5.745502e-02, 3.375392e-02, 1.292756e-02, 7.836037e-03, 1.174968e-02} - {4.089034e-01, 3.941046e-02, 7.198697e-02, 4.792857e-02, 1.674505e-02, 1.091480e-01, 2.421304e-02, 4.359293e-03, 5.291328e-02, 1.975969e-02, 8.711946e-03, 1.840823e-01, 6.447338e-02, 5.458374e-02, 5.252041e-02, 0., 3.773224e-01, 7.441590e-03, 3.285156e-03, 1.496724e-02} - {1.850678e-01, 1.618439e-02, 2.471740e-02, 4.734683e-02, 7.636277e-03, 8.144497e-03, 1.428948e-02, 7.029673e-02, 8.034795e-02, 3.283790e-02, 5.073011e-02, 9.188042e-02, 2.753426e-02, 4.817496e-02, 3.543506e-02, 4.333327e-01, 0., 1.655336e-01, 1.858990e-03, 9.186346e-03} - {2.203919e-01, 2.773107e-02, 2.203648e-03, 1.919060e-02, 3.029885e-02, 4.811385e-03, 2.910739e-03, 4.414083e-02, 1.308917e-02, 1.845755e-01, 4.767561e-02, 3.843339e-03, 1.428590e-02, 9.380981e-03, 1.045859e-02, 6.586013e-03, 1.275657e-01, 0., 2.501667e-03, 9.316084e-03} - {1.563210e-02, 9.484742e-03, 1.734849e-03, 6.097222e-03, 1.137160e-01, 1.684220e-02, 1.460234e-02, 7.593148e-03, 3.527113e-03, 6.716738e-02, 1.748051e-02, 2.083875e-03, 4.583566e-03, 1.053467e-02, 3.632983e-02, 1.666183e-02, 8.209850e-03, 1.433638e-02, 0., 1.177739e-01} - {1.894005e-02, 1.649650e-02, 7.841764e-03, 9.401073e-03, 3.611660e-01, 3.429965e-03, 1.297909e-01, 1.581212e-02, 9.324313e-03, 3.248149e-02, 1.208529e-02, 2.810701e-02, 4.317700e-03, 1.147739e-02, 1.924430e-02, 2.681747e-02, 1.433209e-02, 1.886046e-02, 4.160621e-02, 0.} -}; -fprintf(stdout, "\n\n\nLG"); - -fprintf (stdout, mapMatrixToDict (models.protein.alphabet, models.protein.LG.empirical_Q_normalized), "\n"); - - -fprintf(stdout, "\n\n\nJTT"); -models.protein.JTT.empirical_Q_normalized = { - {0., 1.164983e-02, 4.242226e-02, 6.594220e-02, 5.605013e-03, 1.300142e-01, 5.055569e-03, 1.901336e-02, 2.198075e-02, 2.824504e-02, 1.099041e-02, 2.373925e-02, 9.902243e-02, 2.285967e-02, 2.714587e-02, 2.651969e-01, 2.681624e-01, 1.940882e-01, 1.208940e-03, 4.506008e-03} - {4.415494e-02, 0., 5.415286e-03, 3.332529e-03, 2.749291e-02, 4.082289e-02, 1.666262e-02, 7.914733e-03, 2.915936e-03, 1.499622e-02, 9.581053e-03, 1.333012e-02, 6.248431e-03, 3.749032e-03, 5.207011e-02, 1.470474e-01, 2.749309e-02, 4.123968e-02, 1.582953e-02, 6.831604e-02} - {6.359823e-02, 2.141968e-03, 0., 4.801284e-01, 1.318116e-03, 9.506860e-02, 2.372630e-02, 6.096320e-03, 1.680615e-02, 5.602049e-03, 4.448682e-03, 2.361102e-01, 6.425849e-03, 2.141930e-02, 7.908676e-03, 4.020279e-02, 2.487944e-02, 2.092512e-02, 8.238323e-04, 1.466400e-02} - {8.198614e-02, 1.093179e-03, 3.981835e-01, 0., 1.776388e-03, 8.335330e-02, 5.602518e-03, 5.875793e-03, 1.030317e-01, 8.881953e-03, 4.099415e-03, 2.459647e-02, 9.701838e-03, 1.403343e-01, 1.626078e-02, 2.131682e-02, 1.940362e-02, 3.088188e-02, 1.639765e-03, 2.049689e-03} - {1.062933e-02, 1.375595e-02, 1.667369e-03, 2.709508e-03, 0., 3.751538e-03, 1.042113e-02, 4.085083e-02, 1.458950e-03, 2.278042e-01, 1.021273e-02, 3.126321e-03, 7.503140e-03, 1.875789e-03, 3.334736e-03, 6.440240e-02, 8.128382e-03, 3.939149e-02, 7.711647e-03, 1.771572e-01} - {1.337503e-01, 1.108021e-02, 6.523640e-02, 6.896835e-02, 2.035091e-03, 0., 4.635577e-03, 2.826581e-03, 1.605493e-02, 6.331483e-03, 3.052693e-03, 3.290136e-02, 1.051474e-02, 9.497160e-03, 6.941973e-02, 1.278738e-01, 1.854212e-02, 3.120506e-02, 7.801362e-03, 1.695907e-03} - {1.690710e-02, 1.470223e-02, 5.292712e-02, 1.506973e-02, 1.837743e-02, 1.506951e-02, 0., 9.556410e-03, 3.124215e-02, 4.925195e-02, 7.718657e-03, 1.712807e-01, 5.770555e-02, 2.333939e-01, 1.639271e-01, 5.072240e-02, 2.793385e-02, 8.086143e-03, 1.837774e-03, 1.889208e-01} - {2.779938e-02, 3.053185e-03, 5.945561e-03, 6.909805e-03, 3.149544e-02, 4.017296e-03, 4.178032e-03, 0., 1.205203e-02, 2.127567e-01, 1.131285e-01, 2.089020e-02, 4.981443e-03, 3.213843e-03, 1.221257e-02, 2.763923e-02, 1.494435e-01, 6.328057e-01, 1.928334e-03, 9.802181e-03} - {2.839528e-02, 9.938532e-04, 1.448174e-02, 1.070526e-01, 9.938358e-04, 2.016082e-02, 1.206827e-02, 1.064847e-02, 0., 1.334602e-02, 1.462393e-02, 1.076208e-01, 1.093234e-02, 1.218169e-01, 3.333640e-01, 3.237125e-02, 5.650736e-02, 8.234754e-03, 1.277824e-03, 2.839562e-03} - {2.382744e-02, 3.337780e-03, 3.152324e-03, 6.026521e-03, 1.013368e-01, 5.192023e-03, 1.242394e-02, 1.227558e-01, 8.715324e-03, 0., 9.030557e-02, 5.841096e-03, 5.358937e-02, 2.911240e-02, 1.900652e-02, 4.042405e-02, 1.594697e-02, 1.169137e-01, 7.602722e-03, 7.788057e-03} - {3.607815e-02, 8.298205e-03, 9.741158e-03, 1.082369e-02, 1.767841e-02, 9.741134e-03, 7.576573e-03, 2.539957e-01, 3.716131e-02, 3.514064e-01, 0., 1.407080e-02, 8.298110e-03, 1.876080e-02, 2.200785e-02, 1.948259e-02, 1.237496e-01, 2.016795e-01, 2.886323e-03, 6.133368e-03} - {4.288589e-02, 6.353632e-03, 2.845187e-01, 3.573905e-02, 2.978184e-03, 5.777729e-02, 9.252443e-02, 2.581152e-02, 1.505011e-01, 1.250853e-02, 7.743476e-03, 0., 6.154998e-03, 3.156908e-02, 2.303155e-02, 3.450795e-01, 1.375939e-01, 1.092017e-02, 3.971070e-04, 2.263447e-02} - {1.506166e-01, 2.507558e-03, 6.519569e-03, 1.186906e-02, 6.018014e-03, 1.554655e-02, 2.624568e-02, 5.182250e-03, 1.287209e-02, 9.662355e-02, 3.844928e-03, 5.182271e-03, 0., 6.603124e-02, 3.627542e-02, 1.902389e-01, 6.887338e-02, 1.404214e-02, 1.003017e-03, 3.677695e-03} - {4.279042e-02, 1.851553e-03, 2.674426e-02, 2.112825e-01, 1.851531e-03, 1.728089e-02, 1.306371e-01, 4.114574e-03, 1.765145e-01, 6.459804e-02, 1.069787e-02, 3.271080e-02, 8.126180e-02, 0., 1.542939e-01, 3.744234e-02, 3.065318e-02, 1.193211e-02, 2.468744e-03, 8.229025e-03} - {4.086529e-02, 2.068139e-02, 7.941514e-03, 1.968861e-02, 2.647175e-03, 1.015850e-01, 7.379083e-02, 1.257423e-02, 3.884775e-01, 3.391706e-02, 1.009248e-02, 1.919228e-02, 3.590242e-02, 1.240861e-01, 0., 6.833079e-02, 3.805319e-02, 1.141599e-02, 1.803412e-02, 7.610617e-03} - {2.987659e-01, 4.370794e-02, 3.021117e-02, 1.931559e-02, 3.825913e-02, 1.400361e-01, 1.708689e-02, 2.129669e-02, 2.823048e-02, 5.398425e-02, 6.686193e-03, 2.151960e-01, 1.409036e-01, 2.253455e-02, 5.113617e-02, 0., 2.795782e-01, 2.711589e-02, 4.457448e-03, 2.030591e-02} - {3.522204e-01, 9.527538e-03, 2.179747e-02, 2.049851e-02, 5.629777e-03, 2.367402e-02, 1.097104e-02, 1.342509e-01, 5.745369e-02, 2.482901e-02, 4.951422e-02, 1.000388e-01, 5.947417e-02, 2.150877e-02, 3.320143e-02, 3.259548e-01, 0., 7.593050e-02, 1.154850e-03, 6.495937e-03} - {2.247538e-01, 1.259980e-02, 1.616311e-02, 2.876304e-02, 2.405365e-02, 3.512603e-02, 2.799949e-03, 5.011897e-01, 7.381676e-03, 1.604864e-01, 7.114419e-02, 6.999878e-03, 1.069059e-02, 7.381575e-03, 8.781546e-03, 2.787208e-02, 6.694340e-02, 0., 3.436295e-03, 5.345272e-03} - {6.481609e-03, 2.239168e-02, 2.946223e-03, 7.071031e-03, 2.180197e-02, 4.065785e-02, 2.946258e-03, 7.071054e-03, 5.303293e-03, 4.831833e-02, 4.714033e-03, 1.178523e-03, 3.535470e-03, 7.070948e-03, 6.422769e-02, 2.121298e-02, 4.713974e-03, 1.590964e-02, 0., 2.415905e-02} - {1.072149e-02, 4.288707e-02, 2.327365e-02, 3.922601e-03, 2.222760e-01, 3.922483e-03, 1.344137e-01, 1.595179e-02, 5.230110e-03, 2.196630e-02, 4.445615e-03, 2.981167e-02, 5.753066e-03, 1.046008e-02, 1.202907e-02, 4.288676e-02, 1.176761e-02, 1.098310e-02, 1.072173e-02, 0.} + +HIVWm={ +{0,0.167653,4.43521,5.56325,0.597923,1.8685,0.005,0.005,0.592784,0.16024,0.005,0.617509,1.00981,0.005,0.0744808,8.5942,24.1422,24.8094,0.005,0.005} +{0.167653,0,0.005,0.005,0.362959,0.0489798,0.005,0.005,0.005,0.005,0.005,0.0604932,0.005,0.005,2.86364,1.12195,0.005,0.005,5.49894,8.34835} +{4.43521,0.005,0,12.1233,0.005,10.3969,2.31779,0.145124,0.894313,0.005,0.005,29.4087,0.005,0.005,0.0674539,0.427881,0.630395,2.91786,0.005,2.28154} +{5.56325,0.005,12.1233,0,0.005,14.7801,0.005,0.0390512,23.9626,0.129839,0.005,0.201526,0.005,3.20656,0.0251632,0.005,0.458743,2.19952,0.005,0.005} +{0.597923,0.362959,0.005,0.005,0,0.005,0.005,1.48288,0.005,7.48781,0.005,0.005,0.0342252,0.005,0.005,4.27939,0.114512,2.28,0.005,4.12728} +{1.8685,0.0489798,10.3969,14.7801,0.005,0,0.005,0.005,0.279425,0.0489798,0.0489798,0.0604932,0.005,0.0604932,13.4379,6.27966,0.0489798,2.79622,2.8258,0.005} +{0.005,0.005,2.31779,0.005,0.005,0.005,0,0.005,0.22406,1.76382,0.005,8.59876,13.9444,18.5465,6.84405,0.725157,0.95956,0.827479,0.005,47.4889} +{0.005,0.005,0.145124,0.0390512,1.48288,0.005,0.005,0,0.817481,9.10246,17.3064,0.987028,0.005,0.0342252,1.34069,0.740091,9.36345,24.8231,0.005,0.114512} +{0.592784,0.005,0.894313,23.9626,0.005,0.279425,0.22406,0.817481,0,0.005,4.09564,10.6655,0.111928,13.0705,39.8897,0.005,4.04802,0.128065,0.005,0.005} +{0.16024,0.005,0.005,0.129839,7.48781,0.0489798,1.76382,9.10246,0.005,0,11.3839,0.005,9.83095,2.89048,0.586757,6.14396,0.005,2.95344,1.37031,0.005} +{0.005,0.005,0.005,0.005,0.005,0.0489798,0.005,17.3064,4.09564,11.3839,0,0.201526,0.005,0.005,3.28652,0.392575,7.41313,14.7683,0.005,0.579198} +{0.617509,0.0604932,29.4087,0.201526,0.005,0.0604932,8.59876,0.987028,10.6655,0.005,0.201526,0,0.344848,0.342068,0.16024,14.5699,4.54206,0.0744808,0.005,5.06475} +{1.00981,0.005,0.005,0.005,0.0342252,0.005,13.9444,0.005,0.111928,9.83095,0.005,0.344848,0,3.04502,0.404723,14.249,4.33701,0.005,0.005,0.005} +{0.005,0.005,0.005,3.20656,0.005,0.0604932,18.5465,0.0342252,13.0705,2.89048,0.005,0.342068,3.04502,0,10.6746,0.16024,0.203091,0.005,0.0443298,0.005} +{0.0744808,2.86364,0.0674539,0.0251632,0.005,13.4379,6.84405,1.34069,39.8897,0.586757,3.28652,0.16024,0.404723,10.6746,0,8.35024,0.928203,0.279425,5.96564,0.005} +{8.5942,1.12195,0.427881,0.005,4.27939,6.27966,0.725157,0.740091,0.005,6.14396,0.392575,14.5699,14.249,0.16024,8.35024,0,6.34079,0.862637,1.10156,0.933142} +{24.1422,0.005,0.630395,0.458743,0.114512,0.0489798,0.95956,9.36345,4.04802,0.005,7.41313,4.54206,4.33701,0.203091,0.928203,6.34079,0,0.005,0.005,0.490608} +{24.8094,0.005,2.91786,2.19952,2.28,2.79622,0.827479,24.8231,0.128065,2.95344,14.7683,0.0744808,0.005,0.005,0.279425,0.862637,0.005,0,0.005,1.35482} +{0.005,5.49894,0.005,0.005,0.005,2.8258,0.005,0.005,0.005,1.37031,0.005,0.005,0.005,0.0443298,5.96564,1.10156,0.005,0.005,0,0.005} +{0.005,8.34835,2.28154,0.005,4.12728,0.005,47.4889,0.114512,0.005,0.005,0.579198,5.06475,0.005,0.005,0.005,0.933142,0.490608,1.35482,0.005,0} }; -fprintf (stdout, mapMatrixToDict (models.protein.alphabet, models.protein.JTT.empirical_Q_normalized), "\n"); +cpREV={ +{0.0, 0.17842659999999999, 0.092912900000000007, 1.2410950000000001, 0.0089291809999999996, 0.19922690000000001, 0.95218210000000003, 1.8519509999999999, 5.2413160000000003, 0.1140412, 0.069691009999999998, 0.073883550000000006, 0.062992709999999993, 0.11301460000000001, 1.800713, 0.99883580000000005, 2.9123169999999998, 0.079395489999999999, 0.1433528, 3.7744770000000001} +{0.17842659999999999, 0.0, 0.78291299999999997, 0.057953739999999997, 0.1821885, 1.9239010000000001, 0.062738630000000004, 1.0893999999999999, 10.4955, 0.32451750000000001, 0.3932002, 7.5492400000000002, 0.33623259999999999, 0.082086770000000003, 0.3498713, 1.9264349999999999, 1.1352580000000001, 0.57242859999999995, 0.17113149999999999, 0.1366145} +{0.092912900000000007, 0.78291299999999997, 0.0, 7.1851820000000002, 1.3742680000000001e-06, 0.087059890000000001, 0.50383730000000004, 0.48689009999999999, 14.05444, 1.762721, 0.027694420000000001, 6.1903180000000004, 0.039721729999999997, 0.19554460000000001, 0.0073425540000000003, 7.3483460000000003, 2.1471749999999998, 0.00073109369999999996, 2.622763, 0.049312059999999998} +{1.2410950000000001, 0.057953739999999997, 7.1851820000000002, 0.0, 0.023400190000000001, 0.18438560000000001, 7.4266189999999996, 2.1124000000000001, 11.26995, 0.039169990000000002, 0.030205019999999999, 0.066227720000000004, 0.033575769999999998, 0.1031734, 0.1509482, 0.58229880000000001, 0.15168809999999999, 0.01423897, 0.90783380000000002, 0.40760740000000001} +{0.0089291809999999996, 0.1821885, 1.3742680000000001e-06, 0.023400190000000001, 0.0, 1.046446e-08, 7.5192150000000001e-11, 0.058911230000000002, 3.9633880000000001, 0.0006594967, 0.0060792190000000003, 3.7228779999999998e-16, 0.0072131779999999998, 0.1993818, 0.0048783949999999998, 0.26394820000000002, 3.225214e-06, 0.44408330000000001, 0.77416119999999999, 0.022435119999999999} +{0.19922690000000001, 1.9239010000000001, 0.087059890000000001, 0.18438560000000001, 1.046446e-08, 0.0, 3.6916709999999999, 0.055163400000000001, 8.9084339999999997, 6.7127360000000003e-06, 0.6802781, 3.030805, 0.001233336, 0.00149661, 0.74269090000000004, 0.00059064049999999998, 0.12020939999999999, 4.3329830000000002e-05, 0.027370909999999998, 0.0090477370000000001} +{0.95218210000000003, 0.062738630000000004, 0.50383730000000004, 7.4266189999999996, 7.5192150000000001e-11, 3.6916709999999999, 0.0, 1.38937, 7.2980799999999997, 0.00010299589999999999, 0.0012831209999999999, 3.608816, 0.076595659999999996, 0.052886250000000003, 0.028898150000000001, 0.067767090000000002, 0.060166240000000003, 0.02252612, 0.1240642, 0.57954090000000003} +{1.8519509999999999, 1.0893999999999999, 0.48689009999999999, 2.1124000000000001, 0.058911230000000002, 0.055163400000000001, 1.38937, 0.0, 9.1395180000000007, 0.035604820000000002, 0.021579359999999999, 0.055044000000000003, 0.021872639999999999, 0.19847719999999999, 0.079150559999999995, 0.99842149999999996, 0.078627669999999997, 0.13868530000000001, 0.22958419999999999, 0.42281999999999997} +{5.2413160000000003, 10.4955, 14.05444, 11.26995, 3.9633880000000001, 8.9084339999999997, 7.2980799999999997, 9.1395180000000007, 0.0, 4.7065859999999997, 5.8791029999999997, 1.455741, 2.298295, 5.642309, 10.494960000000001, 5.4391160000000003, 3.4432849999999999, 7.01389, 20.55414, 6.890244} +{0.1140412, 0.32451750000000001, 1.762721, 0.039169990000000002, 0.0006594967, 6.7127360000000003e-06, 0.00010299589999999999, 0.035604820000000002, 4.7065859999999997, 0.0, 1.6011230000000001, 0.50597930000000002, 10.96748, 2.7147049999999999, 0.050165679999999997, 0.60076070000000004, 3.0871520000000001, 0.063187480000000004, 0.29031649999999998, 7.9266750000000004} +{0.069691009999999998, 0.3932002, 0.027694420000000001, 0.030205019999999999, 0.0060792190000000003, 0.6802781, 0.0012831209999999999, 0.021579359999999999, 5.8791029999999997, 1.6011230000000001, 0.0, 0.021584510000000001, 5.6479850000000003, 3.3906179999999999, 1.149931, 0.1580539, 0.57027919999999999, 0.3378544, 0.15213199999999999, 3.59531} +{0.073883550000000006, 7.5492400000000002, 6.1903180000000004, 0.066227720000000004, 3.7228779999999998e-16, 3.030805, 3.608816, 0.055044000000000003, 1.455741, 0.50597930000000002, 0.021584510000000001, 0.0, 1.238634, 0.0046490350000000001, 0.0099489939999999992, 0.086884050000000004, 1.0392980000000001, 0.0080242630000000002, 0.07109973, 0.034934399999999997} +{0.062992709999999993, 0.33623259999999999, 0.039721729999999997, 0.033575769999999998, 0.0072131779999999998, 0.001233336, 0.076595659999999996, 0.021872639999999999, 2.298295, 10.96748, 5.6479850000000003, 1.238634, 0.0, 3.94794, 0.074172790000000002, 0.018613540000000001, 1.4156120000000001, 0.10111489999999999, 0.0022467590000000001, 4.3967200000000002} +{0.11301460000000001, 0.082086770000000003, 0.19554460000000001, 0.1031734, 0.1993818, 0.00149661, 0.052886250000000003, 0.19847719999999999, 5.642309, 2.7147049999999999, 3.3906179999999999, 0.0046490350000000001, 3.94794, 0.0, 0.35561979999999999, 0.98130640000000002, 0.036744859999999997, 0.2199856, 7.0744639999999999, 1.6439459999999999} +{1.800713, 0.3498713, 0.0073425540000000003, 0.1509482, 0.0048783949999999998, 0.74269090000000004, 0.028898150000000001, 0.079150559999999995, 10.494960000000001, 0.050165679999999997, 1.149931, 0.0099489939999999992, 0.074172790000000002, 0.35561979999999999, 0.0, 1.284651, 0.90571120000000005, 0.0055160740000000002, 0.19921330000000001, 0.2217442} +{0.99883580000000005, 1.9264349999999999, 7.3483460000000003, 0.58229880000000001, 0.26394820000000002, 0.00059064049999999998, 0.067767090000000002, 0.99842149999999996, 5.4391160000000003, 0.60076070000000004, 0.1580539, 0.086884050000000004, 0.018613540000000001, 0.98130640000000002, 1.284651, 0.0, 3.0585749999999998, 0.1385142, 0.8104751, 0.074770409999999995} +{2.9123169999999998, 1.1352580000000001, 2.1471749999999998, 0.15168809999999999, 3.225214e-06, 0.12020939999999999, 0.060166240000000003, 0.078627669999999997, 3.4432849999999999, 3.0871520000000001, 0.57027919999999999, 1.0392980000000001, 1.4156120000000001, 0.036744859999999997, 0.90571120000000005, 3.0585749999999998, 0.0, 0.01412361, 0.099842550000000002, 0.2166054} +{0.079395489999999999, 0.57242859999999995, 0.00073109369999999996, 0.01423897, 0.44408330000000001, 4.3329830000000002e-05, 0.02252612, 0.13868530000000001, 7.01389, 0.063187480000000004, 0.3378544, 0.0080242630000000002, 0.10111489999999999, 0.2199856, 0.0055160740000000002, 0.1385142, 0.01412361, 0.0, 0.61212840000000002, 0.096635689999999996} +{0.1433528, 0.17113149999999999, 2.622763, 0.90783380000000002, 0.77416119999999999, 0.027370909999999998, 0.1240642, 0.22958419999999999, 20.55414, 0.29031649999999998, 0.15213199999999999, 0.07109973, 0.0022467590000000001, 7.0744639999999999, 0.19921330000000001, 0.8104751, 0.099842550000000002, 0.61212840000000002, 0.0, 0.5010635} +{3.7744770000000001, 0.1366145, 0.049312059999999998, 0.40760740000000001, 0.022435119999999999, 0.0090477370000000001, 0.57954090000000003, 0.42281999999999997, 6.890244, 7.9266750000000004, 3.59531, 0.034934399999999997, 4.3967200000000002, 1.6439459999999999, 0.2217442, 0.074770409999999995, 0.2166054, 0.096635689999999996, 0.5010635, 0.0}}; + + + + + + +fprintf (stdout, mapMatrixToDict (models.protein.alphabet, cpREV), "\n"): diff --git a/res/TemplateBatchFiles/libv3/tasks/estimators.bf b/res/TemplateBatchFiles/libv3/tasks/estimators.bf index 5611e5fe1..05842380b 100644 --- a/res/TemplateBatchFiles/libv3/tasks/estimators.bf +++ b/res/TemplateBatchFiles/libv3/tasks/estimators.bf @@ -371,7 +371,7 @@ lfunction estimators.TraverseLocalParameters (likelihood_function_id, model_desc /** - * @name estimators.ApplyExistingEstimates + * @name * @param {String} likelihood_function_id * @param {Dictionary} model_descriptions * @param {Matrix} initial_values @@ -614,7 +614,7 @@ lfunction estimators.CreateLFObject (context, data_filter, tree, model_template, for (i = 0; i < components; i += 1) { lf_components[2 * i] = filters[i]; DataSetFilter ^ (filters[i]) = CreateFilter( ^ (data_filter[i]), 1); - } + } user_model_id = context + ".user_model"; utility.ExecuteInGlobalNamespace ("`user_model_id` = 0"); diff --git a/src/core/batchlan.cpp b/src/core/batchlan.cpp index 3bc253af5..e0963c2e6 100644 --- a/src/core/batchlan.cpp +++ b/src/core/batchlan.cpp @@ -901,6 +901,11 @@ void KillLFRecord (long lfID, bool completeKill) { } } + // unregister event listeners to release reference counts + + me->UnregisterListeners(); + + if (lfID* nodie, long number , _String nodeNam bool setDef = true; if (autoSolveBranchLengths) { long nodeModelID = cNt.GetModelIndex(); - if (nodeModelID != HY_NO_MODEL) { + if (nodeModelID != HY_NO_MODEL && !IsModelOfExplicitForm(nodeModelID)) { _Formula * expressionToSolveFor = nil; long alreadyConverted = convertedMatrixExpressions.Find ((BaseRef)nodeModelID); if (alreadyConverted < 0) { diff --git a/src/core/global_object_lists.cpp b/src/core/global_object_lists.cpp index 0edf02080..79c130f66 100644 --- a/src/core/global_object_lists.cpp +++ b/src/core/global_object_lists.cpp @@ -260,7 +260,6 @@ namespace hyphy_global_objects { _data_filter_listeners.Insert ((BaseRef)index, (long)current_listeners, false, false); } - if (current_listeners->_SimpleList::Find((long)listener) < 0L) { (*current_listeners) << listener; } diff --git a/src/core/include/likefunc.h b/src/core/include/likefunc.h index ca277115b..0a02eae6c 100644 --- a/src/core/include/likefunc.h +++ b/src/core/include/likefunc.h @@ -314,6 +314,7 @@ class _LikelihoodFunction: public BaseObj */ _AssociativeList*CollectLFAttributes (void) const; + void UnregisterListeners (void); protected: @@ -369,7 +370,7 @@ class _LikelihoodFunction: public BaseObj void LocateTheBump (long,_Parameter , _Parameter& , _Parameter&, _Parameter = -1.); void GradientLocateTheBump (_Parameter, _Parameter&, _Matrix&, _Matrix&); void GradientDescent (_Parameter& , _Matrix& ); - void ConjugateGradientDescent + _Parameter ConjugateGradientDescent (_Parameter , _Matrix& , bool localOnly = false, long = 0x7fffffff,_SimpleList* only_these_parameters = nil, _Parameter check_lf = A_LARGE_NUMBER); _Parameter SetParametersAndCompute diff --git a/src/core/likefunc.cpp b/src/core/likefunc.cpp index f806778a0..7e2f32f2b 100644 --- a/src/core/likefunc.cpp +++ b/src/core/likefunc.cpp @@ -56,7 +56,7 @@ using namespace hyphy_global_objects; - //#define _UBER_VERBOSE_LF_DEBUG 1 +//#define _UBER_VERBOSE_LF_DEBUG 1 //#define _COMPARATIVE_LF_DEBUG_DUMP //#define _COMPARATIVE_LF_DEBUG_CHECK @@ -753,6 +753,15 @@ void _LikelihoodFunction::Rebuild (bool rescan_parameters) { //_______________________________________________________________________________________ +void _LikelihoodFunction::UnregisterListeners (void) { + unsigned long partition_count = CountObjects(kLFCountPartitions); + for (unsigned long i = 0UL; i < partition_count; i++) { + UnregisterChangeListenerForDataFilter(theDataFilters.GetElement(i), this); + } +} + +//_______________________________________________________________________________________ + void _LikelihoodFunction::Clear (void) { DeleteCaches (); @@ -760,10 +769,7 @@ void _LikelihoodFunction::Clear (void) unsigned long partition_count = CountObjects(kLFCountPartitions); theTrees.Clear(); - - for (unsigned long i = 0UL; i < partition_count; i++) { - UnregisterChangeListenerForDataFilter(theDataFilters.GetElement(i), this); - } + UnregisterListeners (); theDataFilters.Clear(); theProbabilities.Clear(); @@ -1959,8 +1965,8 @@ _Parameter _LikelihoodFunction::Compute (void) _Variable *v = LocateVar (indexInd.lData[i]); if (v->HasChanged()) { fprintf (stderr, "[CHANGED] "); + fprintf (stderr, "%s = %15.12g\n", v->GetName()->sData, v->theValue); } - fprintf (stderr, "%s = %15.12g\n", v->GetName()->sData, v->theValue); } #endif if (computeMode == 0 || computeMode == 3) { @@ -4234,10 +4240,10 @@ DecideOnDivideBy (this); if (gradientBlocks.lLength) { for (long b = 0; b < gradientBlocks.lLength; b++) { - ConjugateGradientDescent (prec, bestMSoFar,true,10,(_SimpleList*)(gradientBlocks(b)),maxSoFar); + maxSoFar = ConjugateGradientDescent (prec, bestMSoFar,true,10,(_SimpleList*)(gradientBlocks(b)),maxSoFar); } } else { - ConjugateGradientDescent (prec, bestMSoFar,true,10,nil,maxSoFar); + maxSoFar = ConjugateGradientDescent (prec, bestMSoFar,true,10,nil,maxSoFar); } GetAllIndependent (bestMSoFar); @@ -4539,7 +4545,7 @@ DecideOnDivideBy (this); if (!skipCG && loopCounter&& indexInd.lLength>1 && ( (((long)loopCounter)%indexInd.lLength)==0 )) { _Matrix bestMSoFar; GetAllIndependent (bestMSoFar); - ConjugateGradientDescent (currentPrecision, bestMSoFar); + maxSoFar = ConjugateGradientDescent (currentPrecision, bestMSoFar); logLHistory.Store(maxSoFar); } } @@ -4556,7 +4562,7 @@ DecideOnDivideBy (this); if (optMethod == 7) { _Matrix bestMSoFar (indexInd.lLength,1,false,true); GetAllIndependent(bestMSoFar); - ConjugateGradientDescent (currentPrecision*.01, bestMSoFar); + maxSoFar = ConjugateGradientDescent (currentPrecision*.01, bestMSoFar); } DeleteObject (stepHistory); @@ -5013,16 +5019,16 @@ long _LikelihoodFunction::Bracket (long index, _Parameter& left, _Parameter& if (curVar) { if (CheckAndSetIthIndependent(index,middle)) { - CheckAndSetIthIndependent(index,left); + /*CheckAndSetIthIndependent(index,left); _Parameter lc = Compute(); CheckAndSetIthIndependent(index,right); _Parameter rc = Compute(); - CheckAndSetIthIndependent(index,middle); + CheckAndSetIthIndependent(index,middle);*/ middleValue = Compute(); if (verbosityLevel > 100) { char buf [256]; - snprintf (buf, 256, "\n\t[_LikelihoodFunction::Bracket (index %ld) recomputed the value to midpoint: L(%g) = %g [%g/%g:%g%g]]", index, middle, middleValue, leftValue,lc, rightValue,rc); + snprintf (buf, 256, "\n\t[_LikelihoodFunction::Bracket (index %ld) recomputed the value to midpoint: L(%g) = %g [@%g -> %g:@%g -> %g]]", index, middle, middleValue, left, leftValue,right, rightValue); BufferToConsole (buf); //exit (0); } @@ -5741,7 +5747,7 @@ bool _LikelihoodFunction::SniffAround (_Matrix& values, _Parameter& bestSoFar //_______________________________________________________________________________________ -void _LikelihoodFunction::ConjugateGradientDescent (_Parameter precision, _Matrix& bestVal, bool localOnly, long iterationLimit, _SimpleList* only_these_parameters, _Parameter check_value) { +_Parameter _LikelihoodFunction::ConjugateGradientDescent (_Parameter precision, _Matrix& bestVal, bool localOnly, long iterationLimit, _SimpleList* only_these_parameters, _Parameter check_value) { _Parameter gradientStep = STD_GRAD_STEP, temp, @@ -5758,7 +5764,7 @@ void _LikelihoodFunction::ConjugateGradientDescent (_Parameter precision, _Ma ReportWarning (_String ((_String*)optimizatonHistory->toStr())); } WarnError (errorStr); - return; + return check_value; } //return; } @@ -5863,7 +5869,7 @@ void _LikelihoodFunction::ConjugateGradientDescent (_Parameter precision, _Ma } if (terminateExecution) { - return; + return check_value; } } } @@ -5872,11 +5878,14 @@ void _LikelihoodFunction::ConjugateGradientDescent (_Parameter precision, _Ma if (maxSoFar < initial_value && CheckEqual(maxSoFar, initial_value) == false) { WarnError (_String("Internal optimization error in _LikelihoodFunction::ConjugateGradientDescent. Worsened likelihood score from ") & initial_value & " to " & maxSoFar); } + if (vl>1) { BufferToConsole("\n"); } - + + return maxSoFar; + } //_______________________________________________________________________________________ @@ -7767,7 +7776,7 @@ _Parameter _LikelihoodFunction::ComputeBlock (long index, _Parameter* siteRes, (_GrowingVector*)conditionalTerminalNodeLikelihoodCaches(index), overallScalingFactors.lData[index], blockID * sitesPerP, - (1+blockID) * sitesPerP, + (1L+blockID) * sitesPerP, catID, siteRes, scc, @@ -7775,7 +7784,8 @@ _Parameter _LikelihoodFunction::ComputeBlock (long index, _Parameter* siteRes, branchIndex >= 0 ? branchValues->lData: nil); } - + + if (np > 1) { _Parameter correction = 0.; for (blockID = 0; blockID < np; blockID ++) { @@ -7785,12 +7795,34 @@ _Parameter _LikelihoodFunction::ComputeBlock (long index, _Parameter* siteRes, sum = temp_sum; } - } else { sum = thread_results[0]; - } + /*#ifdef _UBER_VERBOSE_LF_DEBUG + static _Parameter previous_results [4096]; + static long previous_scalers [4096]; + + if (likeFuncEvalCallCount > 12050) { + abort (); + } + + if (likeFuncEvalCallCount > 12000) { + fprintf (stderr, "Block sum %g\n", sum); + for (long i = 0L; i < sitesPerP - 1; i++) { + fprintf (stderr, "[%ld] site %ld \t %g (%g) [delta %g]; scaler %ld (%ld)\n", catID, i, siteRes[i], previous_results [catID * sitesPerP + i], fabs(siteRes[i]-previous_results [catID * sitesPerP + i]), scc [i], previous_scalers [catID * sitesPerP + i]); + } + + } + if (likeFuncEvalCallCount >= 12000) { + + for (long i = 0L; i < sitesPerP - 1; i++) { + previous_results [catID * sitesPerP + i] = siteRes[i]; + previous_scalers [catID * sitesPerP + i] = scc [i]; + } + } + #endif*/ + delete [] thread_results; /* #pragma omp parallel for default(shared) schedule(static,1) private(blockID) num_threads (np) reduction(+:sum) if (np>1) diff --git a/src/core/likefunc2.cpp b/src/core/likefunc2.cpp index 01c65b955..5b78b241f 100644 --- a/src/core/likefunc2.cpp +++ b/src/core/likefunc2.cpp @@ -719,18 +719,23 @@ void _LikelihoodFunction::PopulateConditionalProbabilities (long in for (long r1 = lowerBound, r2 = lowerBound2; r1 < upperBound; r1++,r2++) { if (siteCorrectors) { long scv = *siteCorrectors; - - if (scv < scalers.lData[r1]) { // this class has a _smaller_ scaling factor - buffer[r1] = currentRateWeight * buffer[r2] + buffer[r1] * acquireScalerMultiplier (scalers.lData[r1] - scv); + + if (currentRateCombo == 0L) { // first entry + buffer[r1] = currentRateWeight * buffer[r2]; scalers.lData[r1] = scv; } else { - if (scv > scalers.lData[r1]) { // this is a _larger_ scaling factor - buffer[r1] += currentRateWeight * buffer[r2] * acquireScalerMultiplier (scv - scalers.lData[r1]); - } else { // same scaling factors - buffer[r1] += currentRateWeight * buffer[r2]; + if (scv < scalers.lData[r1]) { // this class has a _smaller_ scaling factor + buffer[r1] = currentRateWeight * buffer[r2] + buffer[r1] * acquireScalerMultiplier (scalers.lData[r1] - scv); + scalers.lData[r1] = scv; + } else { + if (scv > scalers.lData[r1]) { // this is a _larger_ scaling factor + buffer[r1] += currentRateWeight * buffer[r2] * acquireScalerMultiplier (scv - scalers.lData[r1]); + } else { // same scaling factors + buffer[r1] += currentRateWeight * buffer[r2]; + } } } - + siteCorrectors++; } else { buffer[r1] += currentRateWeight * buffer[r2]; @@ -1231,7 +1236,8 @@ void _LikelihoodFunction::CleanupParameterMapping (void) parameterTransformationFunction.Clear(); } - +//#define _UBER_VERBOSE_LF_DEBUG 1 +//extern long likeFuncEvalCallCount; //_______________________________________________________________________________________________ @@ -1270,20 +1276,30 @@ _Parameter _LikelihoodFunction::SumUpSiteLikelihoods (long index, const _Paramet WarnError ("Constant-on-partition categories are currently not supported by the evaluation engine"); } else { for (unsigned long patternID = 0UL; patternID < pattern_count; patternID++) { - long patternFrequency = index_filter->GetFrequency(patternID);; + + + long patternFrequency = index_filter->GetFrequency(patternID); if (patternFrequency > 1) { logL += myLog(patternLikelihoods[patternID])*patternFrequency; cumulativeScaler += patternScalers.lData[patternID]*patternFrequency; - } else - // all this to avoid a double*long multiplication - { + } else { + // avoid a double*long multiplication logL += myLog(patternLikelihoods[patternID]); cumulativeScaler += patternScalers.lData[patternID]; } +#ifdef _UBER_VERBOSE_LF_DEBUG + if (likeFuncEvalCallCount > 12000) { + fprintf (stderr, "[%ld] %g (scaler %ld)\n", patternID, logL, cumulativeScaler); + } +#endif } } } - +#ifdef _UBER_VERBOSE_LF_DEBUG + if (likeFuncEvalCallCount > 12050) { + abort(); + } +#endif return logL - cumulativeScaler * _logLFScaler; } diff --git a/src/core/site.cpp b/src/core/site.cpp index 8beacac09..638e6926e 100644 --- a/src/core/site.cpp +++ b/src/core/site.cpp @@ -346,9 +346,9 @@ long _TranslationTable::MultiTokenResolutions (_String const& tokens, long* r *resolution_arrays [3] = {large_store + tokens.sLength, large_store + tokens.sLength + baseLength,large_store + tokens.sLength + 2*baseLength}, resolutions_index = 0L; - for (unsigned long digit1 = 0; digit1 < large_store[0]; digit1 ++) { - for (unsigned long digit2 = 0; digit2 < large_store[1]; digit2 ++) { - for (unsigned long digit3 = 0; digit3 < large_store[2]; digit3 ++) { + for (long digit1 = 0L; digit1 < large_store[0]; digit1 ++) { + for (long digit2 = 0L; digit2 < large_store[1]; digit2 ++) { + for (long digit3 = 0L; digit3 < large_store[2]; digit3 ++) { receptacle[resolutions_index++] = resolution_arrays[0][digit1] * baseLength * baseLength + resolution_arrays[1][digit2] * baseLength + resolution_arrays[2][digit3]; } } @@ -360,8 +360,8 @@ long _TranslationTable::MultiTokenResolutions (_String const& tokens, long* r *resolution_arrays [2] = {large_store + tokens.sLength,large_store + tokens.sLength + baseLength}, resolutions_index = 0L; - for (unsigned long digit1 = 0; digit1 < large_store[0]; digit1 ++) { - for (unsigned long digit2 = 0; digit2 < large_store[1]; digit2 ++) { + for (long digit1 = 0L; digit1 < large_store[0]; digit1 ++) { + for (long digit2 = 0L; digit2 < large_store[1]; digit2 ++) { receptacle[resolutions_index++] = resolution_arrays[0][digit1] * baseLength + resolution_arrays[1][digit2]; } } @@ -2426,7 +2426,7 @@ void _DataSetFilter::SetFilter (_DataSet const * ds, unsigned char unit, _Sim // security checks if (horizontalList.empty() || verticalList.lLengththeNodeMap.lLength : ds->NoOfSpecies(),0,1); } @@ -2865,7 +2865,7 @@ _String* _DataSetFilter::GetExclusions (void) const { //_______________________________________________________________________ unsigned long _DataSetFilter::GetDimension (bool correct) const { - unsigned long result = ComputePower (theData->theTT->baseLength, unitLength); + unsigned long result = ComputePower (theData->theTT->LengthOfAlphabet(), unitLength); return correct ? result - theExclusions.lLength : result; } @@ -6114,7 +6114,7 @@ _Matrix * _DataSet::HarvestFrequencies (unsigned char unit, unsigned char atom, } - _Matrix * out = new _Matrix (ComputePower (theTT->baseLength, atom), + _Matrix * out = new _Matrix (ComputePower (theTT->LengthOfAlphabet(), atom), posSpec?unit/atom:1, false, true); diff --git a/src/mains/unix.cpp b/src/mains/unix.cpp index b47b07ab1..d9a27a397 100644 --- a/src/mains/unix.cpp +++ b/src/mains/unix.cpp @@ -46,6 +46,7 @@ SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. const char hy_usage[] = "usage: HYPHYMP or HYPHYMPI [-h] " +"[-v] " "[-c] " "[-d] " "[-p] " @@ -61,6 +62,7 @@ const char hy_help_message [] = "\n" "optional flags:\n" " -h show this help message and exit\n" +" -v print HyPhy version\n" " -c calculator mode; causes HyPhy to drop into an expression evaluation until 'exit' is typed\n" " -d debug mode; causes HyPhy to drop into an expression evaluation mode upon script error\n" " -p postprocessor mode; drops HyPhy into an interactive mode where general post-processing scripts can be selected\n" @@ -77,6 +79,7 @@ const char hy_help_message [] = " analysis arguments if batch file is present, all remaining positional arguments are interpreted as inputs to analysis prompts\n" ; + #ifdef _MINGW32_MEGA_ #include HANDLE _HY_MEGA_Pipe = INVALID_HANDLE_VALUE; @@ -525,6 +528,12 @@ void ProcessConfigStr (_String& conf) exit (0); } + case 'v': + case 'V': { + fprintf(stderr, "\n%s\n", GetVersionString().getStr() ); + exit (0); + } + case 'p': case 'P': { usePostProcessors = true; diff --git a/tests/hbltests/libv3/LEISR.wbf b/tests/hbltests/libv3/LEISR.wbf index 9d4bbc8d9..b01d3c10a 100644 --- a/tests/hbltests/libv3/LEISR.wbf +++ b/tests/hbltests/libv3/LEISR.wbf @@ -4,46 +4,43 @@ utility.ToggleEnvVariable ("OPTIMIZATION_TIME_HARD_LIMIT", 1); LoadFunctionLibrary("LEISR.bf", { - "0": "Protein", - "1": "WAG", - "2": "Gamma", - "3": "Yes", - "4": PATH_TO_CURRENT_BF + "data/dat.prot", - "5": "Y" + "0": PATH_TO_CURRENT_BF + "data/dat.prot", + "1": "Y", + "2": "Protein", + "3": "WAG", + "4": "Gamma" }); LoadFunctionLibrary("LEISR.bf", { - "0": "Protein", - "1": "WAG", - "2": "GDD", - "3": "No", - "4": PATH_TO_CURRENT_BF + "data/dat.prot", - "5": "Y" + "0": PATH_TO_CURRENT_BF + "data/dat.prot", + "1": "Y", + "2": "Protein", + "3": "WAG", + "4": "GDD" }); LoadFunctionLibrary("LEISR.bf", { - "0": "Protein", - "1": "WAG", - "2": "No", - "3": "No", - "4": PATH_TO_CURRENT_BF + "data/dat.prot", - "5": "Y" + "0": PATH_TO_CURRENT_BF + "data/dat.prot", + "1": "Y", + "2": "Protein", + "3": "WAG", + "4": "No" }); LoadFunctionLibrary("LEISR.bf", { - "0": "Nucleotide", - "1": "GTR", - "2": "Gamma", - "3": PATH_TO_CURRENT_BF + "data/dat.nuc", - "4": "Y" + "0": PATH_TO_CURRENT_BF + "data/dat.nuc", + "1": "Y", + "2": "Nucleotide", + "3": "GTR", + "4": "Gamma" }); LoadFunctionLibrary("LEISR.bf", { - "0": "Nucleotide", - "1": "GTR", - "2": "No", - "3": PATH_TO_CURRENT_BF + "data/dat.nuc", - "4": "Y" + "0": PATH_TO_CURRENT_BF + "data/dat.nuc", + "1": "Y", + "2": "Nucleotide", + "3": "HKY85", + "4": "No" }); \ No newline at end of file diff --git a/tests/hbltests/libv3/WAG_MLE-freqs.bf b/tests/hbltests/libv3/WAG_MLE-freqs.bf new file mode 100644 index 000000000..4fb449218 --- /dev/null +++ b/tests/hbltests/libv3/WAG_MLE-freqs.bf @@ -0,0 +1,35 @@ +LoadFunctionLibrary("libv3/UtilityFunctions.bf"); +LoadFunctionLibrary("libv3/IOFunctions.bf"); +LoadFunctionLibrary("libv3/models/protein/empirical.bf"); +LoadFunctionLibrary("libv3/tasks/alignments.bf"); +LoadFunctionLibrary("libv3/tasks/trees.bf"); +LoadFunctionLibrary("libv3/tasks/estimators.bf"); + +file_name = "`PATH_TO_CURRENT_BF`data/CD2.prot"; +wag_mlf.data = {}; +wag_mlf.filter = {}; + +wag_mlf.alignment_info = alignments.ReadNucleotideAlignment(file_name, "wag_mlf.data", "wag_mlf.filter"); + +wag_mlf.partitions_and_trees = trees.LoadAnnotatedTreeTopology.match_partitions(wag_mlf.alignment_info[terms.data.partitions], {}); + +wag_mlf.partitions = alignments.DefineFiltersForPartitions (wag_mlf.partitions_and_trees, "wag_mlf.data" , "wag_mlf.filter.", wag_mlf.alignment_info); + +trees = utility.Map (wag_mlf.partitions_and_trees, "_partition_", '_partition_[terms.data.tree]'); +filter_names = utility.Map (wag_mlf.partitions, "_partition_", '_partition_[terms.data.name]'); + + + +wag_freq.results = estimators.FitSingleModel_Ext(filter_names, + trees, + "models.protein.WAGF.ModelDescription", + None,None); + +console.log ("WAG + F, Log(L) => " + wag_freq.results[terms.fit.log_likelihood]); + +wag_mlf.results = estimators.FitSingleModel_Ext(filter_names, + trees, + "models.protein.WAGML.ModelDescription", + None,None); + +console.log ("WAG + ML, Log(L) => " + wag_mlf.results[terms.fit.log_likelihood]); diff --git a/tests/hbltests/libv3/data/dat.nuc b/tests/hbltests/libv3/data/dat.nuc index 516ec504c..8a70180f9 100644 --- a/tests/hbltests/libv3/data/dat.nuc +++ b/tests/hbltests/libv3/data/dat.nuc @@ -1,31 +1,31 @@ >t8 -GACAAGTACC +GACAAGTACCCCCC >t9 -AATAGGTTAA +AATAGGTTAATTTT >t6 -AAGCGTAGAC +AAGCGTAGACCCCC >t7 -GGCCATCACA +GGCCATCACATTTT >t4 -GGGCAGTAAT +GGGCAGTAATCCCC >t5 -ACCCTGGCGT +ACCCTGGCGTTTTT >t2 -GCTCCGTACA +GCTCCGTACAAAAA >t3 -AGCATTACGA +AGCATTACGATTTT >t1 -GGGCCGTGAT +GGGCCGTGATTTTT >t14 -GGCGGTGTTC +GGCGGTGTTCAAAA >t15 -TCGATTCCGT +TCGATTCCGTAAAA >t10 -GCTCCGCATA +GCTCCGCATAAAAA >t11 -GGCAATATTA +GGCAATATTACCCC >t12 -AGAGATCCTT +AGAGATCCTTCCCC >t13 -GGCCATCACA +GGCCATCACACCCC (((t12:0.6007418204,(((t13:0.02560783178,t7:0.06241315883):0.3850171629,((t1:0.3803801907,t4:0.2644045493):0.288204158,((t8:0.2050027882,(t10:0.2882345759,t2:0.2495629364):0.4606286683):0.5220576057,t5:0.8048256312):0.05435458641):0.2048651455):0.4364210158,(t6:0.5203341169,t9:0.4957795972):0.9166075331):0.5339687597):0.2889535469,(t15:0.8320546721,t14:0.7063525459):0.2664719955):0.3073638263,(t3:0.4359256458,t11:0.7989921882):0.04232083797);