Hi,

I am trying to simulate positive selection for non coding regions on a specific branch. For this, I use a four category model (as described in Zhang et al. 2005). It does not seem to work. Basically that's how I try to do it:

I simulate independently 4 dataset for each of the four site class with fraction of sites f0,f1,f2a,f2b respectively and I concatenate them. Here is the relevant part of my code:


ACCEPT_BRANCH_LENGTHS = 1;
characters = {{"A","C","G","T"}
                      {"1","","",""}};

/*set up parameters*/
treeString="(A:0.0070826,B:0.0117252,C:0.0623285)";
global kappa_inv=0.22;
obsFreqs={{0.5}{0.25}{0.25}{0.5}};
branch_name="A";
zeta0=0;
zeta1=1;
zeta2=1.5;
f0=500;
f1=300;
f2a=200;
f2b=250;

/*define model*/
HKY85RateMatrix = {{*,t*kappa_inv,t,t*kappa_inv}
                  {t*kappa_inv,*,t*kappa_inv,t}
                  {t,t*kappa_inv,*,t*kappa_inv}
                  {t*kappa_inv,t,t*kappa_inv,*}};
Model HKY85       = (HKY85RateMatrix, obsFreqs);

useModel(HKY85);

/*create neutral tree"
Tree dummy = treeString;
branchNames     = BranchName   (dummy,-1);
branchLengthsNow  = BranchLength (dummy,-1);
tToBranchLengthFactor:= 2*((obsFreqs[0]*obsFreqs[2]+obsFreqs[1]*obsFreqs[3]) + kappa_inv*(obsFreqs[0]*obsFreqs[1]+obsFreqs[0]*obsFreqs[3]+obsFreqs[1]*obsFreqs[2]+obsFreqs[2]*obsFreqs[3]));
Tree neutralTree=treeString;
for (k = 0; k < Columns (branchNames)-1; k=k+1)  {ExecuteCommands ("neutralTree."+branchNames[k]+".t:="+branchLengthsNow[k]+"/tToBranchLengthFactor;");}

Tree myTree=treeString;

/* zeta0 on foreground and background branches */
ClearConstraints (myTree);
ReplicateConstraint ("this1.?.t:=zeta0*this2.?.t",myTree,neutralTree);
/*simulate dataset*/
DataSet simulatedData = Simulate (myTree,obsFreqs,characters,f0,0);

/* zeta1 on foreground and background branches */
ClearConstraints (myTree);
ReplicateConstraint ("this1.?.t:=zeta1*this2.?.t",myTree,neutralTree);
/*simulate and concatenate dataset*/
DataSet sim = Simulate (myTree,obsFreqs,characters,f1,0);
DataSet simulatedData = Concatenate (simulatedData, sim);

/* zeta0 on background branches and zeta2 on foreground branch*/
ClearConstraints (myTree);
ReplicateConstraint("this1.?.t :=zeta0*this2.?.t", myTree,neutralTree);
ExecuteCommands("newTree."+branch_name+".t := zeta2*neutralTree."+branch_name+".t;");
/*simulate and concatenate dataset*/
DataSet sim = Simulate (myTree,obsFreqs,characters,f2a,0);
DataSet simulatedData = Concatenate (simulatedData, sim);

/* zeta1 on background branches and zeta2 on foreground branch*/
ClearConstraints (myTree);
ReplicateConstraint("this1.?.t :=zeta1*this2.?.t", myTree,neutralTree);
ExecuteCommands("newTree."+branch_name+".t := zeta2*neutralTree."+branch_name+".t;"); 
/*simulate and concatenate dataset*/
DataSet sim = Simulate (myTree,obsFreqs,characters,f2b,0);
DataSet simulatedData = Concatenate (simulatedData, sim);       

/*save concatenated dataset*/
DataSetFilter      simulatedDataFilter = CreateFilter (simulatedData,1);
DATA_FILE_PRINT_FORMAT = 2;
fprintf ("myoutput.phy", CLEAR_FILE,simulatedDataFilter);


Is it the right way to do it?

Thanks a lot,

Olivier

Thanks Sergei!

I am still confused about the branch length scaling. My understanding was that, for instance under HKY model, branch lengths are defined are function of the frequency and  trst trsv  rates:

BranchLength (trvs, trst) = 2*trst*(f_A*f_G+f_C*f_T) + 2*trsv (f_A*f_C+f_A*f_T+f_C*f_G+f_G*f_T)

In your code, in order to calculate the scaling factor you do not account for base freq and kappa_inv. Am I missing something?

thanks again,

Olivier