eliot
YaBB Newbies
Offline
I love YaBB 1G - SP1!
Posts: 9
|
hi,
You guys have built some very nice tools, thanks.
I have a particular problem I'm writing a batch file for, and I'd like to check if I'm doing what I intend to be doing.
I have two corresponding datasets from noncoding sequence. One is putatively neutral, and the other putatively functional. (for example they could be an ancestral repeat, and a putative regulatory element which are nearby each other). I want to test if the substitution rate for a particular branch (branch A) in the functional is higher than expected relative to the neutral. I define "expected" based on information external to this data--say the ratio of subsitution rates genome wide.
For a 3 species tree ((A,B)C), I put each type of sequence in its own partition, and set up a joint likelyhood function.
LikelihoodFunction theLnLik = (neutFilt, neutTree, consFilt, consTree);
Then to do a likelyhood ratio test, I get the likelyhood under a constrained model. (I'm using HKY85 by the way).
global ss=0; consTree.A.trst:=ss*0.26*neutTree.A.trst; consTree.A.trvs:=ss*0.24*neutTree.A.trvs; consTree.B.trst:=ss*0.68*neutTree.B.trst; consTree.B.trvs:=ss*0.69*neutTree.B.trvs; consTree.C.trst:=ss*0.36*neutTree.C.trst; consTree.C.trvs:=ss*0.35*neutTree.C.trvs;
In the second statement, the number 0.26 is a branch specific scaling factor (obtained externally to this data) giving the average ratio between neutral and conserved for branch A transitions. ss is a scaling factor for the level of constraint of this particular element in all branches of the tree.
I then compare the likelihood for this model with that for a model with one additional parameter. This parameter allows branch A to go more quickly.
global ss=0; global RR=1; RR>1; consTree.A.trst:=ss*RR*0.26*neutTree.A.trst; consTree.A.trvs:=ss*RR*0.24*neutTree.A.trvs; consTree.B.trst:=ss*0.68*neutTree.B.trst; consTree.B.trvs:=ss*0.69*neutTree.B.trvs; consTree.C.trst:=ss*0.36*neutTree.C.trst; consTree.C.trvs:=ss*0.35*neutTree.C.trvs;
Is there some problem with my setting RR>1 like this? Are there other problems with how I have written it?
What I want to be testing is whether there has been a relaxation of contraint in branch A.
thanks much, Eliot
|