hyphy-common 2.2.7+dfsg-1 / usr / share / hyphy / TemplateBatchFiles / AddABias.ibf

pThresh = 0.001;

/*--------------------------------------------------------------------------------------------*/

function AddABias (ModelMatrixName&, ModelMatrixName2&, biasedBase)
{
	ModelMatrixName2 = {20,20};
	
	t = 1;
	c = 1;
	_numericRateMatrix = ModelMatrixName;
	
	for (ri = 0; ri < 20; ri = ri+1)
	{
		for (ci = ri+1; ci < 20; ci = ci+1)
		{
			ModelMatrixName2[ri][ci] := _numericRateMatrix__[ri__][ci__] * t;
			ModelMatrixName2[ci][ri] := _numericRateMatrix__[ri__][ci__] * t;
		
		}
	}

	if (biasedBase < 20)
	{
		global rateBiasTo 	 = 2;
		global rateBiasFrom	:= 1/rateBiasTo;
		
		rateBiasTo   :>1;
		for (ri = 0; ri < 20; ri = ri+1)
		{
			if (ri != biasedBase)
			{
				ModelMatrixName2[ri][biasedBase] := _numericRateMatrix__[ri__][biasedBase__] * t * rateBiasTo;
				ModelMatrixName2[biasedBase][ri] := _numericRateMatrix__[ri__][biasedBase__] * t * rateBiasFrom;
			}
		}
	}

	return 1;
}

/*--------------------------------------------------------------------------------------------*/

function AddABiasREL (ModelMatrixName&, ModelMatrixName2&, biasedBase)
{
	ModelMatrixName2 = {20,20};
	
	t = 1;	/* branch length, local parameter */
	c = 1;	/* rate variation */
	_numericRateMatrix = ModelMatrixName;
	
	/* the probability that site is undergoing biased substitution rates */
	global	  P_bias = 0.1;  P_bias :< 0.5;
	
	
	category catVar = (2,{{1-P_bias,P_bias}},MEAN,,{{0,1}},0,1);
	
	for (ri = 0; ri < 20; ri = ri+1)
	{
		for (ci = ri+1; ci < 20; ci = ci+1)
		{
			ModelMatrixName2[ri][ci] := _numericRateMatrix__[ri__][ci__] * t * c;
			ModelMatrixName2[ci][ri] := _numericRateMatrix__[ri__][ci__] * t * c;
		
		}
	}

	if (biasedBase < 20)
	{
		global rateBiasTo 	  = 5.0;
		global rateBiasFrom	 := 1/rateBiasTo;
			
		rateBiasTo    :>1;
		relBias       :>1;	/* UNUSED ?!? */
		for (ri = 0; ri < 20; ri = ri+1)
		{
			if (ri != biasedBase)
			{
				ModelMatrixName2[ri][biasedBase] := _numericRateMatrix__[ri__][biasedBase__] * t * c * ((catVar==1)*rateBiasTo+(catVar==0));
				ModelMatrixName2[biasedBase][ri] := _numericRateMatrix__[ri__][biasedBase__] * t * c * ((catVar==1)*rateBiasFrom+(catVar==0));
			}
		}
	}

	return 1;
}

/*--------------------------------------------------------------------------------------------*/

function runAFit (refClass,ancClass,mySiteID)
{
	singletonScale 			= 1;
	DataSetFilter			aSite = CreateFilter (ds,1,siteIndex==mySiteID,speciesIndex<ds.species-2);
	deltaDistr	   			= {20,1}["_MATRIX_ELEMENT_ROW_=="+ancClass];
	LikelihoodFunction3	    singletonLF = (aSite, singletonTree, deltaDistr);
	Optimize 				(resS, singletonLF);
	AddABias 				(modelNameString,"biasedMatrix", refClass);
	Model biasedModel 		= (biasedMatrix, vectorOfFrequencies, 1);
	Tree biasedTree 		= treeString;
	global biasedScale	 	= 1;
	ReplicateConstraint 	("this1.?.t:=biasedScale*this2.?.t__",biasedTree,givenTree);
	LikelihoodFunction3	    biasedLF = (aSite, biasedTree, deltaDistr);
	Optimize 				(resB, biasedLF);
	
	checkIfDone = {};
	checkIfDone ["ref"] 	  = AAString[refClass];
	checkIfDone ["anc"] 	  = AAString[ancClass];
	checkIfDone ["p"]		  = 1-CChi2(2(resB[1][0]-resS[1][0]),1);
	checkIfDone ["bias"]	  = rateBiasTo;
	checkIfDone ["LL1"]	  	  = resS[1][0];
	checkIfDone ["LL2"]	  	  = resB[1][0];
	checkIfDone ["s1"]	  	  = singletonScale;
	checkIfDone ["s2"]	  	  = biasedScale;
	
	return checkIfDone;
}

