Hello Sergei,

I have a set of noncoding regions and want to look for elements that are accelerated in a particular lineage.  To do this I have taken regions surrounding each element as a reference and want to constrain all element branches to be proportional to the corresponding reference branch under the null model.  Under the alternative model, I want to add another variable to the lineage I am interested in that allows the element branch to be longer.  

I am able to constrain two branches to be proportional and I add a constant of proportionality that I have calculated from the global data set.  It looks something like this:  Element_tree.branchA.t:={constant}*Reference_tree.branchA.t

My questions are as follows:
1)      Under the alternative model, how do I estimate (or allow for) a variable, say G, that is > 1 and allows the element branch to be longer?  The formula would look something like this:  Element_tree.branchA.t:=G*{constant}*Reference_tree.branchA.t

2)      In publications doing similar analyses (I’m thinking about Bush&Lahn 2008), an additional scaling factor, S, that is a tree wide scaling factor which accounts for the fact that elements vary in their general level of conservation was estimated.  The formula for constrained branches would be as follows with this factor: Element_tree.branchA.t:=S*{constant}*Reference_tree.branchA.t.  How do I add such a factor into the analysis?

3)      Is there a way to weight the element and reference branches differently in the likelihood?  


Thanks very much in advance for your help.  

Rachel

Thanks so much for the fast response.  I will be running the analysis in batch mode on many alignments.

Hello Sergei,

I've formalized the list of inputs and outputs and I am hoping if you can tell me whether or not I am on the right track or if you can see any problems with my setup.  In short, I want to determine if the species1 branch is evolving at an accelerated rate.

Below is a snip from a test batch file.  The reference partition contains putatively neutral sequence, the element partition contains putatively functional noncoding sequence.   I define the expected rate of evolution through calculating a constant of proportionality external to the data.  S is a tree-wide scaling factor that accounts for varying levels of conservation and G is a variable that allows the species1 terminal branch to be longer under the alternative hypothesis.  

/*Create the partitions*/
DataSetFilter ReferenceFilter = CreateFilter(myData, 1, "0-429, 1119-1718");
DataSetFilter ElementFilter = CreateFilter(myData, 1, "430-1118");

/*Create vector of observed frequencies*/
HarvestFrequencies(obsFreqsRef, ReferenceFilter,1,1,1);
HarvestFrequencies(obsFreqsEle, ElementFilter, 1,1,1);

/*Define models and trees*/
HKY85RateMatrix = {{*,b,a,b}{b,*,b,a}{a,b,*,b}{b,a,b,*}};

Model HKY85Ref = (HKY85RateMatrix, obsFreqsRef);
Tree myTreeRef = ((species1,species2),species3);

Model HKY85Ele = (HKY85RateMatrix, obsFreqsEle);
Tree myTreeEle = ((species1,species2),species3);

LikelihoodFunction theLnLik = (ReferenceFilter, myTreeRef, ElementFilter, myTreeEle);

/* Initialize S, the tree-wide scaling factor*/
global S;

/*Set the constraints*/
myTreeEle.species1.a:=S*0.4*myTreeRef.species1.a;
myTreeEle.species1.b:=S*0.35*myTreeRef.species1.b;
myTreeEle.species2.a:=S*0.2*myTreeRef.species2.a;
myTreeEle.species2.b:=S*0.23*myTreeRef.species2.b;
myTreeEle.species3.a:=S*0.33*myTreeRef.species3.a;
myTreeEle.species3.b:=S*0.35*myTreeRef.species3.b;

Optimize (paramValues, theLnLik);
nullLnLik = paramValues[1][0];

/*Initialize G, the variable that allows the species1 terminal branch to be longer*/
global G; G:>1;
myTreeEle.species1.a:=G*S*0.4*myTreeRef.species1.a;
myTreeEle.species1.b:=G*S*0.35*myTreeRef.species1.b;
myTreeEle.species2.a:=S*0.2*myTreeRef.species2.a;
myTreeEle.species2.b:=S*0.23*myTreeRef.species2.b;
myTreeEle.species3.a:=S*0.33*myTreeRef.species3.a;
myTreeEle.species3.b:=S*0.35*myTreeRef.species3.b;

Optimize (paramValues, theLnLik);

lnlikDelta = 2(nullLnLik-paramValues[1][0]);


Many thanks.

Rachel

Hello Sergei,

I've set up some test cases and am almost ready to run the analysis.  My current plan is to write a perl wrapper to run the analysis for each alignment rather than use the batch analysis scaffold because it is very easy for me to do so.  Is there an advantage to using the hyphy batch scaffold or will I do just as well with my wrapper?  

The putatively neutral reference partitions are substantially larger than the element partitions.  I am wondering if it is possible and/or necessary to weight the element and reference bases differently.  

The constants (0.4, 0.35 etc) in the constraints (called constants of proportionality)  are based on the whole genome ratio between elements and neutral reference for a particular branch.  These are meant to account for lineage specific differences in constraint (which are genome-wide).  I have not yet determined if adding these is the best course of action.    

Thanks again for your help.

Cheers,

Rachel

Sergei wrote on Mar 9^th, 2009 at 11:51am:


Thanks for the link.  Should be quite helpful.  

What about weighting the element and reference bases differently due the differences in length.  Is that possible/necessary?  

Cheers,

Rachel

Just wanted to say that I think I've got everything worked out.  Thanks so much for your help. 

Rachel