ofedrigo
YaBB Newbies
Offline
I love YaBB 1G - SP1!
Posts: 14
|
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
|