/*--------------------------------------------------------------------------------------------*/

function runAFit2 (ancClass,mySiteID)
{
	DataSetFilter			aSite = CreateFilter (ds,1,siteIndex==mySiteID,speciesIndex<ds.species-2);
	COUNT_GAPS_IN_FREQUENCIES     = 0;
	HarvestFrequencies			  (siteFreqs,aSite,1,1,1);
	singletonScale 			= 1;

	DataSetFilter			aSite = CreateFilter (ds,1,siteIndex==mySiteID,speciesIndex<ds.species-2);
	deltaDistr	   			= {20,1}["_MATRIX_ELEMENT_ROW_=="+ancClass];
	LikelihoodFunction3	    singletonLF = (aSite, singletonTree, deltaDistr);
	Optimize 				(resS, singletonLF);
	
	checkIfDone = {};
	checkIfDone ["anc"] 	  = AAString[ancClass];
	checkIfDone ["LL1"]	  	  = resS[1][0];
	checkIfDone ["s1"]	  	  = singletonScale;
	checkIfDone ["p"]		  = 1;

	minP					  = 1;

	for (_resID = 0; _resID < 20; _resID = _resID + 1)
	{
		if (siteFreqs[_resID])
		{
			AddABias 				(modelNameString,"biasedMatrix", _resID);
			Model biasedModel 		= (biasedMatrix, vectorOfFrequencies, 1);
			Tree biasedTree 		= treeString;
			global biasedScale	 	= 1;
			ReplicateConstraint 	("this1.?.t:=biasedScale*this2.?.t__",biasedTree,givenTree);
			LikelihoodFunction3	    biasedLF = (aSite, biasedTree, deltaDistr);
			Optimize 				(resB, biasedLF);	
			pv						  = Min(20*(1-CChi2(2(resB[1][0]-resS[1][0]),1)),1);
			if (pv < minP)
			{
				minP = pv;
				checkIfDone ["ref"] 	  = AAString[_resID];
				checkIfDone ["bias"]	  = rateBiasTo;
				checkIfDone ["LL2"]	  	  = resB[1][0];
				checkIfDone ["s2"]	  	  = biasedScale;
				checkIfDone ["p"]		  = minP;
			}
		}
	}
	

	return checkIfDone;
}

/*--------------------------------------------------------------------------------------------*/

