/*_____________________________________________________________________ */ function computeScalingFactorB(rateMatrix, baseFreqs) { B = 0; for (n1 = 0; n1 < Rows(rateMatrix); n1 = n1+1) { for (n2 = 0; n2 < Columns(rateMatrix); n2 = n2+1) { if (n2 != n1) { B = B + baseFreqs[n1]*baseFreqs[n2]*rateMatrix[n1][n2]; } } } return B; } /*_____________________________________________________________________ */ NICETY_LEVEL = 3; VERBOSITY_LEVEL = 0; dummy = HYPHY_BASE_DIRECTORY + "TemplateBatchFiles" + DIRECTORY_SEPARATOR + "TemplateModels" + DIRECTORY_SEPARATOR + "chooseGeneticCode.def"; ExecuteCommands ("#include \""+dummy+"\";"); SetDialogPrompt ("Please specify a codon data file:"); DataSet codon_ds = ReadDataFile (PROMPT_FOR_FILE); DataSetFilter codon_dsf = CreateFilter (codon_ds,3,"","",GeneticCodeExclusions); fprintf (stdout,"\n______________READ THE FOLLOWING DATA______________\n", codon_ds); /*-------------------------------------------------------------------------------------------------------------*/ /* generate codon model (MG94customModel) */ ModelMatrixDimension = 0; global R = 1; global AC = 1; global AT = 1; global CG = 1; global CT = 1; global GT = 1; function PopulateModelMatrix (ModelMatrixName&, EFV) { if (!ModelMatrixDimension) { ModelMatrixDimension = 64; for (h = 0 ;h<64; h=h+1) { if (_Genetic_Code[h]==10) { ModelMatrixDimension = ModelMatrixDimension-1; } } } R = R; AT = AT; CG = CG; CT = CT; GT = GT; AC = AC; ModelMatrixName = {ModelMatrixDimension,ModelMatrixDimension}; hshift = 0; if (modelType < 2) { for (h=0; h<64; h=h+1) { if (_Genetic_Code[h]==10) { hshift = hshift+1; continue; } vshift = hshift; for (v = h+1; v<64; v=v+1) { diff = v-h; if (_Genetic_Code[v]==10) { vshift = vshift+1; continue; } nucPosInCodon = 2; if ((h$4==v$4)||((diff%4==0)&&(h$16==v$16))||(diff%16==0)) { if (h$4==v$4) { transition = v%4; transition2= h%4; } else { if(diff%16==0) { transition = v$16; transition2= h$16; nucPosInCodon = 0; } else { transition = v%16$4; transition2= h%16$4; nucPosInCodon = 1; } } if (transition 1) { categoriesUsed = 1; if (modelType == 2) { dummy = HYPHY_BASE_DIRECTORY + "TemplateBatchFiles" + DIRECTORY_SEPARATOR + "TemplateModels" + DIRECTORY_SEPARATOR + "defineGamma.mdl"; ExecuteCommands ("#include \""+dummy+"\";"); /* #include "defineGamma.mdl"; */ } if (modelType == 3) { dummy = HYPHY_BASE_DIRECTORY + "TemplateBatchFiles" + DIRECTORY_SEPARATOR + "TemplateModels" + DIRECTORY_SEPARATOR + "defineHM.mdl"; ExecuteCommands ("#include \""+dummy+"\";"); /* #include "defineHM.mdl"; */ } } if (!SKIP_MODEL_PARAMETER_LIST) { 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; } } } MGCustomRateBiasTerms = {{"AC","1","AT","CG","CT","GT"}}; paramCount = 0; MGCustomModelConstraintString = ""; for (customLoopCounter2=1; customLoopCounter2<6; customLoopCounter2=customLoopCounter2+1) { for (customLoopCounter=0; customLoopCounter