Dear Albert,
There are two steps involved in computing expected syn and non-syn subs for each branch (in the batch language):
1). Based on the genetic code (which is stored in the
_Genetic_Code variable when you are using chooseGeneticCode.def to select a genetic code), one needs to compute a stencil for rate matrix elements: which entries correspond to syn (or non-syn) substitutions. You can check whether this variable is defined by checking if Abs(_Genetic_Code ) is positive) This step should only be done once and can be lifted straight out of 'postdNdS.bf':
Code:/* do once */
branchNames = BranchName (givenTree,-1);
T = Columns (branchNames);
ExecuteCommands(
"GetInformation (aRateMx, givenTree."+branchNames[0]+");");
/* make syn and non-syn template matrices */
nonStopCount = Columns (aRateMx);
vertOnes = {nonStopCount,1};
horOnes = {1,nonStopCount};
for (h1 = 0; h1<nonStopCount; h1=h1+1)
{
vertOnes [h1] = 1;
horOnes [h1] = 1;
}
/* do for each replicate, here simulatedCodonEFV is the vector of
equilibrium codon frequencies for the iterate which is computed to define the model for a replicate, e.g. in simpleBootstrap.bf */
synM = {nonStopCount,nonStopCount};
nonSynM = {nonStopCount,nonStopCount};
hShift = 0;
for (h1 = 0; h1 < 64; h1=h1+1)
{
gc1 = _Genetic_Code[h1];
if (gc1 == 10)
{
hShift = hShift+1;
}
else
{
vShift = hShift;
for (v1 = h1+1; v1 < 64; v1=v1+1)
{
gc2 = _Genetic_Code[v1];
if (gc2 == 10)
{
vShift = vShift + 1;
}
else
{
if (gc1 == gc2)
{
synM [h1-hShift][v1-vShift] = simulatedCodonEFV[h1-hShift];
synM [v1-vShift][h1-hShift] = simulatedCodonEFV[v1-vShift];
}
else
{
nonSynM [h1-hShift][v1-vShift] = simulatedCodonEFV[h1-hShift];
nonSynM [v1-vShift][h1-hShift] = simulatedCodonEFV[v1-vShift];
}
}
}
}
}
2). For every replicate, after it has been fitted, and assuming the tree used for the replicate likelihood function is 'simulatedTree', compute a pair of associative arrays (indexed by branch name) which will be populated with expected number of syn and non-syn subs (per codon here; divide by 3 if you want it by nucleotide):
Code:synSubsAVL = {};
nsSubsAVL = {};
for (h1=0; h1 < T-1; h1=h1+1)
{
abn = branchNames[h1];
ExecuteCommands("GetInformation (aRateMx, givenTree."+abn+");");
synSubs = (horOnes*(aRateMx$synM))*vertOnes;
nsynSubs = (horOnes*(aRateMx$nonSynM))*vertOnes;
synSubsAVL[abn] = synSubs;
nsSubsAVL [abn] = nsynSubs;
}
Then you can do various things with those replicates: write them to file, compute means, etc.
HTH,
Sergei
Home
Pages: 1
dNdS trees - Bootstrap - Var estimates - UI (Read 4877 times)