function promptModel (dummy)
{
	ChoiceList	     (pickAModel,"Subsitution Model",1,SKIP_NONE, "HIV Within","HIV Within",
							 								      "HIV Between","HIV Between",
							 								      "JTT","JTT",
							 								      "Flu H5N1", "Empirical model for H5N1 Influenza",
																  "LG", "Le-Gasquel 2008",
							 								      "REV", "Use general time-reversible model (WARNING: 189 rate parameters will be estimated from your alignment)");
							 								     
							 								    
	if (pickAModel < 0)
	{
		return 0;
	}

	modelSTDINoverload = {};
	modelSTDINoverload["2"] = "Rate variation";
	modelSTDINoverload["3"] = "Beta-Gamma";
	modelSTDINoverload["4"] = "4";
	modelNameString = "_customAAModelMatrix";

	if (pickAModel == 0)
	{
		modelSTDINoverload["0"] = HYPHY_LIB_DIRECTORY + "TemplateBatchFiles" + DIRECTORY_SEPARATOR + "TemplateModels" + DIRECTORY_SEPARATOR + "EmpiricalAA" + 
								  DIRECTORY_SEPARATOR + "HIVWithin";
	}

	if (pickAModel == 1)
	{
		modelSTDINoverload["0"] = HYPHY_LIB_DIRECTORY + "TemplateBatchFiles" + DIRECTORY_SEPARATOR + "TemplateModels" + DIRECTORY_SEPARATOR + "EmpiricalAA" + 
								  DIRECTORY_SEPARATOR + "HIVBetween";
	}

	if (pickAModel == 2)
	{
		modelSTDINoverload["0"] = HYPHY_LIB_DIRECTORY + "TemplateBatchFiles" + DIRECTORY_SEPARATOR + "TemplateModels" + DIRECTORY_SEPARATOR + "EmpiricalAA" + 
								  DIRECTORY_SEPARATOR + "JTT";
	}

	if (pickAModel == 3)
	{
		modelSTDINoverload["0"] = HYPHY_LIB_DIRECTORY + "TemplateBatchFiles" + DIRECTORY_SEPARATOR + "TemplateModels" + DIRECTORY_SEPARATOR + "EmpiricalAA" + 
								  DIRECTORY_SEPARATOR + "H5N1";
	}
	
	if (pickAModel == 4)
	{
		modelSTDINoverload["0"] = HYPHY_LIB_DIRECTORY + "TemplateBatchFiles" + DIRECTORY_SEPARATOR + "TemplateModels" + DIRECTORY_SEPARATOR + "EmpiricalAA" + 
								  DIRECTORY_SEPARATOR + "LG";
	}

	if (pickAModel < 5)
	{
		/* estimate base frequencies as model parameters */
		modelSTDINoverload["1"]	= "ML Estimates";
		modelPath = HYPHY_LIB_DIRECTORY + "TemplateBatchFiles" + DIRECTORY_SEPARATOR + "TemplateModels" + DIRECTORY_SEPARATOR + "Custom_AA_empirical.mdl";
	}
	else
	{
		modelNameString = "mtREVMatrix";
		modelPath = HYPHY_LIB_DIRECTORY + "TemplateBatchFiles" + DIRECTORY_SEPARATOR + "TemplateModels" + DIRECTORY_SEPARATOR + "mtREV.mdl";	
	}
	
	ExecuteAFile (modelPath,modelSTDINoverload);

	if (dummy)
	{
		ChoiceList	     (pickATarget,"Target residue",1,SKIP_NONE, "Fixed","Fixed Sequence",
							 								    "Inferred","Inferred Sequence");
	}						 								    
	return pickAModel;
}

/*--------------------------------------------------------------------------------------------*/

function processAJob (mySiteID)
{
	if (MPI_NODE_COUNT > 1)
	{
		MPIReceive 					(-1,fromNode,mpiResults);
		mySiteID 					= MPINodeStatus[fromNode-1];
		ExecuteCommands 			(mpiResults);
		checkIfDone 			  = _hyphyAssociativeArray;
	}
		
	fprintf (stdout, "\nSite ", mySiteID+1, 
					 "\n\t Character : ", checkIfDone ["anc"],"->",checkIfDone ["ref"],
					 "\n\t Bias term : ", checkIfDone ["bias"],
					 "\n\t LogL std  : ", checkIfDone ["LL1"],
					 "\n\t LogL bias : ", checkIfDone ["LL2"],
					 "\n\t Scale std : ", checkIfDone ["s1"],
					 "\n\t Scale bias: ", checkIfDone ["s2"],
					 "\n\t p-value   : ", checkIfDone ["p"]);
					 
	fprintf (outPath, "\nSite ", mySiteID+1, 
					 "\n\t Character : ", checkIfDone ["anc"],"->",checkIfDone ["ref"],
					 "\n\t Bias term : ", checkIfDone ["bias"],
					 "\n\t LogL std  : ", checkIfDone ["LL1"],
					 "\n\t LogL bias : ", checkIfDone ["LL2"],
					 "\n\t Scale std : ", checkIfDone ["s1"],
					 "\n\t Scale bias: ", checkIfDone ["s2"],
					 "\n\t p-value   : ", checkIfDone ["p"]);

	if (checkIfDone ["p"] < pThresh)
	{
		mx2 = smx;
		strongEvidence[Abs(strongEvidence)] = mx2;
	}
	return fromNode-1;
}